Solid Mechanics

Elastomechanics

Elastomechanics describes how a solid body \( \Omega \subset \mathbb{R}^d \) deforms under given external forces acting on it. The primal variable is the displacement field \( u \colon \Omega \rightarrow \mathbb{R}^d \), which maps each material point \( x \in \Omega \) of the body to a displacement vector \( u(x) \). The displacement is consideres as an element of a Hilbert space, typically \((H^1(\Omega))^d\). The deformation of the body is given by \( \phi(x) = x + u(x) \). The strain tensor \( \epsilon \in \mathbb{R}^{d \times d} \) relates the strain to the displacement gradient \( \nabla u\) and is defind by \[ \epsilon(u) = \frac{1}{2}( \nabla u + \nabla u^T + \nabla u^T \nabla u). \] In the case of small deformations we can disregard the nonlinear part and work with the linearized strain tensor \( \epsilon(u) \approx \frac{1}{2}( \nabla u + \nabla u^T ) \).

The stress tensor \( \sigma \in \mathbb{R}^{d \times d} \) and its relationship to the strain tensor are also crucial to the problems of elastomechanics. The specific relation between the stress and the strain is called consecutive law or material model and it’s most prominent representative is Hook’s law, which is a linear relation between stress and strain: \[ \sigma(x) = C \colon \epsilon(u(x)), \qquad \text{with} \ C \colon \mathbb{R}^{d\times d} \rightarrow \mathbb{R}^{d\times d}. \] Here \(C\) is a fourth order tensor called stiffness tensor. The displacement field \(u\) is a minimizer of the elastic energy :

\[ \mathcal{E}(u) = \frac{1}{2} \int_\Omega \epsilon(u) \colon \sigma(u) \ dx \ + \underbrace{ \int _{\Omega} f \cdot u \ dx + \int _{\partial \Omega} h \cdot u \ ds(x)} _{\text{external forces}} \quad \longrightarrow \ \min ! \]

The corresponding PDE in strong form is obtained by forming the first order optimality condition of the energy minimization, i.e. by \( \mathcal{E}’(u)[v] = 0 \quad \forall v \in H^1(\Omega)\), see Weiser2016.

In the case of Hook’s law we can introduce two parameters \( \lambda \) and \( \mu \) called Lamé coefficients, which fully describe the material law. In this case we get \( \frac{1}{2} \epsilon \colon \sigma = \frac{\lambda}{2} (\text{tr} \ \epsilon)^2 + \mu \ \epsilon \colon \epsilon \) and the elastic energy can be writte as \[ \mathcal{E}(u) = \int_\Omega \frac{\lambda}{2} (\text{tr} \ \epsilon)^2 + \mu \ \epsilon \colon \epsilon \ dx \ + \ \text{external forces}. \]

Cantilever Beam

In this tutorial we study an elastic beam clamped at one end and subjected to a boundary force. The boundary conditions are given on \(\partial\Omega = \Gamma_0 \cup \Gamma_{h} \cup \Gamma_{h_0} \). We apply homogeneous Dirichlet conditions on one side of the bar to model the clamping and inhomogeneous Neumann boundary conditions to express the boundary force \(h\). On the remaining boundary we apply homogeneous Neumann boundary conditions:

\[ u = 0 \quad \text{on} \ \ \Gamma_0, \qquad n^T \nabla u = h \quad \text{on} \ \ \Gamma_{h} \quad \text{and} \quad n^T \nabla u = 0 \quad \text{on} \ \ \Gamma_{h_0} \qquad \]

Having layed the groundwork, we can take a look at the code.

Domain and boundary caches

With Kaskade 7, differential equations or energy functions describing the elastic transformation of the beam materials are provided by the software and recognized with the material parameter(s). Therefore, we will simply call d0, d1 and d2 under functional.moduli in the domain cache.

On the other hand, we will need to specify boundary conditions for our specific problem, namely the fact that the left hand side of the beam is fixed and that the upper side of the beam is subjected to a downward force.

\[ \text{d0: } G(u) = \alpha (u-u_0)^2 - \beta u. \]

The first term with parameter \(\alpha\) is a penalty term handling the Dirichlet boundary condition for the left hand side of the beam. The second term with parameter \(\beta\) is associated with handling the downward force acted on the beam.

