Loading

Solution of an Elastodynamics problem

We present a program to solve the time dependent linear elasticity
equations (*i.e.* Elastodynamics) using
the P_{1} finite element method for space discretization and
the Newmark scheme for time discretization.

- We start, as usual, by including the main header
file of the library, and the header for Solid Mechanics problems.
The program has, as argument the mesh file and the time step.
We construct a Mesh
instance with the flag that enables eliminating imposed boundary conditions.
Note that the final time value is also set by using the OFELI global
variable theFinalTime. Next, we declare the
vectors to use in the calculations: u and v
will contain the initial conditions for displacement and velocity, and
u will contain at each time step the
corresponding displacement vector. The initial displacement is set to
the constant vector (0,-1)
^{T}, and the initial velocity to zero by default.#include "OFELI.h" #include "Solid.h" using namespace OFELI; int main(int argc, char *argv[]) { theFinalTime = 1.; Mesh ms(argv[1],true); theTimeStep = atof(argv[2]); Vect<double> u(ms), v(ms), bc(ms); u.set("-0.1",2);

- Now we instantiate the equation class
Elas2DT3 (Linearized
planar elasticity using the P
_{1}triangle). We also say that we the used terms in the equation are lumped mass, stiffness term (deviatoric and dialational). No body force is applied. We transmit all these data to the TimeStepping instance and choose as linear solver the conjugate gradient, preconditioned with the diagonal ILU.

Since we intend to store for plotting displacement at all time steps, we store the initial displacements in a vtk file.Elas2DT3 eq(ms); eq.setTerms(LUMPED_MASS|DEVIATORIC|DILATATION); TimeStepping ts(NEWMARK,theTimeStep,theFinalTime); ts.setPDE(eq); ts.setInitial(u,v); ts.setLinearSolver(CG_SOLVER,DILU_PREC); saveField(u,ms,"u-0.vtk",VTK);

- We can now loop over time steps. For each one, we
set boundary conditions (Here 0). We run the
time step and store the resulted displacement vector in a vtk file.
We can also store a deformed mesh in another file (More used in Solid
Mechanics).
TimeLoop { bc.setTime(theTime); ts.setBC(bc); ts.runOneTimeStep(); saveField(u,ms,"u-"+itos(theStep)+".vtk",VTK); Mesh dm(ms); DeformMesh(dm,u,1.); dm.save("mm-"+itos(theStep)+".m"); } return 0; }

Copyright © 1998-2016 Rachid Touzani