PDF Index| SDT-nlsim
Contents
Functions
PDF Index |
Purpose
Element property function for contact. This function is used by SDTools for internal developments and is not distributed with the base SDT. Integrated calls for high level handling are available using ctc_utils, this section aims at documenting the implementation details.
Syntax
il=p_contact('dbval 1')
Description
For theory see section 3.1, for implementation details, see section 3.2. Contact is defined in SDT using contact elements. For 3D volumes, contact elements are then 3D shell elements. Since contact is non-linear, these elements are used in combination with the function nl_contact for computation of the contact force. See nl_spring, nl_solve to find out about non-linearities handling in the SDT/non-linear vibrations toolbox.
p_contact handles the contact formulation for integration, the il data is a structure that must be stored as a pro entry in the model Stack to comply with usual handling of non linear springs in SDT (through nl_spring).
Stack entry pro,Contact then follows the format
il = struct('type','p_contact',.... 'il', il,... 'NLdata',NLdata,... 'MAP',Map));
fields il, MAP, and NLdata are defined below.
il = [ProID fe_mat('p_contact','SI',1) ... Match Integ ContactLaw SlaveId FrictionLaw Euler TangStick LocalFric1stDir TOIM]
Physical parameters are stored in p_contact NLdata. The formulation options are as follows:
| ProID | Element property identification number |
| type | Identifier obtained with fe_mat('p_contact','SI',1) |
| Match | Matching strategy (tolerance and choice of contact points). Format is i.val. |
| i (integer part) select the matching strategy. 0 (default): Gauss integration point and matching to the volume underlying the slave surface, 1: Gauss point matching to the slave surface face mesh (quicker), -2 matches using the element nodes as contact point (warning: relevant only when Integ is set at -2000 meaning nodal pressures). | |
| val (digits after the comma) the tolerance token. Set to 0 uses the default matching radius of 1. With x.1 lets p_contact define the matching radius based on the slave surface mesh characteristic length (e.g set Match to 1.1 to use Gauss contact points matched to the surface with a radius based on the mesh). To customize the radius to r define val as val= 10*r+9. Thus Match=0.209 to use Gauss contact points matched to the slave surface underlying volume with a radius of 20. For radii less than 1 the same approach stands, as 1.029 with be used to define a radius of 0.2. | |
| Integ | Contact elements integration rule (see integrules for selection). It will determine the Gauss points used for contact. Adding 1000 to the integrule identifier allows activating a CPress operator converting contact resultants at node to contact pressures at nodes or contact points. This is currently used by nl_contactl. Adding .1 to Integ will generate the observation based on the deformed mesh according to a deformation curved stored in the model stack as curve,initContact. |
| ContactLaw | contact law selection : 1 bilinear, 2 exponential pressure, 3 tabular, 4 power law, 9 external contact law function (user defined). Lagrange contact is defined in nl_contactl. |
| SlaveId | Identification of the set specifying the slave surface, by Id. The slave set data must then contain and ID field with a relevant value. For generation, see ctc_utils Generate or feutil AddSet. |
| FrictionLaw | Friction law selection : 0 none, 1 Coulomb unidirectional, 2 Regularized Coulomb unidirectional with a penalty slope, 3 Regularized Coulomb unidirectional with an arctangent function, 9 external friction law function (user defined). These are all 1D friction laws, to have their 2D equivalent, add 10 to the number: 11, 12, 13 and 19 are supported. |
| Euler | Euler formulation for the sliding velocity computation, to be used with friction. This adds a sliding velocity to the observation on the velocity field. This is useful to simulate displacements on fixed meshes. Use 1 for rotating bodies (this will add the Euler deformation of the meshes as well) or 2 for translating bodies. The Euler formulation assumes sliding is happening in the first friction direction, see MAP definitions below and p_shell. |
| TangStick | Sticking tangent state, for mode computation. If contact is stuck, you may want to add a symmetric stiffness coupling representing the impossibility for the system to slide in this direction for small vibrations. 0 none, 1 used. |
| LocalFric1stDir | Ordering of local friction directions. For each contact point, a local basis is defined, with the 3rd direction as the normal. This option allows customizing the ordering of the local friction direction. If set to 2 the friction directions are permuted and the new first friction direction orientation is made the same regardless of the normal direction. Is is possible to reverse the orientation by setting this value to -2. Set this to 0 to deactivate it. |
| TOIM | Topology Offset Induced Momentum, default to 0. If a clearance exists in the mesh topology between contact surface candidates, one can choose whether this has to be accounted for to define relative movement observation. If this is considered the observation will consider the clearance as a rigid section thus transmitting moments between the spaced contacting points. This feature is mostly useful when considering contact between shells, where one must account for the shell thickness. Set 0 to ignore (the default for solids), or 1 for use. |
The contact normal is defined by the shell element normal. By default, friction directions are defined element-wisely depending on their orientation. If you want to control these directions, or even use customized local frames, a map needs to be defined, and added to the il data structure. These orientation maps are detailed in p_shell(of_mk('BuildNDN') integration rule 23).
For p_contact, the map third direction v3, (v3x,v3y,v3z) is fixed in the contact normal direction (of the element). The first two directions v1 and v2 concern friction directions and can be modified. It typically suffices to define the first direction, the second one being completed to form a direct orthogonal frame. The first vector v1x,v1y,v1z is built in the direction of r lines (element edge defined) and v2x,v2y,v2z is tangent to the surface and orthogonal to v1.
For example, to fix the first sliding direction in the y direction of the global frame, and keep the third direction as the local normal, the map can be defined as
MAP=struct('dir',{{0,1,0,'v3x','v3y','v3z'}},... 'lab',{{'v1x','v1y','v1z','v3x','v3y','v3z'}})
To define the first sliding direction at 45 degrees of the global frames x,z direction,
MAP=struct('dir',{{1/sqrt(2),1/sqrt(2),0,'v3x','v3y','v3z'}},... 'lab',{{'v1x','v1y','v1z','v3x','v3y','v3z'}})
It is possible to define a field Orig in the format [x0 y0 z0] in the local contact MAP for use with the definition of rotation in the fixed frame Euler=1.
The contact and friction laws (see section 3.3 for details) are handled by the function nl_contact. This is a specific non linear function, see nllist, and nl_fun for more details. The NLdata is a structure with fields
NLdata=struct('type','nl_contact',... 'Fu','Kc 1e12',... 'Fv','Mu 0 KtLin 50',... 'cqoff',0,... 'Vel', 0,... 'KcLin', 1, ... 'Sel' , 0 ,... );
| type | fixed value defining the non linearity type to be nl_contact |
| Fu | the contact law parameters, in a string format (defined below) |
| Fv | the friction law parameters, in a string format (defined below) |
| cqoff | offset for gap computations. This is a scalar added to the gap vector. To allow taking into account topology offsets (e.g mesh showing open contact), one can set cqoff either to +Inf to generate an offset vector based on the gap observation of the mesh with zero deformation, or to -Inf to generate an offset vector based on the positive only gap observation of the mesh with zero deformation. |
| Vel | sliding velocity. This is a scalar added to the sliding velocity vector when using Euler friction formulations. |
| KcLin | fixed Jacobian contact stiffness. This is a scalar defining a linear bilateral contact for the Jacobian computation in time simulations using the NLNewmark scheme of fe_time, to improve convergence. |
| Sel | to output the selection generating the gauss point mesh for contact field observation |
| nodeVel | an optional field to provide a restricted list of NodeId subjected to fixed frame motion (option Euler) if needed. The restriction is then applied to the contact points whose gap observation involves DOF relative to the nodes specified. |
The contact law parameters Fu are defined in a string interpreted by nl_contact depending on the contact law specified in the il field.
The contact laws implemented deal with gap/pressure relationships to conform with the Signorini law.
The string is a succession of names and values representing the reserved parameter name and its associated value. The parameters interpreted depend on the contact law and are defined below. It is also possible to specify the internal laws parameters at low level by specifying the data vector specified.
| 1 linear | Parameter Kc giving the contact stiffness. E.g. NLdata.Fu='Kc 1e12' - or vector [Kc]. |
| 2 exponential law | Parameter pz defines the contact pressure at zero gap and lambda defines the exponential coefficient. E.g. NLdata.Fu='pz 1e2 lambda 75e4' - or vector [pz lambda]. |
| 3 tabular contact law | Parameter CurveId gives the curve ID (stacked in the model) giving the behavior. The curve is of the standard SDT format (sdtweb curve). With field X giving overclosures (opposite of the gap) and field Y giving the corresponding contact pressures. E.g. NLdata.Fu='CurveId 1' - or vector [CurveId]. |
| 4: power contact law | Parameter pz defines the contact pressure at zero gap and parameter power defined the power used on the gap. The local contact pressure is then computed as fN = pz · (−g)power if g < 0 and 0 else (which will induce a positive contact force). The tangent state uses the derivative of the given power function. E.g. NLdata.Fu='pz 1 power 3' - or vector [pz power] . |
| 9: external/user defined | The data string starts with the name of the non linear function called, followed by the parameters defined for the said nl_function, see nl_fun. E.g. NLdata.Fu='nl_mycontact Kc 1e4'. Note that the keywords Fu, Fc, cqoff, KcLin are reserved and must not be used in subfunctions. Parameter interpretation, may fail in such case. |
| 9: Lagrange contact | This is an external non linear function,nl_contactl. This law does not take any parameter, such that the .Fu field need to be set as NLdata.Fu='nl_contactl'. Only a static model (based on the status method) has been implemented, to be solved with nl_solve Uzawa. |
The bilinear law allows controlling the contact surface non-linear behavior for simplified studies. These variations are available with the string input format of the contact law definition. One can then use command option SymForceSetval to define a linearized behavior for transient and mode computation. The impacted elements can be restricted to a list defined by an EltId element set named SymForceSet stacked in the model.
Eventually, using contact law 1000 will use the static state stacked in the model to generate a bilinear contact behavior for points into contact and no contact for opened points. This is mainly used to simplify the behavior in transient simulations, mostly for verification purposes.
The friction law parameters Fv are defined in a string interpreted by nl_contact depending on the friction law specified in the il field. The string is a succession of names and values representing the reserved parameter name and its associated value. The parameters interpreted depend on the friction law and are defined below. It is also possible to specify the internal laws parameters at low level by specifying the data vector specified.
| 0: no friction | the string can be empty (no interpretation) E.g. NLdata.Fv=''. |
| 1, 11: Coulomb law | Parameter Mu gives the friction coefficient, optional parameter Nu gives the tangent sticking coefficient property, see section 3.3.E.g. NLdata.Fv='Mu 0.2' - or vector [Mu Nu]. |
| 2, 12: Regularized Coulomb , linear | Parameter Mu gives the friction coefficient, parameter KtLin gives the slope a zero sliding velocity, optional parameter Nu gives the tangent sticking coefficient property, see section 3.3.E.g. NLdata.Fv='Mu 0.2 KtLin 50 Nu 0.1' - or vector [Mu KtLin Nu]. |
| 3, 13: Regularized Coulomb , arctangent | Parameter Mu gives the friction coefficient, parameter KtLin gives the slope a zero sliding velocity, optional parameter Nu gives the tangent sticking coefficient property, section 3.3. NLdata.Fv='Mu 0.2 KtLin 50 Nu 0.1' - or vector [Mu KtLin Mu]. |
| 4, 14: Regularized Coulomb, scaled linear | Parameter Mu gives the friction coefficient, parameter KtLin gives the slope a zero sliding velocity, optional parameter Nu gives the tangent sticking coefficient property, see section 3.3.E.g. NLdata.Fv='Mu 0.2 KtLin 50 Nu 0.1' - or vector [Mu KtLin Nu]. |
| 9, 19 external/user defined | The data string starts with the name of the non linear function called, followed by the parameters defined for the said nl_function, see nl_fun. E.g. NLdata.Fu='nl_myfriction Mu 0.4'. The same reserved keywords than for external contact law definition apply. |
| 9: Lagrange friction | This is an external non linear function, nl_frictionl. Parameter Mu gives the friction coefficient. The .Fu field needs to be set as NLdata.Fu='nl_frictionl Mu 0.5'. Only a static model (based on the status method) has been implemented, to be solved with nl_solve Uzawa. |
It is possible to obtain default values for NLdata or p_contact data, using
% Standard il field, for ProId 1 il = p_contact('dbval 1 -struct') % Standard il field with MAP defining first sliding direction along global x il = p_contact('dbval 1 tan1dx -struct') % 1D friction il = p_contact('dbval 1 tan2dx -struct') % 2D friction % Standard il field with MAP defining first sliding direction along global y il = p_contact('dbval 1 tan1dy -struct') % 1D friction il = p_contact('dbval 1 tan2dy -struct') % 2D friction % Standard il field with a polar MAP il = p_contact('dbval 1 polar -struct') % Standard il field with a polar MAP set at a custom origin [20 0 0] il = p_contact('dbval 1 polar o 20 0 0 -struct')
Since the standardized default output of dbval commands is a il line vector, command option -struct has been added to allow pro entries under the struct format, to set in the model stack.
Default MAPs are exploitable in the dbval command that can be used by specifying the MAP name in the dbval command. These MAPS are 2D maps in the friction plane,
It is possible to obtain NLdata fields using nl_contact('db'). Without other information, a default field will be output, but it is also possible to give the wanted parameters. In this latter case, please note that all mandatory parameters must be specified.
% default field (ContactLaw=1, FrictionLaw=0); NLdata=nl_contact('db'); % Parameters for a tabular contact law (3) % and regularized Coulomb friction law (12) NLdata=nl_contact(... 'db cqoff 0 Vel 0 KcLin 1e4 Fu''CurveId 18'' Fv''Mu 0.3 KtLin 50 Nu 0.25''');
The il, NLdata and MAP fields can be edited using ctc_utilsset command.
A contact implementation requires the definition of a slave surface and a master surface. For the latter one, contact elements must be defined based on the underlying surface. Integrated calls can used with ctc_utilsGenerateContactPair (see example below). The steps necessary to define contact and the corresponding lower level calls are detailed in this section.
The slave surface must be stored in the model stack, with an explicit ID, as explained in the il definition. It can be generated using feutilb:
data=feutilb('FaceSet',model,EltFaceSel); % generate a FaceId set data.ID=1; model=stack_set(model,'set','slave',data); % stack in model
model is a standard SDT model and EltFaceSel is an element selector string (see sdtweb findelt), that must generate a surface selection. Command option -ID allows directly specifying a set ID.
The generation of the contact elements is handled by nl_contact, with call build:
model=nl_contact('build',model,EltFaceSel);
model is a standard SDT model and EltFaceSel is an element selector string (see sdtweb findelt), that must generate a surface selection. By default the contact elements ProId is incremented from the model. Note that no MatId is associated to contact elements.
Command option ProIdVal allows specifying the contact elements ProId to Val.
The non-linear contact forces are stored in the def.FNL field, see nl_spring. The non-linear contact DOF are defined using defined DOF extensions, .5, 51, .52, .53, .54, .55, respectively corresponding to the contact pressure F_N, the friction pressure in the first friction direction F_T1, the friction pressure in the second friction direction F_T2, the Eulerian deformation Utt (geometric deformation in the first friction direction when Euler sliding is applied), the gap G, the sliding velocity norm W. Not all instances are always present depending on the chosen contact formulation. .5, .53 and .54 are at least always output.
This is a case where the same buffer is used for .unl and .snl. .snl contents are FN, FT1, FT2, Utt. The interlacing is also optimized for MATLAB operations and thus uses g at all Gauss points, the uT1 at all Gauss, ...
The preparation is done in an Init and Exit phase. At the end of the Exit phase, model.DOF=Case.DOF and the original model DOFs are stored in Case.mDOF assemble -fetime -matdes 2 3 1 -se load -initfcn "model=nl_spring("InitTimeProp",model);"
model=demosdt('demogartfe') cf=feplot; % open FEPLOT and define a pointer CF to the figure cf.model=model;
Assembling strategy xxx
Additional fields used for non-linear time integration
model.FNLDOF, mo1.FNLlab should be moved to output generation . Current FNLDOF is included in out.DOF.
sdtweb t_nlspring('of_time') tests move to separate out.FNL which should become the reference.
Generated by nl_fun('init') called by nl_spring('nl -storefnl'). xxx -storefnl should be done by non-linearity.
The following example creates a contact linearity between two cubes, computes a static solution and performs some contact field visualizations
First with an integrated call using ctc_utils
model=d_contact('Cubes -offset0'); % Load model model=ctc_utils('GenerateContactPair',model,... struct('ProId',201,'name','contactPro',... 'master','Group2 & SelFace & InNode{z==1}',... 'slave','Group1 & SelFace & InNode{z==1}',... 'contact','linear','Fu','Kc=1e11',... 'friction',0)); % compute a static solution q0=nl_solve('static tol1e-9',model); % perform some displays cf=feplot(model,q0); % show pressure field ctc_utils('showPn-reset',cf.mdl,q0); % show gap field in mm ctc_utils('showG-unit"MM"inMesh',cf.mdl,q0); % show pressure field in the mesh containing the top cube ctc_utils('showPn inMeshEltSel"withnode{z>1.001}"',cf.mdl,q0)
Then with lower level calls
model=d_contact('cubes -offset0'); % Load model % Generate slave surface (ID 301) model=feutil('AddSetFaceId -ID301',model,'slave','Group1 & SelFace & InNode{z==1}'); % Generate contact elements (master), ProId 201 model=ctc_utils('build proid201',model,'Group2 & SelFace & InNode{z==1}'); % set pro, default model=stack_set(model,'pro','contact',... p_contact('dbval 201 tan1dx slaveid 301 -struct')); % compute a static solution q0=nl_solve('static tol1e-9',model);
A usual example is the simulation of an impact between two solids. The behavior of such model is studied to gauge the integration scheme precision. The modified non linear Newmark scheme set here tackles down the problems commonly encountered with the basic Newmark scheme regarding the precision of the solution.
model=d_contact('cubes -cBuild'); % define the intial condition, with NL DOFs [model,Case]=fe_case(nl_spring('AssembleCall'),model); [i1,r1]=feutil('findnode group2',model); q0=struct('DOF',Case.mDOF,'def',[zeros(size(Case.mDOF)) zeros(size(Case.mDOF))]); q0.def(fe_c(Case.mDOF,[i1+.03],'ind'),2)=-.1; model=stack_set(model,'curve','q0',q0); % Modify the TimeOpt Newmark parameters Opt=nl_solve('timeoptnlnewmark'); Opt.Opt=[.25 .5 0 2e-3 1e2]; Opt.RelTol=-1e-4; model=stack_set(model,'info','TimeOpt',Opt); % launch simulation and post treat contact data def=fe_simul('time',model,model.Stack{end,3}); feplot(model,def); fecom colordataa figure(4); subplot(211);plot(def.data,fe_c(def.DOF,[70.03])*def.def); xlabel('Time [s]'); ylabel('Node 65 disp.') subplot(212);plot(def.data,fe_c(def.FNL.DOF,[.5])*def.FNL.def); xlabel('Time [s]'); ylabel('fn at contact points [Pa]')
See also
nl_solve nllist nl_spring fe_timectc_utils, d_contact