The Fréchet derivatives are computed accordingly.

\[ \begin{aligned} &\text{d1: } G’(u)[v] = \alpha (u-u_0)v - \beta v \\ &\text{d2: } G’’(u)[v,w] = \alpha vw \end{aligned} \]

Initial setup and parameters definition

The program files will be simple_beam.cpp and simple_beam.hh with simple_beam.hh containing the definition of the functional class ElasticityFunctional. We will deal with the user-inputs, set up grids and variable space; and output the computation results in the main script simple_beam.cpp.

The beginning of simple_beam.cpp should contain a list of include directives specifying programs from Kaskade 7 to be used:

#include <iostream>

#include "dune/grid/config.h"
#include "dune/grid/uggrid.hh"
#include <dune/istl/operators.hh>

#include "fem/assemble.hh"
#include "fem/spaces.hh"
#include "fem/diffops/elasto.hh"
#include "io/vtk.hh"
#include "linalg/direct.hh"
#include "linalg/apcg.hh"                   // Pcg preconditioned CG method
#include "mg/additiveMultigrid.hh"          // makeBPX MG preconditioner
#include "linalg/threadedMatrix.hh"
#include "linalg/triplet.hh"
#include "utilities/enums.hh"
#include "utilities/gridGeneration.hh"
#include "utilities/kaskopt.hh"
#include "utilities/memory.hh"             // for mmoveUnique 
#include "utilities/timing.hh"

using namespace Kaskade;
#include "simple_beam.hh"

To make it easier for the testers to run the computations with different parameters, the program allows users to modify certain parameters through the command line. The first part of main() in simple_beam.cpp will parse the user command and define and initiate all the necessary parameters, including the material definition.

using namespace boost::fusion;


Timings& timer=Timings::instance();
timer.start("total computing time");
std::cout << "Start elastic beam tutorial program." << std::endl;

int coarseGridSize, maxit, order, refinements, skipEntries, skipDimension, solver, verbose, threads,materialList;
bool additive, direct, vtk, linearStrain;
double atol, L, W;
std::string material, storageScheme("A");

if (getKaskadeOptions(argc,argv,Options
("L",                L,                3.,          "length of the elastic beam")
("W",                W,                0.2,        "width/ of the elastic beam")
("refinements",      refinements,      2,          "number of uniform grid refinements")
("order",            order,            1,          "finite element ansatz order")
("material",         material,         "steel",    "type of material")
("linearStrain",     linearStrain,     true,       "true= linear strain tensor, false=nonlinear strain tensor")
("direct",           direct,           true,       "if true, use a direct solver")
("solver",           solver,           0,          "0=UMFPACK, 1=PARDISO 2=MUMPS 3=SUPERLU 4=UMFPACK32/64 5=UMFPACK64")
("verbosity",        verbose,          0,          "amount of reported details")
("vtk",              vtk,              true,       "write solution to VTK file")
("atol",             atol,             1e-8,       "absolute energy error tolerance for iterative solver")
("maxit",            maxit,            100,        "maximum number of iterations")
("MaterialList",     materialList,      1,         "Report the known materials in kaskade")
))
return 1;

if(materialList>0){
auto materials = Elastomechanics::ElasticModulus::materials();
printMaterial(materials);
}

To facilitate easy examination of materials supported by Kaskade 7, running the program without modifying the materialList parameter will cause the program to print a list of materials with their important parameters.

Grid and variables setup

We will construct a grid reflecting the beam using the previously defined width and length, then allow the built-in functionality Gridmanager to refine our grid according to our user-adjustable refinement settings.

std::cout << "refinements of original mesh : " << refinements << std::endl;
std::cout << "discretization order         : " << order << std::endl;

constexpr int DIM = 3;
using Grid = Dune::UGGrid<DIM>;
using Spaces = boost::fusion::vector<H1Space<Grid> const*>;

Dune::FieldVector<double,DIM> x0(0.0), dx(W); dx[0] = L; 
Dune::FieldVector<double,DIM> dh(0.2); dh[0] = 0.6;
GridManager<Grid> gridManager( createCuboid<Grid>(x0, dx, dh, true) );
gridManager.globalRefine(refinements);
gridManager.enforceConcurrentReads(true);

The variable \(u\) will be defined and represent the point-wise displacement caused by the deformation in the beam. The goal of the computation is to find the target displacement function at equilibrium state.

