Demo problem: Compression of 2D circular disk

In this example we study the compression of a 2D circular disk, loaded by an external pressure. We also demonstrate:

- how to "upgrade" a
`Mesh`

to a`SolidMesh`

- why it is necessary to use "undeformed MacroElements" to ensure that the numerical results converge to the correct solution under mesh refinement if the domain has curvilinear boundaries.
- how to switch between different constitutive equations
- how to incorporate isotropic growth into the model

We validate the numerical results by comparing them against the analytical solution of the equations of linear elasticity which are valid for small deflections.

The figure below shows a sketch of the basic problem: A 2D circular disk of radius is loaded by a uniform pressure . We wish to compute the disk's deformation for a variety of constitutive equations.

The next sketch shows a variant of the problem: We assume that the material undergoes isotropic growth (e.g. via a biological growth process or thermal expansion, say) with a constant growth factor . We refer to the the document "Solid mechanics: Theory and implementation" for a detailed discussion of the theory of isotropic growth. Briefly, the growth factor defines the relative increase in the volume of an infinitesimal material element, relative to its volume in the stress-free reference configuration. If the growth factor is spatially uniform, isotropic growth leads to a uniform expansion of the material. For a circular disk, uniform growth increases the disk's radius from to without inducing any internal stresses. This uniformly expanded disk may then be regarded as the stress-free reference configuration upon which the external pressure acts.

The animation shows the disk's deformation when subjected to uniform growth of and loaded by a pressure that ranges from negative to positive values. All lengths were scaled on the disks initial radius (i.e. its radius in the absence of growth and without any external load).

The figure below illustrates the disk's load-displacement characteristics by plotting the disk's non-dimensional radius as function of the non-dimensional pressure, , where is the characteristic stiffness of the material, for a variety of constitutive equations.

The blue, dash-dotted line corresponds to `oomph-lib's`

generalisation of Hooke's law (with Young's modulus and Poisson ratio ) in which the dimensionless second Piola Kirchhoff stress tensor (non-dimensionalised with the material's Young's modulus , so that ) is given by

Here is Green's strain tensor, formed from the difference between the deformed and undeformed metric tensors, and , respectively. The three different markers identify the results obtained with the two forms of the principle of virtual displacement, employing the displacement formulation (squares), and a pressure/displacement formulation with a continuous (delta) and a discontinuous (nabla) pressure interpolation.

For zero pressure the disk's non-dimensional radius is equal to the uniformly grown radius For small pressures the load-displacement curve follows the linear approximation

We note that the generalised Hooke's law leads to strain softening behaviour under compression (the pressure required to reduce the disk's radius to a given value increases more rapidly than predicted by the linear approximation) whereas under expansion (for negative external pressures) the behaviour is strain softening.

The red, dashed line illustrates the behaviour when Fung & Tong's generalisation of the Mooney-Rivlin law (with Young's modulus, , Poisson ratio and Mooney-Rivlin parameter ) is used as the constitutive equation. For this constitutive law, the non-dimensional strain energy function , where the characteristic stress is given by Young's modulus, i.e. , is given by

where

is the shear modulus, and and are the three invariants of Green's strain tensor. See "Solid mechanics: Theory and implementation" for a detailed discussion of strain energy functions. The figure shows that for small deflections, the disk's behaviour is again well approximated by linear elasticity. However, in the large-displacement regime the Mooney-Rivlin is strain hardening under extension and softening under compression when compared to the linear elastic behaviour.

As usual we define the global problem parameters in a namespace. We provide pointers to the constitutive equations and strain energy functions to be explored, and define the associated constitutive parameters.

//============namespace_for_problem_parameters=====================

/// Global variables

//=================================================================

namespace Global_Physical_Variables

{

/// Pointer to strain energy function

StrainEnergyFunction* Strain_energy_function_pt;

/// Pointer to constitutive law

ConstitutiveLaw* Constitutive_law_pt;

/// Elastic modulus

double E=1.0;

/// Poisson's ratio

double Nu=0.3;

/// "Mooney Rivlin" coefficient for generalised Mooney Rivlin law

double C1=1.3;

StrainEnergyFunction * Strain_energy_function_pt

Pointer to strain energy function.

double C1

"Mooney Rivlin" coefficient for generalised Mooney Rivlin law

Next we define the pressure load, using the general interface defined in the `SolidTractionElement`

class. The arguments of the function reflect that the load on a solid may be a function of the Lagrangian and Eulerian coordinates, and the external unit normal on the solid. Here we apply a spatially constant external pressure of magnitude `P`

which acts in the direction of the negative outer unit normal on the solid.

/// Uniform pressure

double P = 0.0;

/// Constant pressure load

const Vector<double> &n, Vector<double> &traction)

