Demo problem: Compressible and incompressible behaviour

The purpose of this tutorial is to discuss the use of compressible and incompressible constitutive equations for solid mechanics problems. As discussed in the Solid Mechanics Theory Tutorial, problems in which the solid is truly incompressible require a solution based on the pressure/displacement formulation of the principle of virtual displacements (PVD). Mathematically, this is because the pressure acts as the Lagrange multiplier that enforces incompressibility. The pressure/displacement form of the PVD may be discretised with finite elements that employ continuous (e.g. Taylor-Hood-type) or discontinuous (e.g. Crouzeix-Raviart-type) elements.

Some constitutive equations allow for compressible and incompressible behaviour, depending on their parameters. Problems involving such constitutive equations may be solved with large number of formulations. As an example, consider `oomph-lib's`

generalised Hookean constitutive equation (see the disclaimer below).

** I. The displacement form **

`oomph-lib's`

generalised Hookean constitutive law assumes that the 2nd Piola Kirchhoff stress (non-dimensionalised on Young's modulus ) is related to Green's strain tensor via

Here and are the metric tensors associated with the undeformed and deformed configurations. This constitutive law reduces to the the classical version of Hooke's law for small strains when In the above form, the constitutive law can be used directly in the displacement-based form of the PVD, unless – the case that corresponds to incompressible behaviour in the small-displacement regime. (While the above formulation only breaks down completely when numerical solutions based on the above form of the constitutive equation become badly-behaved as approaches that value.)

** II. The pressure-displacement form for compressible and incompressible behaviour**

To avoid the problems that arise as the constitutive equation may be rewritten in the form

where the deviatoric stress is given by

To remain consistent with (1), the pressure must then be determined via the equation

where

is the generalised dilatation (which reduces to the actual dilatation in the small-strain limit). The inverse bulk modulus is defined as

and tends to zero as The alternative form of the constitutive equation can therefore be used for any value of . If the constraint (2) enforces i.e incompressible behaviour (at least in the small displacement regime – for large deflections the constraint (2) simply ensures that the generalised dilatation vanishes; this may not be a physically meaningful constraint).

** III. Truly incompressible behaviour**

Finally, we may retain the decomposition of the stress into its deviatoric and non-deviatoric parts but determine the pressure from the actual incompressibility constraint, i.e. replace (2) by

In this case, the deviatoric part of the stress is determined by the material's constitutive parameters, i.e. its Young's modulus and its Poisson ratio while the pressure (acting as the Lagrange multiplier for the constraint (3)) enforces true incompressibility.

Disclaimer |

With this disclaimer in mind, we will now demonstrate the use of the various combinations outlined above in a simple test problem for which an exact (linearised) solution is available.

Here is a sketch of the problem: Three faces of a square elastic body are surrounded by "slippery" rigid walls that allow the body to slide freely along them. The body is subject to a vertical, gravitational body force, acting in the negative -direction.

We choose the height of the square as the reference length for the non-dimensionalisation by setting and use the characteristic stiffness associated with the body's constitutive equation to scale the stresses and the body forces. For instance, for linear elastic behaviour, we choose the reference stress to be the solid's Young's modulus, thus setting The non-dimensional body force is then given by where

indicates the magnitude of the gravitational load relative to the body's stiffness; see the Solid Mechanics Theory Tutorial for full details on the non-dimensionalisation of the governing equations.

Assuming weak loading, i.e. , the body will undergo small deflections and its deformation will be governed by the equations of linear elasticity. Given that the walls bounding the solid are slippery, it is easy to show that a displacement field has the form . Inserting this ansatz into the Navier-Lame equations shows that the non-dimensional displacement in the vertical direction is given by

The stresses are given by

This shows that, as the solid approaches incompressibility, the displacements are suppressed and the stress distribution becomes hydrostatic with

Here is a plot of the three non-zero quantities (the vertical displacement , and the horizontal and vertical stresses, and , respectively) for and The shaded surface shows the exact solution while the mesh shows the finite element solution obtained with the displacement formulation of the problem. The solutions obtained with formulations II and III are graphically indistinguishable. All other quantities (the transverse displacement and the shear stress ) are of order

The driver code listed below also computes solutions for various other values of (including the incompressible case ), and demonstrates how to enforce "true" incompressibility via formulation III. Furthermore, it computes a solution based on the incompressible Mooney-Rivlin constitutive law. The agreement between analytical and computed results is as good as the one shown in the above plot. In particular in all cases where incompressibility is enforced, the vertical displacement is suppressed and the stress state becomes hydrostatic, as predicted by the analytical solution.

Since we wish to compare the stresses and displacements against the analytical solution, we modify the `SolidElement's`

output function by means of a wrapper class that overloads the `SolidElement::output(...)`

.

//=================start_wrapper==================================

/// Wrapper class for solid element to modify the output

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

template <class ELEMENT>

class MySolidElement : public virtual ELEMENT

{

public:

/// Constructor: Call constructor of underlying element

MySolidElement() : ELEMENT() {}

/// Overload output function

void output(std::ostream &outfile, const unsigned &n_plot)

{

Vector<double> s(2);

Vector<double> x(2);

Vector<double> xi(2);

DenseMatrix<double> sigma(2,2);

//Tecplot header info

outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

//Loop over plot points

for(unsigned l2=0;l2<n_plot;l2++)

{

s[1] = -1.0 + l2*2.0/(n_plot-1);

for(unsigned l1=0;l1<n_plot;l1++)

{

s[0] = -1.0 + l1*2.0/(n_plot-1);

// Get Eulerian and Lagrangian coordinates and the stress

this->interpolated_x(s,x);

this->interpolated_xi(s,xi);

this->get_stress(s,sigma);

//Output the x,y coordinates

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

{outfile << x[i] << " ";}

// Output displacements, the difference between Eulerian and Lagrangian

// coordinates

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

{outfile << x[i]-xi[i] << " ";}

//Output stress

outfile << sigma(0,0) << " "

<< sigma(1,0) << " "

<< sigma(1,1) << " "

<< std::endl;

}

}

}

};

As usual we define the various problem parameters in a global namespace. We prepare a pointer to a constitutive equation and define Poisson's ratio for use with the generalised Hookean constitutive law.

//=======start_namespace==========================================

/// Global variables

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

namespace Global_Physical_Variables

{

/// Pointer to constitutive law

ConstitutiveLaw* Constitutive_law_pt=0;

/// Poisson's ratio for Hooke's law

double Nu=0.45;

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

ConstitutiveLaw * Constitutive_law_pt

Pointer to constitutive law.

We also prepare a pointer to a strain energy function and define the coefficients for the Mooney Rivlin law:

/// Pointer to strain energy function

StrainEnergyFunction* Strain_energy_function_pt=0;

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

double C1=1.3;

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

double C2=1.3;

StrainEnergyFunction * Strain_energy_function_pt

Pointer to strain energy function.

double C1

First "Mooney Rivlin" coefficient for generalised Mooney Rivlin law.

double C2

Second "Mooney Rivlin" coefficient for generalised Mooney Rivlin law.

Finally, we define the gravitational body force.

/// Non-dim gravity

double Gravity=0.0;

/// Non-dimensional gravity as body force

const Vector<double> &xi,

Vector<double> &b)

{

b[0]=0.0;

b[1]=-Gravity;

}

} //end namespace