// construction of finite element space for the vector solution u.
H1Space<Grid> h1Space(gridManager,gridManager.grid().leafGridView(),order);

auto varSetDesc = makeVariableSetDescription(makeSpaceList(&h1Space),
                                            make_vector(Variable<SpaceIndex<0>,Components<DIM>>("u")));
using VarSetDesc = decltype(varSetDesc);

Functional definition and obtaining solution

To impose the differential equations governing the displacement variable, we will define a functional \(F\), the type of which is implemented in simple_beam.hh. In order to construct such functional variable, it is necessary to specify the material using ElasticModulus::material(), together with the name of the material given in command line, which is steel by default when it is not specified in command line.

using StrainTensor = LinearizedGreenLagrangeTensor<double,DIM>;
if(!linearStrain)
using StrainTensor = GreenLagrangeTensor<double,DIM>;

using Functional = ElasticityFunctional<VarSetDesc,StrainTensor>;
using Assembler = VariationalFunctionalAssembler<LinearizationAt<Functional> >;
using CoefficientVectors = VarSetDesc::CoefficientVectorRepresentation<0,1>::type;

// Create the variational functional.
Functional F(ElasticModulus::material(material));

A system of linear equations are then constructed, representing the discretized conditions from the functional. Matrix \(A\) (stored in the Assembler variable) and the rhs are assembled in preparation of the processing by a solver.

// construct Galerkin representation
Assembler assembler(varSetDesc.spaces);
VarSetDesc::VariableSet x(varSetDesc);

assembler.assemble(linearization(F,x));

CoefficientVectors rhs(assembler.rhs());
CoefficientVectors solution = varSetDesc.zeroCoefficientVector();

By default, the program will opt for a direct solver (UMFPACK).

if (direct)
{
DirectType directType = static_cast<DirectType>(solver);
AssembledGalerkinOperator<Assembler> A(assembler, directType == DirectType::MUMPS || directType == DirectType::PARDISO);
directInverseOperator(A,directType,MatrixProperties::POSITIVEDEFINITE).applyscaleadd(-1.0,rhs,solution);
component<0>(x) = component<0>(solution);
}
else
{

using X = Dune::BlockVector<Dune::FieldVector<double,DIM>>;
DefaultDualPairing<X,X> dp;
using Matrix = NumaBCRSMatrix<Dune::FieldMatrix<double,DIM,DIM>>;
using LinOp = Dune::MatrixAdapter<Matrix,X,X>;
Matrix Amat(assembler.get<0,0>(),true);
LinOp A(Amat);
SymmetricLinearOperatorWrapper<X,X> sa(A,dp);
PCGEnergyErrorTerminationCriterion<double> term(atol,maxit);

Dune::InverseOperatorResult res;
X xi(component<0>(rhs).N());

std::unique_ptr<SymmetricPreconditioner<X,X>> mg;

if (order==1)
    mg = moveUnique(makeBPX(Amat,gridManager));
else
{
    H1Space<Grid> p1Space(gridManager,gridManager.grid().leafGridView(),1);
    if (storageScheme=="A")
    mg = moveUnique(makePBPX(Amat,h1Space,p1Space,DenseInverseStorageTag<double>(),gridManager.grid().maxLevel()));
    else if (storageScheme=="L")
    mg = moveUnique(makePBPX(Amat,h1Space,p1Space,DenseCholeskyStorageTag<double>(),gridManager.grid().maxLevel()));
    else
    {
    std::cerr << "unknown storage scheme provided\n";
    return -1;
    }
}

Pcg<X,X> pcg(sa,*mg,term,verbose);
pcg.apply(xi,component<0>(rhs),res);
std::cout << "PCG iterations: " << res.iterations << "\n";
xi *= -1;
component<0>(x) = xi;
}

We can then deduce normal stress values from the displacement information. We will then combine both displacement and normal stress values in variable set vsd to be written in VTK form.

//Postprocessing 
L2Space<Grid> l2Space(gridManager,gridManager.grid().leafGridView(),0);

