This document provides the theoretical background to oomphlib's
solid mechanics capabilities. We start with a review of the relevant theory to establish the overall framework and notation, and then discuss the implementation of the methodology in oomphlib
.
Here is an overview of the structure of this document:
If you're not keen on theory, you may prefer to start by exploring the solid mechanics tutorials in oomphlib's
list of example driver codes.
All problems discussed so far were formulated in an Eulerian frame, i.e. we assumed all unknowns to be functions of the spatiallyfixed, Eulerian coordinates, , and of time, . Many problems in solid mechanics are formulated more conveniently in terms of bodyattached, Lagrangian coordinates. In this section we will briefly review the essential concepts of nonlinear continuum mechanics and present the principle of virtual displacements which forms the basis for all largedisplacement solid mechanics computations in oomphlib
.
Throughout this section we will use the summation convention that repeated indices are to be summed over the range of the three spatial coordinates and all free indices range from 1 to 3. We will retain the summation signs for all other sums, such as sums over the nodes etc.
The figure below introduces the essential geometrical concepts: We employ a set of Lagrangian coordinates, , to parametrise the (Eulerian) position vector to material points in the body's undeformed position:
The specific choice of the Lagrangian coordinates is irrelevant for the subsequent development. For analytical studies, it is advantageous to choose a bodyfitted Lagrangian coordinate system (as shown in the sketch) because this allows boundary conditions to be applied on isolevels of the coordinates; in computational studies, it is usually preferable to use a coordinate system in which the governing equations are as compact as possible. A Cartesian coordinate system is best suited for this purpose.
We denote the tangent vectors to the coordinate lines in the undeformed configuration by
and define the components of the covariant metric tensor via the inner products
This tensor defines the "metric" because the (square of the) infinitesimal length, , of a line element with coordinate increments is given by
The volume of the infinitesimal parallelepiped formed by the coordinate increments is given by
where
As the body deforms, the Lagrangian coordinates remain "attached" to the same material points. The body's deformation can therefore be described by the vector field that specifies the position vectors to material particles in the deformed configuration,
As in the undeformed coordinate system, we form the tangent vectors to the deformed coordinate lines and denote them by
The inner product of these vectors defines the metric tensor in the deformed configuration
and we have equivalent relations for the lengths of line elements,
and the volume of infinitesimal parallelepipeds,
where
Since the metric tensors and provide a measure of the length of material line elements in the deformed and undeformed configurations, respectively, their difference
provides an objective measure of the strain and is known as the Green strain tensor.
Let us now assume that the body is subjected to
and its displacements are prescribed on the remaining part of the boundary, (where and ), i.e.
for given .
The deformation is governed by the principle of virtual displacements
where is the (symmetric) second PiolaKirchhoff stress tensor, is the density of the undeformed body, and represents the variation of . See, e.g., Green, A.E. & Zerna, W. "Theoretical Elasticity", Dover (1992); or Wempner, G. "Mechanics of Solids with Applications to Thin Bodies", Kluwer (1982) for more details.
Upon choosing the particles' position vector in the deformed configuration, , as the unknown, the variation of the strain tensor becomes
and we have
The 2nd Piola Kirchhoff stress tensor is symmetric, therefore we can write the variation of the strain energy in (2) as
and obtain the displacement form of the principle of virtual displacements:
This must be augmented by a constitutive equation which determines the stress as a function of the body's deformation, (and possibly the history of its deformation). Here we will only consider elastic behaviour, where the stress is only a function of the strain.
For purely elastic behaviour, the stress is only a function of the instantaneous, local strain and the constitutive equation has the form
The functional form of the constitutive equation is different for compressible/incompressible/nearincompressible behaviour:
which ensures that the volume of infinitesimal material elements remains constant during the deformation. This condition is typically enforced by a Lagrange multiplier which plays the role of a pressure. In such cases, the stress tensor has the form
where only the deviatoric part of the stress tensor, depends directly on the strain. The pressure must be determined independently by enforcing the incompressibility constraint (4). Given the deformed and undeformed metric tensors, the computation of the stress tensor for an incompressible material therefore requires the computation of the following quantities:
where the deviatoric part of the stress tensor, depends on the strain. This form of the constitutive law is identical to that of the incompressible case and it involves a pressure which must be determined from an additional equation. In the incompressible case, this equation was given by the incompressibility constraint (4). Here, we must augment the constitutive law for the deviatoric stress by an additional equation for the pressure. Generally this takes the form
where is the "bulk modulus", a material property that must be specified by the constitutive law.
which only involves quantities that remain finite as we approach true incompressibility.
The abstract base class ConstitutiveLaw
provides interfaces for the computation of the stress in all three forms.
A hyperelastic material is one for which the stress can be derived from a potential function , known as the strainenergy function, and
A strainenergy function exists if the elastic deformations are reversible and isothermal, or reversible and isentropic. If the material is homogeneous and isotropic then the strainenergy function depends only on the three strain invariants:
and can be written . It may be shown, see Green & Zerna, that
where
The abstract base class StrainEnergyFunction
provides the interfaces and and should be used as the base class for all strainenergy functions. A class StrainEnergyFunctionConstitutiveLaw
that inherits from ConstitutiveLaw
uses a specified strainenergy function to compute the appropriate stresses.
The principle of virtual displacements (2) is written in dimensional form. We generally prefer to work with nondimensional quantities and will now discuss the nondimensionalisation used in the actual implementation of the equations in oomphlib
. For this purpose we first rewrite equation (2) as
where we have used asterisks to label those dimensional quantities that will have nondimensional equivalents. (Some quantities, such as the strain, are already dimensionless, while others, such as the density will not have any nondimensional counterparts. We do not introduce modifiers for these).
We now nondimensionalise all lengths with a problemspecific lengthscale, , given e.g. the length of the solid body, so that
We use a characteristic stiffness, , (e.g. the material's Young's modulus ) to nondimensionalise the stress and the loads as
and we nondimensionalise time with a problemspecific timescale (e.g. the period of some external forcing), so that
This transforms (9) into
where
is the ratio of system's "intrinsic" timescale, to the time, , used in the nondimensionalisation of the equations. If a given problem has no externally imposed timescale (e.g. in the free vibrations of an elastic body) (or some suitable problemdependent multiple thereof) provides a natural timescale for the nondimensionalisation. Therefore we use as the default value in all solid mechanics equations. If preferred, computations can, of course, be performed with dimensional units, provided all quantities are expressed in consistent units (e.g. from the SI system). In this case represents the dimensional density of the material.
We adopt a similar approach for nondimensionalisation of the constitutive equations. Typically, the constitutive parameters (e.g. Young's modulus and Poisson's ratio for a Hookean constitutive equation) are passed to the ConstitutiveLaw
as arguments to its constructor. Where possible, we select one of these parameters as the reference stress and give it a default value of 1.0. Hence, if a Hookean constitutive law is instantiated with just one argument (the Poisson ratio ), the stress is assumed to have been scaled on Young's modulus. If two arguments are provided, the second argument should be interpreted as the ratio of the material's Young's modulus to the reference stress used in the nondimensionalisation of the equations.
Many solid mechanics problems can be regarded as twodimensional in the sense that the quantities of interest only depend on two spatial coordinates. In such problems it is important to consider what constraints the system is subjected to in the third direction — clearly all real solid bodies are threedimensional! In plane strain problems, the displacements of material points are assumed to be parallel to the 2D plane spanned by the coordinates and , so that any displacements normal to this plane are suppressed. In this case, significant transverse stresses can develop. Conversely, in plane stress problems, it is assumed that no stresses develop in the transverse direction; in this case we must allow material particles to be displaced transversely. Since we have formulated our problem in terms of positions (i.e. in terms of the displacement of material points), our formulation naturally produces a plane strain problem if we reduce the equations to two dimensions: We assume that the transverse displacement vanishes, while the remaining (inplane) displacements are only functions of the inplane coordinates, and .
It is important to remember that the 2D version of all equations must produce planestrain behaviour when new, strainenergybased constitutive equation classes are formulated; when implementing such strainenergy functions, recall that any invariants of metric tensors etc. are the invariants of the full 3D quantities which are not necessarily the same as those of the corresponding 2D quantities.
Many biological tissues undergo growth processes. For instance, the cells that make up a solid tumour divide regularly and as a result the total mass of the tumour increases. If the growth occurs nonuniformly so that certain parts of the tumour grow faster than others, regions that grow more slowly restrain the expansion of their neighbours. This process can induce significant growth stresses. The scenario is similar (but not identical) to that of thermal growth in which a nonuniform temperature distribution in a body creates thermal stresses. (The important difference between these two cases is that in the latter process mass is conserved – thermal expansion only leads to an increase in the volume occupied by the material, whereas biological growth via cell division increases the mass of the tumour).
It is easy to incorporate such growth processes into our theoretical framework. The general idea is sketched in the following figures:
The individual infinitesimal material elements have expanded (or contracted) isotropically and the elements are in a stressfree state. The isotropic growth is spatially uniform – all elements have expanded (or contracted) by the same amount.  All infinitesimal material elements have expanded (or contracted) isotropically and the elements are in a stressfree state. The isotropic growth is spatially nonuniform, so individual elements have grown (or contracted) by different amounts. 
Since the isotropicallygrown infinitesimal material elements have grown (or contracted) by the same amount, they can be (re)assembled to form a continuous, stressfree body.  Since the individual material elements have grown (or contracted) by different amounts they can (in general) not be (re)assembled to form a continuous body without undergoing some deformation. The deformation of the material elements (relative to their stressfree shape in the hypothetical, stressfree state N1) induces internal stresses – the socalled growthstresses. 
When subjected to external tractions and to body forces, the uniformlygrown material elements deform. Their deformation (relative to their stressfree shape in state U1 (or, equivalently, U2), generates internal stresses which
 When subjected to external tractions and to body forces, the material elements deform further. Their deformation (relative to their stressfree shape in state N1), generates internal stresses which

We start our analysis with the stressfree (and "ungrown") reference configuration "0" at the top of the diagram, and initially follow the deformation shown in the left half of the sketch. In a first step, each infinitesimal material element in the body is subjected to the same isotropic growth which changes its mass from to . Assuming that the growth process does not change the density of the material, also specifies the volumetric growth of the material. All material elements grow by the same amount, therefore the individual elements can be (re)assembled to form a continuous body (state U2). In this state, the body is stressfree and, compared to the reference configuration "0", all lengths have increased by a factor of where is the body's spatial dimension (i.e. in the sketch above). The covariant basis vectors in this uniformlygrown, stressfree configuration are therefore given by
the metric tensor is given by
and the volume of an infinitesimal material element has increased from
to
We now subject the stressfree, uniformlygrown body to external loads and determine its deformation, using the principle of virtual displacements. Since the uniformly grown state "U2" is stressfree, we may regard it as the reference state for the subsequent deformation. The strain tensor that describes the deformation from the stressfree (and uniformlygrown) state "U2" to the final equilibrium configuration "UE" must therefore be defined as
Equations (11) and (12) allow us to express the principle of virtual displacements in terms of
Note that this equation does not contain any references to quantities in the intermediate states "U1" and "U2".
We will now consider the case of spatially nonuniform growth, illustrated in the right half of the sketch. If the isotropic growth is spatially nonuniform, , the growth will try to expand all infinitesimal material elements isotropically – but each one by a different amount as illustrated by the hypothetical state N1 in which each material element has expanded to its stressfree state in which its metric tensor is given by
Material elements will only be stressfree if the strain
relative to their isotropically grown shape in state N1 is zero. In general, the displacements induced by such an isotropic expansion will be incompatible and it would be impossible to (re)assemble the individually grown material elements to a continuous body unless the material elements undergo some deformation. The elements' deformation relative to their stressfree shape in N1 will generate internal "growthstresses" (stage N2). When subjected to external loads and body forces the body will undergo further deformations until the the stress (generated by the particles' total deformation relative to their stressfree state N1) balances the applied loads.
It is important to realise that, as in the case of spatially uniform growth, the strain defined by (12) is an intrinsic quantity that provides a measure of each particles' local deformation relative to its stressfree shape in N1. The intermediate (and in the current case hypothetical), isotropically grown state N1 does not appear in the analysis – it only serves to define the stressfree shape for each infinitesimal material element. Equation (13) therefore remains valid.
If the problem does not have any symmetries (e.g. axisymmetry) whose exploitation would reduce the spatial dimension of the problem, the most compact form of the equations is obtained by resolving all vectors into a fixed Cartesian basis so that the undeformed position vector is given by
where the are the basis vectors in the direction of the Cartesian Eulerian coordinate axes.
Similarly, we write
and
We use the Eulerian coordinates in the undeformed configuration as the Lagrangian coordinates so that
With this choice, the tangent vectors to the undeformed coordinate lines are the Cartesian basis vectors
and the undeformed metric tensor is the Kronecker delta (the unit matrix)
We now expand the (unknown) deformed position vector in the same basis,
and derive a finite element approximation for this vector field from the principle of virtual displacements. For this purpose we decompose the undeformed body into a number of finite elements, using the standard mesh generation process described previously. This establishes the Eulerian position of the nodes in the body's undeformed configuration. Since the Eulerian coordinates of material points in the undeformed configuration coincide with their Lagrangian coordinates (see (16)), a finiteelement representation of the Lagrangian coordinates is established by writing
where is the th Lagrangian coordinate of global node , and the are the global finiteelement shape functions. In practice, the are, of course, represented by local shape functions, , so that the Lagrangian coordinate at local coordinate in element is given by
where we use the same notation as in the Introduction to the Finite Element Method document.
We employ the same basis functions to represent the components of the unknown vector field , by writing
and treat the (Eulerian) nodal positions as the unknowns. With this discretisation, the variations in correspond to variations in the nodal positions so that
The principle of virtual displacement (13) therefore becomes
[Note that summation convention enforces the summation over the index .] The displacement boundary condition (1) determines the positions of all nodes that are located on the boundary ,
and equation (3) requires their variations to vanish,
The variations of all other nodal positions are arbitrary (and independent of each other), therefore the terms in the curly brackets in (18) must vanish individually. This provides one (discrete) equation for each unknown ,
These equations can again be assembled in an elementbyelement fashion.
We will now discuss how the discrete equations (19) are implemented in oomphlib
. To facilitate the analysis of multiphysics problems, we introduce generalisations of the Node
, FiniteElement
and Mesh
classes which provide separate storage (and access functions) for all solid mechanics data. The resulting SolidFiniteElements
can be used as standalone elements for the simulation of pure solid mechanics problems. More importantly, however, the design makes it easy to employ multiple inheritance to create more complex elements that solve the equations of solid mechanics together with any other field equations. For instance, if we combine a FiniteElement
that solves the unsteady heat equation with a SolidFiniteElement
that describes the elastic deformations, we obtain an element that can be used to simulate unsteady heat conduction in an elastic body that is subject to largeamplitude deformations, say. This is illustrated in one of oomphlib's
multiphysics example codes.
The SolidNode
class is derived from the Node
class and implements the additional functionality required for solid mechanics problems. The key new feature is that each Node
must store its (fixed) Lagrangian coordinate , while its Eulerian position must be regarded as an unknown. This requires the following changes to member functions of the Node
class:
Node::x(...)
returns the Eulerian coordinates of the Node
. Internally, this function accesses the nodal coordinates via pointers to double precision numbers. In solid mechanics problems we must be able to regard the nodal positions as unknowns. In SolidNodes
the nodal positions are therefore stored as values of a (member) Data
object created during construction of the SolidNode
. (As usual, the values can be either unknown or pinned, and can have time histories). The function Data::Is_pinned
(which is set to 1) if a coordinate is pinned.Data
member function Data::pin(...)
which can be used to pin specific nodal values, we provide the function TimeStepper
of the positional Data
. By default, we use the same TimeStepper
for the nodal Data
and for the nodal positions. In multiphysics problems this may not be appropriate, however. Consider, for instance, solving an unsteady heat equation in a dynamically deforming, elastic body. In this problem the 2nd timederivatives of the nodal position might be evaluated by a Newmark
scheme, acting on the history values of the nodal positions, whereas the timederivatives of the temperature might be determined by a BDF
scheme, operating on the history values of the nodal Data
. In such cases, the default assignment for the two timesteppers can be overwritten with the access functions Data::time_stepper_pt()
.Data
. Similarly, the elements' internal Data
is used to store any discontinuous solid pressures.SolidNodes
overload the function Data::assign_eqn_numbers()
with Data
by calling Data::assign_eqn_numbers()
.The class SolidFiniteElement
is derived from FiniteElement
and implements the additional functionality required for solid mechanics problems. Again, most of the additional (or revised) functionality is related to the presence of the two coordinate systems which requires the following changes to FiniteElement
member functions:
SolidFiniteElements
are SolidNodes
, therefore we overload the function FiniteElements::construct_node(...)
to ensure that a SolidNode
with the appropriate amount of storage is built. As in the case of FiniteElements
, the required spatial dimension of the elements' constituent nodes, their number of nodal values etc. are specified via the FiniteElement's
(pure) virtual member functions FiniteElement::required_ndim(...)
, FiniteElement::required_nvalue(...)
, etc, which must be implemented for all specific elements that are derived from the SolidFiniteElement
base class. As discussed above, the constructor of the SolidNodes
requires additional parameters, such as the number of Lagrangian coordinates. These must be specified by implementing SolidFiniteElement::required_nlagrangian(...)
and similar other functions. As in the case of FiniteElements
, many of these functions are already implemented as virtual (rather than pure virtual) functions which provide sensible default values. Such functions must be overloaded in specific derived elements if the default assignments are not appropriate.FiniteElement::interpolated_x(...)
, can remain unchanged because Node::x(...)
always returns the Eulerian nodal positions. SolidFiniteElements
provide additional functions, such as SolidNode
; the second version automatically computes the suitably constrained Lagrangian coordinates if the SolidNode
is hanging.)Data
that represents the nodal positions. We must ensure that these Data
items are included in the various equation numbering schemes. For this purpose we provide the function SolidFiniteElement::assign_solid_local_eqn_numbers()
which sets up the local equation numbering scheme for all solid Data
associated with an element. This function is called when the SolidFiniteElement's local equation numbers are generated.SolidFiniteElements
now form a suitable basis for all elements whose deformation is determined by the equations of solid mechanics (or some variant thereof). To implement a specific solid mechanics element, we must represent its geometry, its state of stress, etc., in terms of the SolidFiniteElement's
positional and nonpositional Data
. This requires the specification of the shape functions and the functions that compute the element's Jacobian matrix and its residual vector – the latter implementing the element's contribution to the global residual vector defined by the discretised principle of virtual displacements, (19). As for "normal" FiniteElements
it is sensible to construct specific SolidFiniteElements
in a hierarchy which separates between the implementation of the governing equations and the representation of the element geometry. For instance, the SolidQElement
family represents the generalisation of the QElement
family to SolidFiniteElements
, while PVDEquations
implement the principle of virtual displacements (19). The two are combined by multiple inheritance to form the QPVDElement
class.PVDEquations
and the PVDEquationsWithPressure
classes.The SolidMesh
class is a generalisation of the Mesh
class whose key additional features are:
Mesh::node_pt(...)
function with SolidNode
, rather than a "normal" Node
. Equivalent access functions are implemented for all other Mesh
member functions that return pointers to Nodes
.Nodes
to their Lagrangian coordinates, thus turning the current configuration into the stressfree reference configuration. This function greatly facilitates the construction of SolidMeshes
via inheritance from existing Meshes
. If, for instance, SomeMesh
is an existing, fully functional Mesh
, the corresponding SolidMesh
can be constructed with a few lines of code, as in this example: To evaluate the load terms
in the discretised form of the variational principle (19) we employ the same strategy as for most other Neumanntype boundary conditions and attach socalled SolidTractionElements
to the appropriate faces of higherdimensional "bulk" solid mechanics elements. Our default implementation allows the load (specified by the "user" via a function pointer that is passed to the SolidTractionElements
) to depend on the Eulerian and Lagrangian coordinates, and on the outer unit normal to the solid. This interface should be sufficiently general for most cases of interest. If additional dependencies are required, it is easy to create new SolidTractionElements
. The use of the SolidTractionElements
is demonstrated in several solid mechanics tutorials.
In timedependent problems, the boundary value problem defined by the variational principle (10) must be augmented by suitable initial conditions which specify the state of the system at time The initial conditions specify the initial shape of the solid body,
and its initial velocity,
where and are given. The accelerations at follow from the solution of (10) and can therefore not be enforced, unless we wish to initialise the timestepping procedure with a known exact solution . (Only!) in this case are we allowed to assign an initial value for the acceleration via
We will assume that timestepping is performed with the Newmark
method which is our default timestepper for hyperbolic problems. In this case the timederivatives of the nodal positions in (10) are replaced by an approximation which involves the current and three "history values" of the nodal positions. To start the timeintegration, we must assign suitable values to these quantities to ensure that the initial state of the system is represented correctly.
To assign the initial values for the nodal positions, we (temporarily) remove all boundary conditions for the nodal positions and determine their initial values by solving equation (20) in its weak form,
where is the th component of Equation (23) provides equations for the components of the initial nodal positions, (where . To determine the initial nodal velocities, we repeat the same procedure with the prescribed velocities and solve
for the initial nodal velocities, (where . Finally, assuming that we have an exact solution for the accelerations, we solve
for the initial nodal accelerations, (where . Having determined the nodal positions and their first and second timederivatives at , we can use the functions Newmark::assign_initial_data_values_stage1(...)
and Newmark::assign_initial_data_values_stage2(...)
to compute the positional history values which ensure that the Newmark approximations for the initial velocity and acceleration are correct. This procedure is fully implemented in the function
whose arguments are:
Problem::assign_eqn_numbers()
.SolidMesh
on which the initial conditions are assigned.TimeStepper
(which has to be a member of the Newmark
family).Here is a brief outline of the implementation: All SolidFiniteElements
store a pointer to a SolidInitialCondition
object. By default this pointer is set to NULL, indicating that FiniteElement::get_residual(...)
should compute the residuals of the "normal" governing equations. SolidFiniteElements
whose initial conditions are to be set with the above function, must redirect the computation of the residual to
whenever the pointer to the SolidInitialCondition
is nonNULL, as illustrated in this code fragment:
The SolidInitialCondition
object stores a (pointer to a) GeomObject
and a flag that indicates which timederivative (0th, 1st or 2nd) of the GeomObject's
position vector is to be assigned to the nodal coordinates. Based on the value of this flag, the function SolidFiniteElement::get_residuals_for_solid_ic(...)
, is able to compute the residuals corresponding to equations (23), (24) or (25).
This all sounds very complicated (and it is!) but luckily all the hard work has already been done and the relevant procedures are fully implemented. Hence, the actual assignment of the initial conditions is as simple as this:
If we do not know an exact solution to our problem (and in most cases we obviously won't...), we can only use the procedure described above to determine the initial nodal positions and velocities. In that case we solve the equations (19) for the remaining "history value". Since the equations (19) are linear in the accelerations, this is a linear problem whose Jacobian matrix is proportional to the mass matrix
The procedure which determines the initial "history values" from the given initial positions and velocities while ensuring consistency with the governing equation at is implemented in SolidMesh::Solid_IC_problem.set_newmark_initial_condition_consistently(...)
which takes the same arguments as the function that assigns the acceleration directly, but also requires a function pointer to a "multiplier" . If there is no growth, i.e. if in (19), the multiplier is given by the timescale ratio ; if the body is subjected to uniform isotropic growth, , the multiplier is equal to . If the wrong multiplier is specified (or if it is omitted, in which case the default value of 1.0 is used) the residuals (19) will be nonzero (or at least larger than the tolerance specified in SolidICProblem
). In this case a warning is issued and the code execution terminates. This behaviour can be suppressed by increasing the tolerance suitably, but you do this at your own risk!
Important: The above procedures can only handle the assignment of initial conditions in problems that are formulated in terms of displacements and do not involve any additional variables such as solid pressures. We do not believe that it is possible to implement the assignment of initial conditions for such problems without additional knowledge about the precise form of the constitutive equations. Therefore we provide a virtual function SolidFiniteElement::has_internal_solid_data()
whose role it is to return a bool that indicates if a specific SolidFiniteElement
stores such data. By default, the function returns false
and should be overloaded in derived elements which involve unknowns that do not represent nodal positions. If the function returns true
for any element that is used during the automatic assignment of initial conditions the code execution stops with an appropriate warning message.
A pdf version of this document is available.