void gravity(const double &time, const Vector< double > &xi, Vector< double > &b)

Non-dimensional gravity as body force.

The driver code solves the problem with a large number of different formulations and constitutive equations. We start with the generalised Hookean constitutive equation and consider three different values of Poisson's ratio, corresponding to compressible, near-incompressible and incompressible behaviour:

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

/// Driver for compressed square

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

int main()

{

//Flag to indicate if we want the material to be incompressible

bool incompress=false;

// Label for different cases

int case_number=-1;

// Generalised Hookean constitutive equations

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

{

// Nu values

Vector<double> nu_value(3);

nu_value[0]=0.45;

nu_value[1]=0.499999;

nu_value[2]=0.5;

// Loop over nu values

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

{

// Set Poisson ratio

Global_Physical_Variables::Nu=nu_value[i];

std::cout << "===========================================\n";

std::cout << "Doing Nu=" << Global_Physical_Variables::Nu

<< std::endl;

std::cout << "===========================================\n";

// Create constitutive equation

new GeneralisedHookean(&Global_Physical_Variables::Nu);

First, we solve the problem in the pure displacement formulation, using (the wrapped version of the) displacement-based `QPVDElements`

. (As discussed above, the displacement formulation cannot be used for .)

// Displacement-based formulation

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

// (not for nu=1/2, obviously)

if (i!=2)

{

case_number++;

std::cout

<< "Doing Generalised Hookean with displacement formulation: Case "

<< case_number << std::endl;

//Set up the problem with pure displacement based elements

incompress=false;

problem(incompress);

// Run it

problem.run_it(case_number,incompress);

}

Next we consider the pressure/displacement formulation with continuous pressures (Taylor-Hood), using the

`QPVDElementWithContinuousPressure`

element.

// Generalised Hookean with continuous pressure, compressible

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

{

case_number++;

std::cout

<< "Doing Generalised Hookean with displacement/cont pressure "

<< "formulation: Case " << case_number << std::endl;

//Set up the problem with continous pressure/displacement

incompress=false;

CompressedSquareProblem<MySolidElement<

QPVDElementWithContinuousPressure<2> > >

problem(incompress);

// Run it

problem.run_it(case_number,incompress);

}

We suppress the listing of the remaining combinations (see the driver code compressed_square.cc for details):

- Pressure/displacement formulation with discontinuous pressures (Crouzeix-Raviart), using the
`QPVDElementWithPressure`

element. - Pressure/displacement formulation with continuous pressures (Taylor-Hood), using the
`QPVDElementWithContinuousPressure`

element, with true incompressibility enforced via formulation III. - Pressure/displacement formulation with discontinuous pressures (Crouzeix-Raviart), using the
`QPVDElementWithPressure`

element, with true incompressibility enforced via formulation III.

Before the end of the loop over the different values we delete the constitutive equation, allowing it to be re-built with a different Poisson ratio.

// Clean up

}

} // end generalised Hookean

Next we build the strain-energy-based Mooney-Rivlin constitutive law

// Incompressible Mooney Rivlin

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

{

// Create strain energy function

new MooneyRivlin(&Global_Physical_Variables::C1,

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

new IsotropicStrainEnergyFunctionConstitutiveLaw(

and solve the problem with `QPVDElementWithContinuousPressure`

and `QPVDElementWithPressure`

elements, enforcing true incompressibility enforced via formulation III by setting the `incompressible`

flag to true.

// Mooney-Rivlin with continous pressure/displacement;

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

// incompressible

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

{

case_number++;

std::cout

<< "Doing Mooney Rivlin with cont pressure formulation: Case "

<< case_number << std::endl;

//Set up the problem with continous pressure/displacement

incompress=true;

CompressedSquareProblem<MySolidElement<

QPVDElementWithContinuousPressure<2> > >

problem(incompress);

// Run it

problem.run_it(case_number,incompress);

}

// Mooney-Rivlin with discontinous pressure/displacement;

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

// incompressible

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

{

case_number++;

std::cout

<< "Doing Mooney Rivlin with discont pressure formulation: Case "

<< case_number << std::endl;

//Set up the problem with discontinous pressure/displacement

incompress=true;

CompressedSquareProblem<MySolidElement<

QPVDElementWithPressure<2> > >

problem(incompress);

// Run it

problem.run_it(case_number,incompress);

}

// Clean up

}

} //end of main

The `Problem`

class has the usual member functions. The `i_case`

label is used to distinguish different cases while the boolean `incompressible`

indicates if we wish to enforce incompressibility via formulation III.

//=============begin_problem============================================

/// Problem class

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

template<class ELEMENT>

{

public:

/// Constructor: Pass flag that determines if we want to use

/// a true incompressible formulation

/// Update function (empty)

void actions_after_newton_solve() {}

/// Update function (empty)

void actions_before_newton_solve() {}

/// Doc the solution & exact (linear) solution for compressible

/// or incompressible materials

/// Run the job -- doc in RESLTi_case

private:

/// Trace file

ofstream Trace_file;

/// Pointers to node whose position we're tracing

Node* Trace_node_pt;

/// DocInfo object for output

DocInfo Doc_info;

};

Node * Trace_node_pt

Pointers to node whose position we're tracing.

void run_it(const int &i_case, const bool &incompress)

Run the job – doc in RESLTi_case.

void doc_solution(const bool &incompress)

Doc the solution & exact (linear) solution for compressible or incompressible materials.

CompressedSquareProblem(const bool &incompress)

Constructor: Pass flag that determines if we want to use a true incompressible formulation.

We start by building the mesh – the `SolidMesh`

version of the `ElasticRectangularQuadMesh`

.

//===========start_of_constructor=======================================

/// Constructor: Pass flag that determines if we want to enforce

/// incompressibility

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

template<class ELEMENT>

const bool& incompress)

{

// Create the mesh

// # of elements in x-direction

unsigned n_x=5;

// # of elements in y-direction

unsigned n_y=5;

// Domain length in x-direction

double l_x=1.0;

// Domain length in y-direction

double l_y=1.0;

//Now create the mesh

mesh_pt() = new ElasticRectangularQuadMesh<ELEMENT>(n_x,n_y,l_x,l_y);

We complete the build of the elements by specifying the constitutive equation and the gravitational body force.

//Assign the physical properties to the elements

unsigned n_element=mesh_pt()->nelement();

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

{

//Cast to a solid element

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

// Set the constitutive law

el_pt->constitutive_law_pt() =

//Set the body force

el_pt->body_force_fct_pt() = Global_Physical_Variables::gravity;

If the element is based on the pressure/displacement form of the principle of virtual displacements we enforce (true) incompressibility if required. Note that, by default, `oomph-lib's`

pressure/displacement-based solid mechanics elements do not assume incompressibility, i.e. the default is to use formulation II.

// Is the element based on the pressure/displacement formulation?

PVDEquationsWithPressure<2>* test_pt =

dynamic_cast<PVDEquationsWithPressure<2>*>(mesh_pt()->element_pt(i));

if (test_pt!=0)

{

// Do we want true incompressibility (formulation III in the

// associated tutorial) or not (formulation II)

if (incompress)

{

test_pt->set_incompressible();

}

else

{

// Note that this assignment isn't strictly necessary as it's the

// default setting, but it doesn't do any harm to be explicit.

test_pt->set_compressible();

}

}

} // end compressibility

Finally, we pick a control node to document the solid's load-displacement characteristics and apply the boundary conditions (no displacements normal to the "slippery" walls) before setting up the equation numbering scheme.

// Choose a control node: The last node in the solid mesh

unsigned nnod=mesh_pt()->nnode();

Trace_node_pt=mesh_pt()->node_pt(nnod-1);

// Pin the left and right boundaries (1 and 2) in the horizontal directions

for (unsigned b=1;b<4;b+=2)

{

unsigned nnod = mesh_pt()->nboundary_node(b);

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

{

dynamic_cast<SolidNode*>(

mesh_pt()->boundary_node_pt(b,i))->pin_position(0);

}

}

// Pin the bottom boundary (0) in the vertical direction

unsigned b=0;

{

unsigned nnod= mesh_pt()->nboundary_node(b);

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

{

dynamic_cast<SolidNode*>(

mesh_pt()->boundary_node_pt(b,i))->pin_position(1);

}

}

//Assign equation numbers

assign_eqn_numbers();

} //end of constructor

The post-processing routine documents the load-displacement characteristics in a trace file and outputs the deformed domain shape.

//==============start_doc===========================================

/// Doc the solution

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

template<class ELEMENT>

{

ofstream some_file;

char filename[100];

// Number of plot points

unsigned n_plot = 5;

// Output shape of and stress in deformed body

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

Doc_info.number());

some_file.open(filename);

mesh_pt()->output(some_file,n_plot);

some_file.close();

// Write trace file: Load/displacement characteristics

Trace_file << Global_Physical_Variables::Gravity << " "

<< Trace_node_pt->x(0) << " "

<< Trace_node_pt->x(1) << " "

<< std::endl;

We then output the exact solution of the linearised equations using the same format as in the element's overloaded output function.

// Output exact solution for linear elasticity

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

sprintf(filename,"%s/exact_soln%i.dat",Doc_info.directory().c_str(),

Doc_info.number());

some_file.open(filename);

unsigned nelem=mesh_pt()->nelement();

Vector<double> s(2);

Vector<double> x(2);

DenseMatrix<double> sigma(2,2);

// Poisson's ratio

double nu=Global_Physical_Variables::Nu;

if (incompress) nu=0.5;

// Loop over all elements

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

{

//Cast to a solid element

ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));