L2Space<Grid>::Element_t<DIM> normalStress(l2Space);
interpolateGlobally<PlainAverage>(normalStress,makeFunctionView(h1Space, [&] (auto const& evaluator)
{
    auto modulus = ElasticModulus::material(material);
    HyperelasticVariationalFunctional<Elastomechanics::MaterialLaws::StVenantKirchhoff<DIM>,StrainTensor>
    energy(modulus);

    energy.setLinearizationPoint(component<0>(x).derivative(evaluator));
    auto stress = Dune::asVector(energy.cauchyStress());

    return Dune::FieldVector<double,3>{stress[0],stress[4],stress[8]};
}));

auto vsd = makeVariableSetDescription(makeSpaceList(&l2Space, &h1Space),
                                    boost::fusion::make_vector(Variable<SpaceIndex<0>,Components<3>>("NormalStress"),
                                                                Variable<SpaceIndex<1>,Components<3>>("Dispacement")));
auto data = vsd.variableSet();
component<0>(data) = normalStress;
component<1>(data) = component<0>(x);
if (vtk)
{
    ScopedTimingSection ts("computing time for file i/o",timer);
    std::string outStr("elasto_");
    if(linearStrain)
        outStr.append("linear_p=");
    else
        outStr.append("nonlinear_p=");
    writeVTK(data,outStr+paddedString(order,1),IoOptions().setOrder(order).setPrecision(16));
}

Visualization

We can visualize the displacement caused by the deformation using ParaView by importing the corresponding VTK file. One should note that to visualize the displacement as color gradient, ‘Displacement’ must be selected from the dropdown menu, where ‘Solid Color’ is the default selection. It is also possible to view normal stress value from the file written.

Displacement

Also, by default the color gradient would represent the magnitude of the displacement vector. One can opt to visualize a particular component instead.

Choose component

Overview of the header file

The definition of the ElasticityFunctional class which contains essentially the equations governing the properties of the materials, as well as the boundary conditions for our problem, is written in header file simple_beam.hh. We note that for the domain cache, most of the details can be completed by calling the function from HyperelasticVariationalFunctional.

#include <type_traits>
#include "fem/functional_aux.hh"
#include "fem/diffops/elastoVariationalFunctionals.hh"

using namespace boost::fusion;

// Deriving from FunctionalBase introduces default D1 and D2 structures.
template <class VarSet, class StrainTensor>
class ElasticityFunctional: public Kaskade::FunctionalBase<VariationalFunctional>
{
public:
  using Scalar = double;
  using AnsatzVars = VarSet;
  using TestVars = VarSet;
  using OriginVars = VarSet;
  static int const dim = AnsatzVars::Grid::dimension;
  using Vector = Dune::FieldVector<Scalar, dim>;
  using Matrix = Dune::FieldMatrix<Scalar, dim, dim>;
  static int constexpr u_Idx = 0;
  static int constexpr u_Space_Idx = result_of::value_at_c<typename AnsatzVars::Variables, u_Idx>::type::spaceIndex;
  
  using MaterialLaw = MaterialLaws::StVenantKirchhoff<dim,Scalar>;
  using ElasticEnergy = HyperelasticVariationalFunctional<MaterialLaw,StrainTensor>;

  class DomainCache 
  {
  public:
    DomainCache(ElasticityFunctional const& functional, typename AnsatzVars::VariableSet const& vars_, int flags=7)
    : vars(vars_), energy(functional.moduli)
    {}

    template <class Entity>
    void moveTo(Entity const& entity) {}

    template <class Position, class Evaluators>
    void evaluateAt(Position const& x, Evaluators const& evaluators)
    {
      using namespace boost::fusion;
      energy.setLinearizationPoint( at_c<u_Idx>(vars.data).derivative(at_c<u_Space_Idx>(evaluators)) );
    }

    Scalar d0() const
    {
      return energy.d0();
    }

    template<int row>
    Vector d1 (VariationalArg<Scalar,dim> const& arg) const
    {
      return energy.d1(arg);
    }

    template<int row, int col>
    Matrix d2 (VariationalArg<Scalar,dim> const& argTest, VariationalArg<Scalar,dim> const& argAnsatz) const
    {
      return energy.d2(argTest,argAnsatz);
    }

  private:
    typename AnsatzVars::VariableSet const& vars;
    ElasticEnergy energy;
  };

  class BoundaryCache : public CacheBase<ElasticityFunctional,BoundaryCache>
  {
  public:
    using FaceIterator = typename AnsatzVars::Grid::LeafIntersectionIterator;

