Ctrl+K

7.14  Constraint and fixed boundary condition handling

7.14  Constraint and fixed boundary condition handling#

7.14.1  Theory and basic example#

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'mCase.T-model.K{1}))
norm(full(Case.T'kCase.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.

7.14.2  Local coordinates#

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').

This built-in implementation may be impractical in practice as this generates a resolution in the local frame and a response in the global frame. It is thus not conforming to another classical formalism where the system is resolved in the global frame and response provided in the local frame.

In practice it is thus not recommended to exploit DID during analysis. It should only be used as intermediate steps in pre/post procedures.

  • Definition of boundary conditions (linear constraints, imposed displacement, external forces) in a local frame. This allows for user-friendly inputs in some advanced applications. In such application define local frames in the concerned nodes DID and define the associated case entries in the local frame. DofLoad, DofSet, MPC, FixDOF, RBE3, rigid, SensDOF entries are supported. Use function feutilb CaseL2G to resolve the case in the global frame before running simulations. This operation removes DID entries and stores the [cGL] matrix in the model stack.
  • Recovery of displacement response in a local basis. The DID implementation does not allow this. You can build the local to global matrix using basis trans command combined with DID definition in a dedicated Node matrix. If feutilb CaseL2G was previously used, the [cGL] matrix is stored in model stack entry curve,OLDcGL. Simply transform the response by using the stored matrix transpose.

The following structural elements support DID definition, rigid, celasand mass2.
rigidelements are treated as a case entry.
celasand mass2elements are projected in the global frame. Command feutilb CaseL2G will thus assemble them into a coupling superelement prior to removing DID entries. Few elements support DID providing elements matrices directly in the global frame, celasand mass2elements.

7.14.3  Enforced displacement#

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.

7.14.4  Resolution as MPC and penalization transformation#

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.

  • RBE3c provides the resolution for RBE3 constraints.
  • RigidC provides the resolution for rigidconstraints.
  • MPCc provides the resolution for MPC constraints.

The output is of the format struct with fields

  • c the constraint matrix.
  • DOF the DOF vector relative to the constraint.
  • slave slave DOF indices in DOF.

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]);

7.14.5  Low level examples#

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')

Contents