{

unsigned dim = traction.size();

for(unsigned i=0;i<dim;i++)

{

traction[i] = -P*n[i];

}

} // end of pressure load

void constant_pressure(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)

Constant pressure load.

Finally, we define the growth function and impose a spatially uniform expansion that (in the absence of any external load) would increase the disk's volume by 10%.

/// Uniform volumetric expansion

double Uniform_gamma=1.1;

/// Growth function

{

gamma = Uniform_gamma;

}

} // end namespace

void growth_function(const Vector< double > &xi, double &gamma)

Growth function.

The driver code is very short: We store the command line arguments (as usual, we use a non-zero number of command line arguments as an indication that the code is run in self-test mode and reduce the number of steps performed in the parameter study) and create a strain-energy-based constitutive equation: Fung & Tong's generalisation of the Mooney-Rivlin law.

//=====start_of_main====================================================

/// Driver code for disk-compression

//======================================================================

{

// Store command line arguments

CommandLineArgs::setup(argc,argv);

// Define a strain energy function: Generalised Mooney Rivlin

new GeneralisedMooneyRivlin(&Global_Physical_Variables::Nu,

// Define a constitutive law (based on strain energy function)

new IsotropicStrainEnergyFunctionConstitutiveLaw(

int main(int argc, char *argv[])

Driver code for disk-compression.

We build a problem object, using the displacement-based `RefineableQPVDElements`

to discretise the domain, and perform a parameter study, exploring the disk's deformation for a range of external pressures.

// Case 0: No pressure, generalised Mooney Rivlin

//-----------------------------------------------

{

//Set up the problem

cout << "gen. M.R.: RefineableQPVDElement<2,3>" << std::endl;

//Run the simulation

problem.parameter_study(0);

} // done case 0

Uniform compression of a circular disk in a state of plane strain, subject to uniform growth.

void parameter_study(const unsigned &case_number)

Run simulation: Pass case number to label output files.

We repeat the exercise with elements from the `RefineableQPVDElementWithContinuousPressure`

family which discretise the principle of virtual displacements (PVD) in the pressure/displacement formulation, using continuous pressures (Q2Q1; Taylor Hood).

// Case 1: Continuous pressure formulation with generalised Mooney Rivlin

//------------------------------------------------------------------------

{

//Set up the problem

RefineableQPVDElementWithContinuousPressure<2> > problem;

cout << "gen. M.R.: RefineableQPVDElementWithContinuousPressure<2> "

<< std::endl;

//Run the simulation

problem.parameter_study(1);

} // done case 1

The next computation employs `RefineableQPVDElementWithPressure`

elements in which the pressure is interpolated by piecewise linear but globally discontinuous basis functions (Q2Q-1; Crouzeiux-Raviart).

// Case 2: Discontinuous pressure formulation with generalised Mooney Rivlin

//--------------------------------------------------------------------------

{

//Set up the problem

problem;

cout << "gen. M.R.: RefineableQPVDElementWithPressure<2>" << std::endl;

//Run the simulation

problem.parameter_study(2);

} // done case 2

Next, we change the constitutive equation to `oomph-lib's`

generalised Hooke's law,

// Change the consitutive law: Delete the old one

// Create oomph-lib's generalised Hooke's law constitutive equation

new GeneralisedHookean(&Global_Physical_Variables::Nu,

before repeating the parameter studies with the same three element types:

// Case 3: No pressure, generalised Hooke's law

//----------------------------------------------

{

//Set up the problem

cout << "gen. Hooke: RefineableQPVDElement<2,3> " << std::endl;

//Run the simulation

problem.parameter_study(3);

} // done case 3

// Case 4: Continuous pressure formulation with generalised Hooke's law

//---------------------------------------------------------------------

{

//Set up the problem

RefineableQPVDElementWithContinuousPressure<2> > problem;

cout << "gen. Hooke: RefineableQPVDElementWithContinuousPressure<2> "

<< std::endl;

//Run the simulation

problem.parameter_study(4);

} // done case 4

// Case 5: Discontinous pressure formulation with generalised Hooke's law

//------------------------------------------------------------------------

{

//Set up the problem

cout << "gen. Hooke: RefineableQPVDElementWithPressure<2> " << std::endl;

//Run the simulation

problem.parameter_study(5);

} // done case 5

// Clean up

} // end of main

We formulate the problem in cartesian coordinates (ignoring the problem's axisymmetry) but discretise only one quarter of the domain, applying appropriate symmetry conditions along the x and y axes. The computational domain may be discretised with the `RefineableQuarterCircleSectorMesh`

that we already used in many previous examples. To use the mesh in this solid mechanics problem we must first "upgrade" it to a `SolidMesh`

. This is easily done by multiple inheritance:

//==========start_mesh=================================================

/// Elastic quarter circle sector mesh with functionality to

/// attach traction elements to the curved surface. We "upgrade"

/// the RefineableQuarterCircleSectorMesh to become an

/// SolidMesh and equate the Eulerian and Lagrangian coordinates,

/// thus making the domain represented by the mesh the stress-free

/// configuration.

/// \n\n

/// The member function \c make_traction_element_mesh() creates

/// a separate mesh of SolidTractionElements that are attached to the

/// mesh's curved boundary (boundary 1).

//=====================================================================

template <class ELEMENT>

public virtual RefineableQuarterCircleSectorMesh<ELEMENT>,

public virtual SolidMesh

{

//////////////////////////////////////////////////////////////////// ////////////////////////////////...

The constructor calls the constructor of the underlying `RefineableQuarterCircleSectorMesh`

and sets the Lagrangian coordinates of the nodes to their current Eulerian positions, making the initial configuration stress-free.

public:

/// Constructor: Build mesh and copy Eulerian coords to Lagrangian

/// ones so that the initial configuration is the stress-free one.

ElasticRefineableQuarterCircleSectorMesh<ELEMENT>(GeomObject* wall_pt,

const double& xi_lo,

const double& fract_mid,

const double& xi_hi,

TimeStepper* time_stepper_pt=

&Mesh::Default_TimeStepper) :

RefineableQuarterCircleSectorMesh<ELEMENT>(wall_pt,xi_lo,fract_mid,xi_hi,

time_stepper_pt)

{

/// Make the current configuration the undeformed one by

/// setting the nodal Lagrangian coordinates to their current

/// Eulerian ones

set_lagrangian_nodal_coordinates();

}

We also provide a helper function that creates a mesh of `SolidTractionElements`

which are attached to the curved domain boundary (boundary 1). These elements will be used to apply the external pressure load.

/// Function to create mesh made of traction elements

void make_traction_element_mesh(SolidMesh*& traction_mesh_pt)

{

// Make new mesh

traction_mesh_pt = new SolidMesh;

// Loop over all elements on boundary 1:

unsigned b=1;

unsigned n_element = this->nboundary_element(b);

for (unsigned e=0;e<n_element;e++)

{

// The element itself:

FiniteElement* fe_pt = this->boundary_element_pt(b,e);

// Find the index of the face of element e along boundary b

int face_index = this->face_index_at_boundary(b,e);

// Create new element

traction_mesh_pt->add_element_pt(new SolidTractionElement<ELEMENT>

(fe_pt,face_index));

}

}

};

The definition of the Problem class is very straightforward. In addition to the constructor and the (empty) `actions_before_newton_solve()`

and `actions_after_newton_solve()`

functions, we provide the function `parameter_study(...)`

which performs a parameter study, computing the disk's deformation for a range of external pressures. The member data includes pointers to the mesh of "bulk" solid elements, and the mesh of `SolidTractionElements`

that apply the pressure load. The trace file is used to document the disk's load-displacement characteristics by plotting the radial displacement of the nodes on the curvilinear boundary, pointers to which are stored in the vector `Trace_node_pt`

.

//======================================================================

/// Uniform compression of a circular disk in a state of plane strain,

/// subject to uniform growth.

//======================================================================

template<class ELEMENT>

{

public:

/// Constructor:

/// Run simulation: Pass case number to label output files

/// Doc the solution

void doc_solution(DocInfo& doc_info);

/// Update function (empty)

void actions_after_newton_solve() {}

/// Update function (empty)

void actions_before_newton_solve() {}

private:

/// Trace file

ofstream Trace_file;

/// Vector of pointers to nodes whose position we're tracing

Vector<Node*> Trace_node_pt;

/// Pointer to solid mesh

/// Pointer to mesh of traction elements

SolidMesh* Traction_mesh_pt;

};

SolidMesh * Traction_mesh_pt

Pointer to mesh of traction elements.

ElasticRefineableQuarterCircleSectorMesh< ELEMENT > * Solid_mesh_pt

Pointer to solid mesh.

Vector< Node * > Trace_node_pt

Vector of pointers to nodes whose position we're tracing.

We start by constructing the mesh of "bulk" `SolidElements`

, using the `Ellipse`

object to specify the shape of the curvilinear domain boundary.

//======================================================================

/// Constructor:

//======================================================================

template<class ELEMENT>

{

// Build the geometric object that describes the curvilinear

// boundary of the quarter circle domain

Ellipse* curved_boundary_pt = new Ellipse(1.0,1.0);

// The curved boundary of the mesh is defined by the geometric object

// What follows are the start and end coordinates on the geometric object:

double xi_lo=0.0;

double xi_hi=2.0*atan(1.0);

// Fraction along geometric object at which the radial dividing line

// is placed

double fract_mid=0.5;

//Now create the mesh using the geometric object

Solid_mesh_pt = new ElasticRefineableQuarterCircleSectorMesh<ELEMENT>(

curved_boundary_pt,xi_lo,fract_mid,xi_hi);

Next we choose the nodes on the curvilinear domain boundary (boundary 1) as the nodes whose displacement we document in the trace file.

// Setup trace nodes as the nodes on boundary 1 (=curved boundary)

// in the original mesh.

unsigned n_boundary_node = Solid_mesh_pt->nboundary_node(1);

Trace_node_pt.resize(n_boundary_node);

for(unsigned j=0;j<n_boundary_node;j++)

{Trace_node_pt[j]=Solid_mesh_pt->boundary_node_pt(1,j);}

The `QuarterCircleSectorMesh`

that forms the basis of the "bulk" mesh contains only three elements – not enough to expect the solution to be accurate. Therefore we apply one round of uniform mesh refinement before attaching the `SolidTractionElements`

to the mesh boundary 1, using the function `make_traction_element_mesh()`

in the `ElasticRefineableQuarterCircleSectorMesh`

.

// Refine the mesh uniformly

Solid_mesh_pt->refine_uniformly();

// Now construct the traction element mesh

Solid_mesh_pt->make_traction_element_mesh(Traction_mesh_pt);

We add both meshes to the `Problem`

and build a combined global mesh:

// Solid mesh is first sub-mesh

add_sub_mesh(Solid_mesh_pt);

// Traction mesh is second sub-mesh

add_sub_mesh(Traction_mesh_pt);

// Build combined "global" mesh

build_global_mesh();

Symmetry boundary conditions along the horizontal and vertical symmetry lines require that the nodes' vertical position is pinned along boundary 0, while their horizontal position is pinned along boundary 2.

// Pin the left edge in the horizontal direction

unsigned n_side = mesh_pt()->nboundary_node(2);

for(unsigned i=0;i<n_side;i++)

{Solid_mesh_pt->boundary_node_pt(2,i)->pin_position(0);}

// Pin the bottom in the vertical direction

unsigned n_bottom = mesh_pt()->nboundary_node(0);

for(unsigned i=0;i<n_bottom;i++)

{Solid_mesh_pt->boundary_node_pt(0,i)->pin_position(1);}

Since we are using refineable solid elements, we pin any "redundant" pressure degrees of freedom in the "bulk" solid mesh (see the exercises in another tutorial for a more detailed discussion of this issue).

// Pin the redundant solid pressures (if any)

PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(

Solid_mesh_pt->element_pt());

Next, we complete the build of the elements in the "bulk" solid mesh by passing the pointer to the constitutive equation and the pointer to the isotropic-growth function to the elements:

//Complete the build process for elements in "bulk" solid mesh

unsigned n_element =Solid_mesh_pt->nelement();

for(unsigned i=0;i<n_element;i++)

{

//Cast to a solid element

ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Solid_mesh_pt->element_pt(i));

// Set the constitutive law

el_pt->constitutive_law_pt() =

// Set the isotropic growth function pointer

el_pt->isotropic_growth_fct_pt()=Global_Physical_Variables::growth_function;

}

We repeat this exercise for the `SolidTractionElements`

which must be given a pointer to the function that applies the pressure load

// Complete build process for SolidTractionElements

n_element=Traction_mesh_pt->nelement();

for(unsigned i=0;i<n_element;i++)

{

//Cast to a solid traction element

SolidTractionElement<ELEMENT> *el_pt =

dynamic_cast<SolidTractionElement<ELEMENT>*>

(Traction_mesh_pt->element_pt(i));

//Set the traction function

el_pt->traction_fct_pt() = Global_Physical_Variables::constant_pressure;

}

Finally, we set up the equation numbering scheme and report the number of unknowns.

//Set up equation numbering scheme

cout << "Number of equations: " << assign_eqn_numbers() << std::endl;

}

The post-processing function outputs the shape of the deformed disk. We use the trace file to record how the disk's volume (area) and the radii of the control nodes on the curvilinear domain boundary vary with the applied pressure. To facilitate the validation of the results against the analytical solution, we also add the radius predicted by the linear theory to the trace file.

//==================================================================

/// Doc the solution

//==================================================================

template<class ELEMENT>

void StaticDiskCompressionProblem<ELEMENT>::doc_solution(DocInfo& doc_info)

{

ofstream some_file;

char filename[100];

// Number of plot points

unsigned npts = 5;

// Output shape of deformed body

sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),

doc_info.number());

some_file.open(filename);

Solid_mesh_pt->output(some_file,npts);

some_file.close();

//Find number of solid elements

unsigned nelement = Solid_mesh_pt->nelement();

// Work out volume

double volume = 0.0;

for(unsigned e=0;e<nelement;e++)

{volume+= Solid_mesh_pt->finite_element_pt(e)->size();}

// Exact outer radius for linear elasticity

double nu=Global_Physical_Variables::Nu;

double exact_r=sqrt(Global_Physical_Variables::Uniform_gamma)*

*((1.0+nu)*(1.0-2.0*nu)));

// Write trace file: Problem parameters

Trace_file << Global_Physical_Variables::P << " "

<< Global_Physical_Variables::Uniform_gamma << " "

<< volume << " "

<< exact_r << " ";

// Write radii of trace nodes

unsigned ntrace_node=Trace_node_pt.size();

for (unsigned j=0;j<ntrace_node;j++)

{

Trace_file << sqrt(pow(Trace_node_pt[j]->x(0),2)+

pow(Trace_node_pt[j]->x(1),2)) << " ";

}

Trace_file << std::endl;

} // end of doc_solution