    BoundaryCache(ElasticityFunctional const& f_, typename AnsatzVars::VariableSet const& vars_, int flags=7)
    : vars(vars_)
    {}

    void moveTo(FaceIterator const& face)
    {
      faceIt = &face;
    }

    template <class Evaluators>
    void evaluateAt(Dune::FieldVector<typename AnsatzVars::Grid::ctype,dim-1> const& x, Evaluators const& evaluators)
    {
      u0 = 0;
      u = at_c<u_Idx>(vars.data).value(at_c<u_Space_Idx>(evaluators));
      auto xglob  = (*faceIt)->geometry().global(x);
      
      Vector left(0.); left[0] = -1.;                   // unit left pointing vector
      Vector up(0.); up[dim-1] = 1.;                      // unit upward pointing vector
      auto n = (*faceIt)->centerUnitOuterNormal();      // unit outer normal of local face
      
      Vector force(0.); force[dim-1] = -1.e7;
      
      if (n*left > 0.9)           // clamp beam on the left
      {
        alpha = 1e14;
        beta  = 0.;
      }
      else if( n*up > 0.9 )       // apply downward force on top
      {
        alpha = 0.;
        beta  = force;
      }
      else                        // homogeneous Neumann BCs on the remaining boundary
      {
        alpha = 0.;
        beta  = 0.;
      }
    }

    Scalar
    d0() const
    {
      return alpha*((u-u0)*(u-u0))/2 - beta*u;
    }

    template<int row>
    Scalar d1_impl (VariationalArg<Scalar,dim,dim> const& arg) const
    {
      return alpha*((u-u0)*arg.value) - beta*arg.value;
    }

    template<int row, int col>
    Scalar d2_impl (VariationalArg<Scalar,dim,dim> const &arg1, VariationalArg<Scalar,dim,dim> const &arg2) const
    {
      return alpha*(arg1.value*arg2.value);
    }

  private:
    typename AnsatzVars::VariableSet const& vars;
    Vector u, u0, beta;
    Scalar alpha;
    FaceIterator const* faceIt;
  };

  ElasticityFunctional(Kaskade::ElasticModulus const& moduli_): moduli(moduli_)
  {
  }


  template <class Cell>
  int integrationOrder(Cell const& /* cell */, int shapeFunctionOrder, bool boundary) const
  {
    if (boundary)
      return 2*shapeFunctionOrder;      // mass term u*u on boundary
    else
      return 2*(shapeFunctionOrder-1);  // energy term "u_x * u_x" in interior
  }

  Kaskade::ElasticModulus moduli;
};

#endif /* ELASTOMECHANICS_HH_ */

For our boundary cache, we implemented the d0, d1 and d1 computations using the penalty parameters alpha and beta, which has either zero or a constant non-zero value, depending on where the position of interest is. At a boundary point to the left of the beam, alpha will have nonzero value, activating the dirichlet boundary condition. On the other hand, at a boundary point close to the top of the beam, beta will have a nonzero value, activating the condition corresponding to the downward force applied on the top surface of the beam. Such implementations were carried out in the evaluateAt() function.

template <class Evaluators>
void evaluateAt(Dune::FieldVector<typename AnsatzVars::Grid::ctype,dim-1> const& x, Evaluators const& evaluators)
{
    u0 = 0;
    u = at_c<u_Idx>(vars.data).value(at_c<u_Space_Idx>(evaluators));
    auto xglob  = (*faceIt)->geometry().global(x);
    
    Vector left(0.); left[0] = -1.;                   // unit left pointing vector
    Vector up(0.); up[dim-1] = 1.;                      // unit upward pointing vector
    auto n = (*faceIt)->centerUnitOuterNormal();      // unit outer normal of local face
    
    Vector force(0.); force[dim-1] = -1.e7;
    
    if (n*left > 0.9)           // clamp beam on the left
    {
    alpha = 1e14;
    beta  = 0.;
    }
    else if( n*up > 0.9 )       // apply downward force on top
    {
    alpha = 0.;
    beta  = force;
    }
    else                        // homogeneous Neumann BCs on the remaining boundary
    {
    alpha = 0.;
    beta  = 0.;
    }
}
Page last modified: 2022-01-26 15:49:59 +0100 CET