Constraint and fixed boundary condition handling

Contents Functions Index

PDF

7.13 Constraint and fixed boundary condition handling

rigid links, FixDof and KeepDOF entries, symmetry conditions, continuity constraints in CMS applications, ... all lead to problems of the form

[Ms²+Cs+K]{q(s)}=[b] {u(s)}

{y(s)}= [c] {q(s)}

[c_int]{q(s)}=0

(7.1)

The linear constraints [c_int]{q(s)}=0 can be integrated into the problem using Lagrange multipliers but the preferred approach here is to eliminate these constraints. This is done by building a basis T for the kernel of the constraint equations

range([T]_{N× (N-NC)})=ker([c_int]_{NS × N}) (7.2)

and solving problem

[T^TMTs²+T^TCTs+T^TKT]{q_R(s)}=[T^Tb] {u(s)}

{y(s)}= [cT] {q_R(s)}

which is strictly equivalent to solving (7.1).

The basis T is generated using [Case,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 tt model.DOF describes the full list of DOFs.

The assembly of unconstrained M, ... or constrained T^TMT matrices can be controlled with appropriate options in fe_mk, fe_load, ...

When defining local displacement bases (non zero value of DID in node column 3), master DOFs are defined in the local coordinate system. As a result, M is expected to be define in the global response system while the projected matrix T^TMT is defined in local coordinates. mpc constraints are defined using the local basis.

For the two bay truss example, can be written as follows :

 femesh('reset');
 model2 = femesh('test 2bay');
 model2=fe_case(model,'SetCase1', ...         % 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',model2)  % Notice the size of T and 
 fe_c(Case.DOF)                % display the list of active DOFs
 model2 = fe_mknl(model2)
 
 % Now reassemble unconstrained matrices and verify the equality
 % of projected matrices
 [m,k,mdof]=fe_mk(model2,'options',[0 2 2]);

 norm(full(Case.T'*m*Case.T-model2.K{1}))
 norm(full(Case.T'*k*Case.T-model2.K{2}))

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

 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 this is supported by the rigid function. The following illustrates restitution of a constrained solution on all DOFs

% Example of a plate with a rigid edge
femesh('reset');
femesh('testquad4;divide 10 10;addsel');
% 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;
model.pl=m_elastic('dbval 100 steel');
model.il=p_shell('dbval 110 Mindlin 5e-2');

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