The function `parameter_study(...)`

computes the disk's deformation for a range of external pressures and outputs the results. The output directory is labelled by the `unsigned`

function argument. This ensures that parameter studies performed with different constitutive equations are written into different directories.

//==================================================================

/// Run the paramter study

//==================================================================

template<class ELEMENT>

const unsigned& case_number)

{

// Output

DocInfo doc_info;

char dirname[100];

sprintf(dirname,"RESLT%i",case_number);

// Set output directory

doc_info.set_directory(dirname);

// Step number

doc_info.number()=0;

// Open trace file

char filename[100];

sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());

Trace_file.open(filename);

//Parameter incrementation

double delta_p=0.0125;

unsigned nstep=21;

// Perform fewer steps if run as self-test (indicated by nonzero number

// of command line arguments)

if (CommandLineArgs::Argc!=1)

{

nstep=3;

}

// Offset external pressure so that the computation sweeps

// over a range of positive and negative pressures

Global_Physical_Variables::P =-delta_p*double(nstep-1)*0.5;

// Do the parameter study

for(unsigned i=0;i<nstep;i++)

{

//Solve the problem for current load

newton_solve();

// Doc solution

doc_solution(doc_info);

doc_info.number()++;

// Increment pressure load

Global_Physical_Variables::P += delta_p;

}

} // end of parameter study

