The aim of this tutorial is to demonstrate the solution of the axisymmetric equations of linear elasticity in cylindrical polar coordinates.
This implementation of the equations and the documentation were developed jointly with Matthew Russell with financial support to Chris Bertram from the Chiari and Syringomyelia Foundation.
Consider a three-dimensional, axisymmetric body (of density , Young's modulus , and Poisson's ratio ), occupying the region whose boundary is . Using cylindrical coordinates , the equations of linear elasticity can be written as
where , is the stress tensor, is the body force and is the displacement field.
Note that, despite the fact that none of the above physical quantities depend on the azimuthal angle , each can have a non-zero component. Also note that variables written with a superscript asterisk are dimensional, and their non-dimensional counterparts will be written without an asterisk. (The coordinate is, by definition, non-dimensional, so it will always be written without an asterisk.)
A boundary traction and boundary displacement are imposed along the boundaries and respectively, where so that
where is the outer unit normal vector.
The constitutive equations relating the stresses to the displacements are
where is the identity tensor and superscript denotes the transpose. In cylindrical coordinates, the matrix representation of the tensor is
and is equal to the trace of this matrix.
We non-dimensionalise the equations, using a problem specific reference length, , and a timescale , and use Young's modulus to non-dimensionalise the body force and the stress tensor:
The non-dimensional form of the axisymmetric linear elasticity equations is then given by
and the non-dimensional parameter
is the ratio of the elastic body's intrinsic timescale, , to the problem-specific timescale, , that we use for time-dependent problems. The boundary conditions are
We must also specify initial conditions:
oomph-lib, the non-dimensional version of the axisymmetric linear elasticity equations (1) combined with the constitutive equations (2) are implemented in the
AxisymmetricLinearElasticityEquations class. This class implements the equations in a way which is general with respect to the specific element geometry. To obtain a fully functioning element class, we must combine the equations class with a specific geometric element class, as discussed in the (Not-So-)Quick Guide. For example, we will combine the
AxisymmetricLinearElasticityEquations class with
QElement<2,3> elements, which are 9-node quadrilateral elements, in our example problem. As usual, the mapping between local and global (Eulerian) coordinates within an element is given by
where ; is the number of nodes in the element. is the -th global (Eulerian) coordinate of the -th node in the element and the are the element's shape functions, which are specific to each type of geometric element.
We store the three components of the displacement vector as nodal data in the order and use the shape functions to interpolate the displacements as
where is the -th displacement component at the -th node in the element, i.e., .
The solution of time dependent problems requires the specification of a
TimeStepper that is capable of approximating second time derivatives. In the example problem below we use the Newmark timestepper.
As a test problem we consider forced oscillations of the circular cylinder shown in the sketch below:
It is difficult to find non-trivial exact solutions of the governing equations (1), (2), so we manufacture a time-harmonic solution:
This is an exact solution if we set the body force to
where and are the nondimensional Lamé parameters. We impose the displacement along the boundaries according to (4), and impose the traction
along the boundary .
The problem we are solving then consists of equations (1), (2) along with the body force, (5) and boundary traction (6). The initial conditions for the problem are the exact displacement, velocity (and acceleration; see below) according to the solution (4).
The animation below shows the time dependent deformation of the cylinder in the - plane, while the colour contours indicate the azimuthal displacement component. The animation is for , since the time scale is nondimensionalised on the reciprocal of the angular frequency of the oscillations.
As usual, we define all non-dimensional parameters in a namespace. In this namespace, we also define the body force, the traction to be applied on the boundary , and the exact solution. Note that, consistent with the enumeration of the unknowns, discussed above, the order of the components in the functions that specify the body force and the surface traction is .
In addition, the namespace includes the necessary machinery for providing the time dependent equations with their initial data from the exact solution. There are 9 functions, one for each of the components of displacement, velocity and acceleration, and a helper function. For brevity, we list only one of these functions; the others are similar.
We start by creating a
DocInfo object that will be used to output the solution, and then build the problem.
Next we set up a timestepper and assign the initial conditions.
We calculate the number of timesteps to perform - if we are validating, just do small number of timesteps; otherwise do a full period the time-harmonic oscillation.
Finally we perform a time dependent simulation.
The problem class is very simple, and similarly to other problems with Neumann conditions, there are separate meshes for the "bulk" elements and the "face" elements that apply the traction boundary conditions. The function
assign_traction_elements() attaches the traction elements to the appropriate bulk elements.
The problem constructor creates the mesh objects (which in turn create the elements), pins the appropriate boundary nodes and assigns the boundary conditions according to the functions defined in the
Then the physical parameters are set for each element in the bulk mesh.
We then loop over the traction elements and set the applied traction.
Finally, we add the two meshes as sub-meshes, build a global mesh from these and assign the equation numbers.
We create the face elements that apply the traction to the boundary .
The time integration in this problem is performed using the Newmark scheme which, in addition to the standard initial conditions (3), requires an initial value for the acceleration. Since we will be solving a test case in which the exact solution is known, we can use the exact solution to provide the complete set of initial data required. For the details of the Newmark scheme, see the tutorial on the linear wave equation.
If we're doing an impulsive start, set the displacement, velocity and acceleration to zero, and fill in the time history to be consistent with this.
If we are not doing an impulsive start, we must provide the timestepper with time history values for the displacement, velocity and acceleration. Each component of the these vectors is represented by a function pointer, and in this case, the function pointers return values based on the exact solution.
Then we loop over all nodes in the bulk mesh and set the initial data values from the exact solution.
This member function documents the computed solution to file and calculates the error between the computed solution and the exact solution.
AxisymmetricLinearElasticity::youngs_modulus_pt(). The explanation for this is that this function specifies the ratio of the material's actual Young's modulus to the Young's modulus used in the non-dimensionalisation of the equations. The capability to specify such ratios is important in problems where the elastic body is made of multiple materials with different constitutive properties. If the body is made of a single, homogeneous material, the specification of the non-dimensional Young's modulus is not required – it defaults to 1.0.
AxisymmetricLinearElasticityProblem::set_initial_conditionsfunction and compare the system's evolution to that obtained when a "smooth" start from the exact solution is performed.
AxisymmetricLinearElasticityProblem::set_initial_conditionssets the correct initial values for the displacement, velocity and acceleration. (Hint: the relevant code is already contained in the driver code, but was omitted in the code listings shown above.)
A pdf version of this document is available.