Demo: A unilateral contact problem in elastostatics

This code solves an elastostatic problem involving unilateral contact. The contact is taken into account using a penalty method.

  • 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:
    • We make use in this code of some global variables of to monitor iterations. These variables are MaxNbIterations, theIteration, MaxNbIterations, theTolerance and theDiscrepancy with obvious meaning. Note that all these variables have default values. 

         MaxNbIterations = data.getNbIter();
   theTolerance = data.getTolerance();

    • The iteration loop can be drived by the macro IterationLoop which acts like a loop and which updates iteration index and stops iterations when the value of theDiscrepancy is under theTolerance or when theIteration reaches the value of MaxNbIterations. 

         IterationLoop {
      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(v,bc);
      A.solve(v);
                              
      theDiscrepancy = Discrepancy(u,v,2,1);
      if (verbose > 0)
         cout << "Iteration: " << it << ". Discrepancy: " << theDiscrepancy << endl;
   }

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

    • The contact condition 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
    • We have used global variables of to monitor iterations. These variables are theIteration
      

    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.