Recall how `oomph-lib`

employs `MacroElements`

to represent the exact domain shapes in adaptive computations involving problems with curvilinear boundaries. When an element is refined, the (Eulerian) position of any newly-created nodes is based on the element's `MacroElement`

counterpart, rather than being determined by finite-element interpolation from the "father element". This ensures that **(i)** newly-created nodes on curvilinear domain boundaries are placed exactly onto those boundaries and **(ii)** that newly-created nodes in the interior are placed at positions that match smoothly onto the boundary.

This strategy is adapted slightly for solid mechanics problems:

- The Eulerian position of newly-created
`SolidNodes`

is determined by finite element interpolation from the "father element",**unless**the newly-created`SolidNode`

is located on a domain boundary and its position is pinned by displacement boundary conditions.

- The same procedure is employed to determine the Lagrangian coordinates of newly-created
`SolidNodes`

.

These modifications ensure that, as before, newly-created nodes on curvilinear domain boundaries are placed exactly onto those boundaries if their positions are pinned by displacement boundary conditions. (If the nodal positions are not pinned, the node's Eulerian position will be determined as part of the solution.) The use of finite-element interpolation from the "father element" in the interior of the domain for both Lagrangian and Eulerian coordinates ensures that the creation of new nodes does not induce any stresses into a previously computed solution.

- Our discretisation of the problem in cartesian coordinates did not exploit the problem's axisymmetry. Examine the trace file to assess to which extent the computation retained the axisymmetry.

A pdf version of this document is available.