//Tecplot header info

some_file << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

//Loop over plot points

for(unsigned l2=0;l2<n_plot;l2++)

{

s[1] = -1.0 + l2*2.0/(n_plot-1);

for(unsigned l1=0;l1<n_plot;l1++)

{

s[0] = -1.0 + l1*2.0/(n_plot-1);

// Get Lagrangian coordinates

el_pt->interpolated_x(s,x);

// Output the x,y,..

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

{some_file << x[i] << " ";}

// Exact vertical displacement

double v_exact=Global_Physical_Variables::Gravity*

(1.0+nu)*(1.0-2*nu)/(1.0-nu)*(0.5*x[1]*x[1]-x[1]);

// x and y displacement

some_file << "0.0 " << v_exact << " ";

// Stresses

sigma(0,0)=nu/(1.0-nu)*Global_Physical_Variables::Gravity*(x[1]-1.0);

sigma(1,0)=0.0;

sigma(1,1)=Global_Physical_Variables::Gravity*(x[1]-1.0);

// Output linear stress tensor

some_file << sigma(0,0) << " "

<< sigma(1,0) << " "

<< sigma(1,1) << " "

<< std::endl;

}

}

}

some_file.close();

// Increment label for output files

Doc_info.number()++;

} //end doc

As usual, `oomph-lib`

provides self-tests that assess if the enforcement incompressibility (or the lack thereof) is consistent:

- The compiler will not allow the user to enforce incompressibility on elements that are based on the displacement form of the principle of virtual displacements. This is because the displacement-based elements do not have a member function
`incompressible()`

.

- Certain constitutive laws, such as the Mooney-Rivlin law used in the present example require an incompressible formulation. If
`oomph-lib`

is compiled with the`PARANOID`

flag, an error is thrown if such a constitutive law is used by an element for which incompressibility has not been requested.

**Recall that the default setting is not to enforce incompressibility!**

If the library is compiled without the`PARANOID`

flag no warning will be issued but the results will be "wrong" at least in the sense that the material does not behave like an incompressible Mooney-Rivlin solid. In fact, it is likely that the Newton solver will diverge. Anyway, as we keep saying, without the`PARANOID`

flag, you're on your own!

You should experiment with different combinations of constitutive laws and element types to familiarise yourself with these issues.

- The source files for this tutorial are located in the directory:

demo_drivers/solid/compressed_square/

- The driver code is:

demo_drivers/solid/compressed_square/compressed_square.cc

A pdf version of this document is available.