Contents     Functions         Previous Next     PDF Index

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

 
[Ms2+Cs+K]{q(s)}=[b] {u(s)}
   {y(s)}= [c] {q(s)} 
   [cint]{q(s)}=0 
    (7.1)

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

range([T]N× (NNC))=ker([cint]NS × N)     (7.2)

Solving problem

 
[TTMTs2+TTCTs+TTKT]{qR(s)}=[TTb] {u(s)}
   {y(s)}= [cT] {qR(s)} 
    (7.3)

is then strictly equivalent to solving (7.1).

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

[ca   ce]N× Ne{
qa
qe
}=0  ⇔   [−(−ce−1ca)   I]{
qa
qe
}=[−G  I] {q} = 0     (7.4)

The subspace with DOFs eliminated is spanned by

[T]N× (NNe)= [
I(NNe)× (NNe) 
GNe× (NNe)
] =  [
I 
ce−1ca
]     (7.5)

The problem that verifies constraints can also be written using Lagrange multipliers, which leads to

[
[Z(s)]
[
G 
I 
]  
  [−G  I]
] {
q 
Fc
} =   {
F 
}     (7.6)

The response can be computed using elimination (equation (7.3)) 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

Fc = [[Zea(s)]+Zee(s) [G]] {q} −Fe = [TeT Z(sT]{qa}−TeTF     (7.7)

A common approximation is to ignore viscous and inertia terms in the resultant, that is assume TeT Z(s) TTeT 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

  {qall,global} = [cGL] {qall,local}

Constraints (mpc, rigid, ...) are defined in local coordinates, that is they correspond to

  {qall,local} = [Tlocal] {qmaster,local}

with qmaster,local master DOFs (DOFs in Case.DOF) defined in the local coordinate system and the Case.T corresponding to

  {qall,global}  = [T]{qmaster,local} = [cGL][Tlocal] {qmaster,local}

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

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.

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

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

©1991-2017 by SDTools
Previous Up Next