PDF Index|
Contents
Functions
PDF Index |
rigid links, FixDof, MPC entries, symmetry conditions, continuity constraints in CMS applications, ... all lead to problems of the form
| (7.2) |
The linear constraints [cint]{q(s)}=0 can be integrated into the problem using Lagrange multipliers or constraint elimination. Elimination is done by building a basis T for the kernel of the constraint equations, that is such that
| (7.3) |
Solving problem
| (7.4) |
is then strictly equivalent to solving (7.2).
The basis T is generated using [Case,NNode,model.DOF]=fe_case(model,'gett') where Case.T gives the T basis and Case.DOF describes the active or master DOFs (associated with the columns of T), while model.DOF or the Case.mDOF field when it exists, describe the full list of DOFs.
The NoT command option controls the need to return matrices, loads, ... in the full of unconstrained DOFs [M], {b} ... or constrained TTMT, TTb in fe_mknl, fe_load, ... .
For the two bay truss example, can be written as follows :
model = femesh('test 2bay'); model2=fe_case(model, ... % defines a new case 'FixDof','2-D motion',[.03 .04 .05]', ... % 2-D motion 'FixDof','Clamp edge',[1 2]'); % clamp edge Case=fe_case('gett',model) % Notice the size of T and fe_c(Case.DOF) % display the list of active DOFs model = fe_mknl(model) % Now reassemble unconstrained matrices and verify the equality % of projected matrices [m,k,mdof]=fe_mknl(model,'NoT'); norm(full(Case.T'*m*Case.T-model.K{1})) norm(full(Case.T'*k*Case.T-model.K{2}))
To compute resultants associated with constraint forces further details are needed. One separates active DOF qa which will be kept and slave DOF that will be eliminated qe so that the constraint is given by
| (7.5) |
The subspace with DOFs eliminated is spanned by
| (7.6) |
The problem that verifies constraints can also be written using Lagrange multipliers, which leads to
| (7.7) |
The response can be computed using elimination (equation (7.4)) and the forces needed to verify the constraints (resultant forces) can be assumed to be point forces associated with the eliminated DOF qe which leads to
| (7.8) |
A common approximation is to ignore viscous and inertia terms in the resultant, that is assume TeT Z(s) T≈ TeT K T.
In the presence of local coordinate systems (non zero value of DID in node column 3), the Case.cGL matrix built during the gett command, gives a local to global coordinate transformation
| (7.9) |
Constraints (mpc, rigid, ...) are defined in local coordinates, that is they correspond to
| (7.10) |
with qmaster,local master DOFs (DOFs in Case.DOF) defined in the local coordinate system and the Case.T corresponding to
| (7.11) |
As a result, model matrices before constraint elimination (with NoT) are expected to be defined in the global response system, while the projected matrix TTMT are defined in local coordinates.
celas use local coordinate information for their definition. cbush are defined in global coordinates but allow definition of orientation through the element CID.
An example of rigid links in local coordinates can be found in se_gimbal('ScriptCgl').
For a DofSet entry, one defines the enforced motion in Case.TIn and associated DOFs in Case.DofIn. The DOFs specified in Case.DofIn are then fixed in Case.T.
Whatever the constraint formulation it requires a transformation into an explicit multiple point constraint during the resolution. This transformation is accessible for RBE3 and rigidconstraints, a cleaned resolution of MPC constraints is also accessible using fe_mpc.
The output is of the format struct with fields
Such format allows the user to transform a constraint into a penalization using the constraint matrix as an observation matrix. One can indeed introduce for each constraint equation a force penalizing its violation through a coefficient kc so that {f}penal = kc [c]Nc × N {q}N × 1. This can be written by means of a symmetric stiffness matrix [kpenal]N × N = kc [c]T [I]Nc × Nc [c]Nc × N added to the system stiffness.
% Transformation of a constraint into a penalty % Generation of a screw model example model=demosdt('demoscrew layer 1 40 20 3 3 space .2 layer 2 40 20 4'); % Model a screw connection with a RBE3 constraint % see sdtweb fe_case.html#ConnectionScrew r1=struct('Origin',[20 10 0],'axis',[0 0 1],'radius',3, ... 'planes',[0 0 111 1 0;3 0 111 1 0; % [z0 type ProId zTol rTol] 5.2 0 112 1 6; 7.2 0 112 1 6], ... 'MatProId',[101 101],'rigid',[Inf abs('rigid')],'NewNode',0); r1.planes(:,2)=1; % RBE3 mo2=fe_caseg('ConnectionScrew',model,'screw1',r1); % display the connection in feplot cf=feplot(mo2);fecom('colordatamat -alpha .1'); % Replace RBE3 by a penalized coupling % Get the constraint matrix r1=fe_mpc('rbe3c',mo2,'screw1'); % remove the RBE3 constraint mo2=fe_case(mo2,'reset'); % Generate the penalization stiffness with default kc kc=sdtdef('kcelas'); SE=struct('DOF',r1.DOF,'Opt',[1;1],... 'K',{{feutilb('tkt',r1.c,kc*speye(length(r1.slave)))}}); % Instance the superelement in the model mo2=fesuper('seadd -unique 1 1 screw1',mo2,SE,[1 1]); % Compute the system modes def=fe_eig(cf.mdl,[5 20 1e3]);
A number of low level commands (feutil GetDof, FindNode, ...) and functions fe_c can be used to operate similar manipulations to what fe_case GetT does, but things become rapidly complex. For example
% Low level handling of constraints femesh('reset'); model = femesh('test 2bay'); [m,k,mdof]=fe_mknl(model) i1 = femesh('findnode x==0'); adof1 = fe_c(mdof,i1,'dof',1); % clamp edge adof2 = fe_c(mdof,[.03 .04 .05]','dof',1); % 2-D motion adof = fe_c(mdof,[adof1;adof2],'dof',2); ind = fe_c(model.DOF,adof,'ind'); mdof=mdof(ind); tmt=m(ind,ind); tkt=k(ind,ind);
Handling multiple point constraints (rigid links, ...) really requires to build a basis T for the constraint kernel. For rigid links the obsolete rigid function supports some constraint handling. The following illustrates restitution of a constrained solution on all DOFs
% Example of a plate with a rigid edge model=femesh('testquad4 divide 10 10');femesh(model) % select the rigid edge and set its properties femesh(';selelt group1 & seledge & innode {x==0};addsel'); femesh('setgroup2 name rigid'); FEelt(femesh('findelt group2'),3)=123456; FEelt(femesh('findelt group2'),4)=0; model=femesh; % Assemble model.DOF=feutil('getdof',model);% full list of DOFs [tmt,tkt,mdof] = fe_mknl(model); % assemble constrained matrices Case=fe_case(model,'gett'); % Obtain the transformation matrix [md1,f1]=fe_eig(tmt,tkt,[5 10 1e3]); % compute modes on master DOF def=struct('def',Case.T*md1,'DOF',model.DOF) % display on all DOFs feplot(model,def); fecom(';view3;ch7')
