Solid and Structural Mechanics Demo 3

Solid and Structural Mechanics Demo 3: A unilateral contact problem

This code solves an elastostatic problem involving unilateral contact.

The code essentially looks like lelas2d. The only difference is that since the problem is nonlinear, due to the contact condition, an iterative procedure must be supplied. We describe here the part of the code that performs iterations:

int it = 0; while (it < it_max && err > tolerance) { it++; v = 0.; A = 0.; MeshElements(ms) { Elas2DT3 eq(theElement); eq.Deviator(); eq.Dilatation(); eq.BodyRHS(body_f); eq.ElementAssembly(A); eq.ElementAssembly(v); } MeshSides(ms) { Elas2DT3 eq(theSide,u); eq.BoundaryRHS(bound_f); eq.SignoriniContact(bound_f,penal); eq.SideAssembly(A); eq.SideAssembly(v); } A.Prescribe(ms,v,bc); A.Factor(); A.Solve(v); err = Discrepancy(u,v,2,1); if (verbose > 0) cout << "Iteration: " << it << ". Discrepancy: " << err << endl; }

Let us give some comments on the source code here above:

The contact conduction is implemented via a penalty formulation as a boundary traction. For this, we make a loop over sides. The vector bound_f contains the gap at nodes (the distance between boundary nodes and the obstacle. Note that the only concerned sides are those who have a code equal to CONTACT_A or CONTACT_B

An example

Let us run this program with the data presented in the following project file:

<?xml version="1.0" encoding="ISO-8859-1" ?> <OFELI_File> <info> <title></title> <date></date> <author></author> </info> <Project name="beam"> <verbose>1</verbose> <output>1</output> <save>1</save> <plot>1</plot> <init>1</init> <mesh_file>beam.m</mesh_file> <plot_file>beam.d</plot_file> <prescription_file>beam.pr</prescription_file> </Project> </OFELI_File>

The only thing that may differ from other seen project files is that prescription data is given in a separate file: beam.pr. Let us give this one:

<?xml version="1.0" encoding="ISO-8859-1" ?> <OFELI_File> <info> <title></title> <date></date> <author></author> </info> <Prescription> <BoundaryCondition code="1">0.0</BoundaryCondition> <Traction code="1" dof="2">-1000.0</Traction> </Prescription> </OFELI_File>

Clearly, we have prescribed by this file a (Dirichlet) boundary condition that enforces a value of 0.0 for degrees of freedom with code 1. In addition, we have prescribed a traction of value -1000 for the second degrees of freedom (y-traction) having code 1.