1. Computational Solid Mechanics Toolbox

1.1. Notations

Notation

Quantity

Unit

\(\rho_s\)

density

\(kg.m^{-3}\)

\(\lambda_s\)

first Lamé coefficients

\(N.m^{-2}\)

\(\mu_s\)

second Lamé coefficients

\(N.m^{-2}\)

\(E_s\)

Young modulus

\(kg.m^{-1}.s^{-2}\)

\(\nu_s\)

Poisson’s ratio

dimensionless

  • strain tensor \(\boldsymbol{F}_s = \boldsymbol{I} + \nabla \boldsymbol{\eta}_s\)

  • Cauchy-Green tensor \(\boldsymbol{C}_s = \boldsymbol{F}_s^{T} \boldsymbol{F}_s\)

  • Green-Lagrange tensor

\begin{align} \boldsymbol{E}_s &= \frac{1}{2} \left( \boldsymbol{C}_s - \boldsymbol{I} \right) \\ &= \underbrace{\frac{1}{2} \left( \nabla \boldsymbol{\eta}_s + \left(\nabla \boldsymbol{\eta}_s\right)^{T} \right)}_{\boldsymbol{\epsilon}_s} + \underbrace{\frac{1}{2} \left(\left(\nabla \boldsymbol{\eta}_s\right)^{T} \nabla \boldsymbol{\eta}_s \right)}_{\boldsymbol{\gamma}_s} \end{align}

1.2. Models

\[ \rho^*_{s} \frac{\partial^2 \boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot \left(\boldsymbol{F}_s \boldsymbol{\Sigma}_s\right) = \boldsymbol{f}^t_s\]

1.2.1. Linear elasticity

\[\begin{align} \boldsymbol{F}_s &= \text{Identity} \\ \boldsymbol{\Sigma}_s &=\lambda_s tr( \boldsymbol{\epsilon}_s)\boldsymbol{I} + 2\mu_s\boldsymbol{\epsilon}_s \end{align}\]

1.2.2. Hyper-elasticity

Saint-Venant-Kirchhoff
\[\boldsymbol{\Sigma}_s=\lambda_s tr( \boldsymbol{E}_s)\boldsymbol{I} + 2\mu_s\boldsymbol{E}_s\]
Neo-Hookean
\[\boldsymbol{\Sigma}_s= \mu_s J^{-2/3}(\boldsymbol{I} - \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{-1})\]
\[\boldsymbol{\Sigma}_s^ = \boldsymbol{\Sigma}_s^\text{iso} + \boldsymbol{\Sigma}_s^\text{vol}\]
Isochoric part : \(\)\boldsymbol{\Sigma}_s^\text{iso}\(\)
Table 1. Isochoric law
Name \(\mathcal{W}_S(J_s)\) \(\boldsymbol{\Sigma}_s^{\text{iso}}\)

Neo-Hookean

\(\mu_s J^{-2/3}(\boldsymbol{I} - \frac{1}{3} \text{tr}(\boldsymbol{C}) \ \boldsymbol{C}^{-1}) \)

Volumetric part : \(\boldsymbol{\Sigma}_s^\text{vol}\)
Table 2. Volumetric law
Name \(\mathcal{W}_S(J_s)\) \(\boldsymbol{\Sigma}_s^\text{vol}\)

classic

\(\frac{\kappa}{2} \left( J_s - 1 \right)^2\)

simo1985

\(\frac{\kappa}{2} \left( ln(J_s) \right)\)

The solid mechanics model can be selected in json file :

Listing : select solid model
"Model":"Hyper-Elasticity"
Table 3. Table of Models for model option
Model Name in json

Linear Elasticity

Elasticity

Hyper Elasticity

Hyper-Elasticity

When materials are closed to incompressibility formulation in displacement/pressure are available.

Table 4. Table of Models for material_law with hyper elasticity model
Model Name Volumic law

Saint-Venant-Kirchhoff

SaintVenantKirchhoff

classic, simo1985

NeoHookean

NeoHookean

classic, simo1985

option: mechanicalproperties.compressible.volumic_law

1.3. Materials

Materials section
"Materials":
{
    "<marker>":
    {
        "name":"solid",
        "E":"1.4e6",
        "nu":"0.4",
        "rho":"1e3"
    }
}

where E stands for the Young’s modulus in \(\)Pa\(\), nu the Poisson’s ratio ( dimensionless ) and rho the density in \(\)kg\cdot m^{-3}\(\).

1.4. Boundary Conditions

Table 5. Boundary conditions
Name Options Type

Dirichlet

faces, edges and component-wise

"Dirichlet"

Neumann

scalar, vectorial

"Neumann_scalar" or "Neumann_vectorial"

Pressure follower ,

Nonlinear boundary condition set in deformed domain

TODO

Robin

TODO

TODO

1.5. Volumic forces

Table 6. Volumic forces
Name Options Type

Expression

Vectorial

"VolumicForces"

1.6. Post Process

1.6.1. Exports for visualisation

The fields allowed to be exported in the Fields section are:

  • displacement

  • velocity

  • acceleration

  • stress or normal-stress

  • pressure

  • material-properties

  • pid

  • fsi

  • Von-Mises

  • Tresca

  • principal-stresses

  • all

1.6.2. Measures

  • Points

  • Maximum

  • Minimum

  • VolumeVariation

Points

Same syntax as FluidMechanics with available Fields :

  • displacement

  • velocity

  • acceleration

  • pressure

  • principal-stress-0

  • principal-stress-1

  • principal-stress-2

  • sigma_xx, sigma_xy, …​

Maximum/Minimum

The Maximum and minimum can be evaluated and save on .csv file. User need to define (i) <Type> ("Maximum" or "Minimum"), (ii) "<tag>" representing this data in the .csv file, (iii) "<marker>" representing the name of marked entities and (iv) the field where extremum is computed.

"<Type>":
{
    "<tag>":
    {
        "markers":"marker>",
        "fields":["displacement","velocity"]
    }
}
VolumeVariation
"VolumeVariation":<marker>

1.7. Action

Let’s finish with a simple example in order to show how this works and how to use them. We will interest us to the deformation of an elastic structure.

1.7.1. Feel++ code

First at all, we define our model type with

typedef FeelModels::SolidMechanics< Simplex<FEELPP_DIM,1>,
                                    Lagrange<OrderDisp, Vectorial,Continuous,PointSetFekete> > model_type;

We choose here \(\)\mathbb{P}_1\(\) space for displacement order. This definition allows us to create our fluid model object SM like this

auto SM = model_type::New("solid");

The method New retrieve all data from the configuration and json files, as well build a mesh if need.

SM->isStationary()

will determine if our model is stationary or not.

If it isn’t, our model is time reliant, and a loop on time is necessary. We then solve our problem and export the results at each time step.

    {
        SM->init();
        SM->printAndSaveInfo();

        for ( ; !SM->timeStepBase()->isFinished(); SM->updateTimeStep() )
        {
            SM->solve();
            SM->exportResults();
        }
    }

If it is stationary, we need to check if we are in quasi static mode or not.

bool algoQuasiStatic = boption(_name="solve-quasi-static");

If not, we save and print our model and solvers. Then the system is solve and we can export the results.

if ( !algoQuasiStatic )
        {
            SM->init();
            SM->printAndSaveInfo();
            SM->solve();
            SM->exportResults();
        }
Code

Here is the code

{% include "../Examples/solid_model.cpp" %}

1.7.2. Configuration file

The config file is used to define options linked to our case we would have the possibility to change at will. It can be, for example, files paths as follows

[solid]
filename=$top_srcdir/applications/models/solid/TurekHron/csm3.json

# precondtioner config
geofile=$top_srcdir/applications/models/solid/TurekHron/csm.geo

[exporter]
directory=applications/models/solid/TurekHron/csm3/$solid_tag

It can also be resolution dependent parameters such as mesh elements size, methods used to define our problem and solvers.

[solid]

material_law=StVenantKirchhoff# StVenantKirchhoff, NeoHookean

# use density and material coeff cst in appli
jacobian-linear-update=false
linearsystem-cst-update=false

# snes and ksp config
#reuse-prec=true#false
#reuse-jac=true#false
reuse-jac.rebuild-at-first-newton-step=true
reuse-prec.rebuild-at-first-newton-step=true
snes-maxit=500
snes-maxit-reuse=10
snes-ksp-maxit=1000
snes-ksp-maxit-reuse=100

# precondtioner config
pc-type=lu #lu,gasm,ml
ksp-converged-reason=1

In this case, we use the Saint-Venant-Kirchhoff model to define our problem, we set the update of linear system constant and jacobian linear as "no update", we discretize values associated to solvers ( SNES and KSP ), and finally we choose LU as the preconditioner method.

Code
{% include "../Examples/csm3.cfg" %}

1.7.3. Json file

First at all, we define some general information like the name ( and short name ) and the model we would like to use

"Name": "Solid Mechanics ",
"ShortName":"Solid",
"Model":"Hyper-Elasticity",

Then we define parameters we will need to solve our problem. Here we define a gravitational constant.

"Parameters":
    {
        "gravity":
        {
            "value":"2"
        }
    },

After that, we define the material properties. In our case, we define the solid we will study, named beam here, by \(\)E\(\), \(\)\nu\(\) and \(\)\rho\(\), respectively its Young’s modulus ( \(\)kg/ms^2\(\) ), its Poisson’s ratio ( dimensionless ) and its density (in \(\)kg/m^3\(\))

Materials":
    {
        "beam":{
            "name":"solid",
            "E":"1.4e6",
            "nu":"0.4",
            "rho":"1e3"
        }
    },

The boundary conditions are the next aspect we define. Here, we impose on the displacement several conditions :

  • A Dirichlet condition on fixed wall

  • A Neumann condition on free wall

  • A volumic force, represent here by the action of the gravity on the solid.

BoundaryConditions":
    {
        "displacement":
        {
            "Dirichlet":
            {
                "fixed-wall":
                {
                    "expr":"{0,0}"
                }
            },
            "Neumann_scalar":
            {
                "free-wall":
                {
                    "expr":"0"
                }
            },
            "VolumicForces":
            {
                "":
                {
                    "expr":"{0,-gravity*1e3}:gravity"
                }
            }
        }
    },

The post process aspect is the last one to define. We want to export displacement values as well as measure displacement and velocity on point \(\)A\(\) along with the maximum of this values on all the free wall.

PostProcess":
    {
        "Fields":["displacement"],
        "Measures":
        {
            "Points":
            {
                "pointA":
                {
                    "coord":"{0.6,0.2,0}",
                    "fields":["displacement","velocity"]
                }
            },
            "Maximum":
            {
                "free-wall":
                {
                    "markers":"free-wall",
                    "fields":["displacement","velocity"]
                }
            }
        }
    }
Code
Unresolved directive in 03-modeling/Solid/README.adoc - include::03-modeling/Examples/csm3.json[]

1.7.4. Compilation/Execution

Once you’ve a build dir, you just have to realise the command make at

{buildir}/applications/models/solid

This will generate executables named feelpp_application_solid_*. To execute it, you need to give the path of the cfg file associated to your case, with --config-file.

For example

mpirun -np 4 feelpp_application_fluid_2d --config-file={sourcedir}/applications/models/solid/TurekHron/csm3.cfg

is how to execute the case ahead on 4 processors.

The result files are then stored by default in

feel/applications/models/solid/{case_name}/{OrderDis}{Geometric_order}/{processor_used}

If we return once again at the example, they are in

feel/applications/models/solid/TurekHron/csm3/P1G1/np_8

2. Computational Solid Mechanics

Computational solid mechanics ( or CSM ) is a part of mechanics that studies solid deformations or motions under applied forces or other parameters.

2.1. Second Newton’s law

The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows

\[\rho^*_s \frac{\partial^2\boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot (\boldsymbol{F}_s\boldsymbol{\Sigma}_s) = \boldsymbol{f}^t_s\]

It’s define here into a Lagrangian frame.

2.1.1. Variables, symbols and units

Notation

Quantity

Unit

\(\rho_s^*\)

strucure density

\(kg/m^3\)

\(\boldsymbol{\eta}_s\)

displacement

\(m\)

\(\boldsymbol{F}_s\)

deformation gradient

\(\boldsymbol{\Sigma}_s\)

second Piola-Kirchhoff tensor

\(N/m^2\)

\(f_s^t\)

body force

\(N/m^2\)

2.2. Lamé coefficients

The Lamé coefficients are deducing from the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material we work on and can be express

\[\lambda_s = \frac{E_s\nu_s}{(1+\nu_s)(1-2\nu_s)} \hspace{0.5 cm} , \hspace{0.5 cm} \mu_s = \frac{E_s}{2(1+\nu_s)}\]

2.3. Axisymmetric reduced model

We interest us here to a 1D reduced model, named generalized string.

The axisymmetric form, which will interest us here, is a tube of length \(L\) and radius \(R_0\). It is oriented following the axis \(z\) and \(r\) represent the radial axis. The reduced domain, named \(\Omega_s^*\) is represented by the dotted line. So the domain, where radial displacement \(\eta_s\) is calculated, is \(\Omega_s^*=\lbrack0,L\rbrack\).

We introduce then \(\Omega_s^{'*}\), where we also need to estimate a radial displacement as before. The unique variance is this displacement direction.

Reduced Model Geometry
Figure 1 : Geometry of the reduce model

The mathematical problem associated to this reduce model can be describe as

\[ \rho^*_s h \frac{\partial^2 \eta_s}{\partial t^2} - k G_s h \frac{\partial^2 \eta_s}{\partial x^2} + \frac{E_s h}{1-\nu_s^2} \frac{\eta_s}{R_0^2} - \gamma_v \frac{\partial^3 \eta}{\partial x^2 \partial t} = f_s.\]

where \(\eta_s\), the radial displacement that satisfy this equation, \(k\) is the Timoshenko’s correction factor, and \(\gamma_v\) is a viscoelasticity parameter. The material is defined by its density \(\rho_s^*\), its Young’s modulus \(E_s\), its Poisson’s ratio \(\nu_s\) and its shear modulus \(G_s\)

At the end, we take \( \eta_s=0\text{ on }\partial\Omega_s^*\) as a boundary condition, which will fixed the wall to its extremities.

Benchmarks

3. Turek Hron CSM Benchmark

In order to validate our fluid-structure interaction solver, we realize here a benchmark on the deformation of an elastic structure, initially proposed by Turek and Hron.

Computer codes, used for the acquisition of results, are from Vincent Chabannes.

This benchmark is linked to the Turek Hron CFD and Turek Hron FSI benchmarks.

3.1. Problem Description

We consider a solid structure, composed of a hyperelastic bar, bound to one of his extremity \(\Gamma_F^*\) to a rigid stationary circular structure. We denote \(\Gamma_{L}^*=\partial\Omega_s^* \backslash \Gamma_F^*\) the other boundaries. The geometry can be represented as follows

CSM Geometry
Figure 1 : Geometry of the Turek & Hron CSM benchmark.

\(\Omega_s^*\) represent the initial domain, before any deformations. We denote \(\Omega^t_s\) the domain obtained during the deformations application, at the time t.

By the Newton’s second law, we can describe the fundamental equation of the solid mechanic in a Lagrangian frame. Furthermore, we will suppose that the hyperelastic material follows a compressible Saint-Venant-Kirchhoff model. The Lamé coefficients, used in this model, are deducing from the Young’s modulus \(E_s\) and the Poisson’s ratio \(\nu_s\) of the material.

We will observe, during this simulation, the displacement of \(A\), on the \(x\) and \(y\) axis, when the elastic structure is subjected to its own weight, and compare our results to the reference ones given Turek and Hron.

In the first two cases, CSM1 and CSM2, we want to determine the steady state condition. To find it, a quasi-static algorithm is used, which increase at each step the gravity parameter.

For the third cases, CSM3, we realise a transient simulation, where we will observe the comportment of \(A\) subjected to its weight during a given time span.

3.1.1. Boundary conditions

We set

  • on \(\Gamma_{F}^*\), a condition that imposes this boundary to be fixed : \(\boldsymbol{\eta}_s=0\)

  • on \(\Gamma_{L}^*\), a condition that lets these boundaries be free from constraints : \((\boldsymbol{F}_s\boldsymbol{\Sigma}_s)\boldsymbol{ n }^*_s=\boldsymbol{0}\)

where \(\boldsymbol{n}_s^*\) is the outer unit normal vector from \(\partial \Omega_s^*\).

3.1.2. Initial conditions

The source term \(\boldsymbol{f}_s\), chosen in order to set in motion the elastic structure, is define by

\[\boldsymbol{f}_s=(0,-\rho^*_s g)\]

where \(g\) is a gravitational constant.

The reference document ( Turek and Hron ) don’t specify the time interval used to obtain their results. In our case, it’s chosen as \(\lbrack0,10\rbrack\).

3.2. Inputs

The following table displays the various fixed and variables parameters of this test-case.

Table 7. Fixed and Variable Input Parameters
Name Description Nominal Value Units

\(g\)

gravitational constant

2

\(m / s^2\)

\(l\)

elastic structure length

\(0.35\)

\(m\)

\(h\)

elastic structure height

\(0.02\)

\(m\)

\(r\)

cylinder radius

\(0.05\)

\(m\)

\(C\)

cylinder center coordinates

\((0.2,0.2)\)

\(m\)

\(A\)

control point coordinates

\((0.2,0.2)\)

\(m\)

\(B\)

point coordinates

\((0.15,0.2)\)

\(m\)

\(E_s\)

Young’s modulus

\(1.4\times 10^6\)

\(kg / ms^2\)

\(\nu_s\)

Poisson’s ratio

\(0.4\)

dimensionless

\(\rho^*_s\)

density

\(1000\)

\(kg/ m^3\)

As for solvers we used, Newton’s method is chosen for the non-linear part and a direct method based on LU decomposition is selected for the linear part.

3.3. Outputs

As described before, we have

\[\rho^*_s \frac{\partial^2\boldsymbol{\eta}_s}{\partial t^2} - \nabla \cdot (\boldsymbol{F}_s\boldsymbol{\Sigma}_s) = \boldsymbol{f}^t_s\]

We search the displacement \(\boldsymbol{\eta}_s\), on \(\Omega_s^*\), which will satisfy this equation.

In particular, the displacement of the point \(A\) is the one that interests us.

3.4. Discretization

To realize these tests, we made the choice to used Finite Elements Method, with Lagrangian elements of order \(N\) to discretize space.

Newmark-beta method, presented into Chabannes papers, is the one we used for the time discretization. We used this method with \(\gamma=0.5\) and \(\beta=0.25\).

3.4.1. Solvers

Here are the different solvers ( linear and non-linear ) used during results acquisition.

Table 8. KSP configuration

type

gmres

relative tolerance

1e-13

max iteration

1000

reuse preconditioner

true

Table 9. SNES configuration

relative tolerance

1e-8

steps tolerance

1e-8

max iteration

500

max iteration with reuse

10

reuse jacobian

false

reuse jacobian rebuild at first Newton step

true

Table 10. KSP in SNES configuration

relative tolerance

1e-5

max iteration

500

reuse preconditioner

CSM1/CSM2 : false | CSM3 : true

reuse preconditioner rebuild at first Newton step

true

Table 11. Preconditioner configuration

type

lu

package

mumps

3.5. Implementation

To realize the acquisition of the benchmark results, code files contained and using the Feel++ library will be used. Here is a quick look to the different location of them.

First at all, the main code can be found in

    feelpp/applications/models/solid

The configuration file for the CSM3 case, the only one we work on, is located at

    feelpp/applications/models/solid/TurekHron

The result files are then stored by default in

    feel/applications/models/solid/TurekHron/csm3/"OrderDisp""Geometric_order"/"processor_used"

Like that, for the CSM3 case executed on 8 processors, with a \(P_1\) displacement approximation space and a geometric order of 1, the path is

    feel/applications/models/solid/TurekHron/csm3/P1G1/np_8

At least, to retrieve results that interested us for the benchmark and to generate graphs, we use a Python script located at

    feelpp-benchmarking-book/CFD/Turek-Hron/postprocess_cfd.py

3.6. Results

3.6.1. CSM1

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

Reference TurekHron

-7.187

-66.10

1061

4620 (\(P_2\))

-7.039

-65.32

4199

17540 (\(P_2\))

-7.047

-65.37

16495

67464 (\(P_2\))

-7.048

-65.37

1061

10112 (\(P_3\))

-7.046

-65.36

1906

17900 (\(P_3\))

-7.049

-65.37

1061

17726 (\(P_4\))

-7.048

-65.37

All the files used for this case can be found in this rep [ geo file, config file, json file ]

3.6.2. CSM2

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

Reference TurekHron

-0.4690

-16.97

1061

4620 (\(P_2\))

-0.459

-16.77

4201

17548 (\(P_2\))

-0.459

-16.77

16495

67464 (\(P_2\))

-0.459

-16.78

1061

10112 (\(P_3\))

-0.4594

-16.78

16475

150500 (\(P_3\))

-0.460

-16.78

1061

17726 (\(P_4\))

-0.460

-16.78

All the files used for this case can be found in this rep [geo file, config file, json file].

3.6.3. CSM3

The results of the CSM3 benchmark are detailed below.

Table 12. Results for CSM3

\(\Delta t\)

\(N_{elt}\)

\(N_{dof}\)

\(x\) displacement \(\lbrack\times 10^{-3}\rbrack\)

\(y\) displacement \(\lbrack\times 10^{-3}\rbrack\)

/

Reference TurekHron

−14.305 ± 14.305 [1.0995]

−63.607 ± 65.160 [1.0995]

0.02

4199

17536(\(P_2\))

-14.585 ± 14.590 [1.0953]

-63.981 ± 65.521 [1.0930]

4199

38900(\(P_3\))

-14.589 ± 14.594 [1.0953]

-63.998 ± 65.522 [1.0930]

1043

17536(\(P_4\))

-14.591 ± 14.596 [1.0953]

-64.009 ± 65.521 [1.0930]

4199

68662(\(P_4\))

-14.590 ± 14.595 [1.0953]

-64.003 ± 65.522 [1.0930]

0.01

4199

17536(\(P_2\))

-14.636 ± 14.640 [1.0969]

-63.937 ± 65.761 [1.0945]

4199

38900(\(P_3\))

-14.642 ± 14.646 [1.0969]

-63.949 ± 65.771 [1.0945]

1043

17536(\(P_4\))

-14.645 ± 14.649 [1.0961]

-63.955 ± 65.778 [1.0945]

4199

68662(\(P_4\))

-14.627 ± 14.629 [1.0947]

-63.916 ± 65.739 [1.0947]

0.005

4199

17536(\(P_2\))

-14.645 ± 14.645 [1.0966]

-64.083 ± 65.521 [1.0951]

4199

38900(\(P_3\))

-14.649 ± 14.650 [1.0966]

-64.092 ± 65.637 [1.0951]

1043

17536(\(P_4\))

-14.652 ± 14.653 [1.0966]

-64.099 ± 65.645 [1.0951]

4199

68662(\(P_4\))

-14.650 ± 14.651 [1.0966]

-64.095 ± 65.640 [1.0951]

fullviewCSM

\(\) \text{Figure 2: x and y displacements} \(\)

All the files used for this case can be found in this rep [ geo file, config file, json file ]

3.6.4. Conclusion

To obtain these data, we used several different mesh refinements and different polynomial approximations for the displacement on the time interval \(\lbrack 0,10 \rbrack\).

Our results are pretty similar to those from Turek and Hron, despite a small gap. This gap can be caused by the difference between our time interval and the one used for the reference acquisitions.

3.7. Bibliography

References for this benchmark
  • [TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow, Lecture Notes in Computational Science and Engineering, 2006.

  • [Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Université de Grenoble, 2013.

4. Solenoïd Benchmark

An analytical solution of the elasticity equation exists in an infinitely long solenoid. Such a geometry has the advantage to be axisymetrical, we can first expresse the solution with cylindrical coordinates.

4.1. Definition

4.1.1. Equilibrium equation in cylindrical coordinates

The analytical solution comes from the electromagnetism equations, indeed this case of a soleinoid crossed by an electric current can model the behaviour of a resistive magnet. The volumic forces \(\mathbf{f}\) of the equilibrium equation are consequently dependent on the current density \(\mathbf{J}\) and the induced magnetic field \(\mathbf{B}\), which gives :

\[\nabla\cdot\bar{\bar{\sigma}} + \mathbf{J}\times\mathbf{B} = \mathbf{0}\]

We will use the following notations:

\[\mathbf{u} = \begin{pmatrix} u \\ v \\ w \\ \end{pmatrix} \quad \mathbf{J} = \begin{pmatrix} J_{r} \\ J_{\theta} \\ J_{z} \\ \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} B_{r} \\ B_{\theta} \\ B_{z} \\ \end{pmatrix}, \quad \text{and} \quad \bar{\bar{\sigma}} = \begin{pmatrix} \sigma_{r r} & \sigma_{r \theta} & \sigma_{r z} \\ \sigma_{\theta r} & \sigma_{\theta \theta} & \sigma_{\theta z} \\ \sigma_{z r} & \sigma_{z \theta} & \sigma_{z z} \\ \end{pmatrix}\]

We consider here that the current distribution in the soleinoid is uniform. That means that \(J_{\theta}\) is constant, and \(J_r = J_z = 0\). The volumic forces \(\mathbf{J} \times \mathbf{B}\) becomes :

\[\mathbf{J}\times\mathbf{B} = \begin{pmatrix} J_{\theta} B_z \\ 0 \\ - J_{\theta} B_r \end{pmatrix}\]

The revwriting of the equilibrium equation in cylindrical coordinates gives :

\[\left\{ \begin{array}{l} \displaystyle{ \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} + \frac{1}{r}( \sigma_{rr} - \sigma_{\theta \theta} ) + J_{\theta} B_z = 0 }\\[0.4cm] \displaystyle{ \frac{\partial \sigma_{\theta r}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{\theta \theta}}{\partial \theta} + \frac{\partial \sigma_{\theta z}}{\partial z} + \frac{2}{r} \sigma_{\theta r} = 0 }\\[0.4cm] \displaystyle{ \frac{\partial \sigma_{zr}}{\partial r} + \frac{1}{r}\frac{\partial \sigma_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{zz}}{\partial z} + \frac{1}{r} \sigma_{rz} - J_{\theta} B_r = 0 }\\ \end{array} \right.\]

The axisymetrical properties of this geometries means that the displacement is invariant with respect to \(\theta\).
Furthermore, the soleinoid we consider has the particularity to be infinitely long, so there is no displacement along \(z\) axis.
We can consequently get rid of all derivatives \(\frac{\partial \cdot}{\partial \theta}\) and \(\frac{\partial \cdot}{\partial z}\).
Finally, we shall note that components \(\sigma_{\theta r}\) and \(\sigma_{zr}\) of the stress tensor \(\bar{\bar{\sigma}}\) are expressed from Hooke’s law only from the components \(v\) and \(w\) of the displacement vector \(\mathbf{u}\), which nullify the two last equations.

Thus, the equilibrium equation is reduced to:

\[-\sigma_{\theta \theta} + \frac{\partial}{\partial r}(r \sigma_{rr}) = -r J_{\theta} B_z\]

We need to express this equation in terms of the displacement \(\mathbf{u}\), This can be done using Hooke’s law which links the stress tensor to the tensor of small deformation by:

\[\bar{\bar{\sigma}}=2\mu\bar{\bar{\varepsilon}} + \lambda Tr(\bar{\bar{\varepsilon}})Id\]

where \(\mu\) and \(\lambda\) are the Lamé coefficients, and the tensor of small deformation is given in cylindrical coordinates by:

\[\bar{\bar{\varepsilon}} = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) = \begin{pmatrix} \frac{\partial u}{\partial r} & \frac{1}{2}\left( \frac{1}{r}\frac{\partial u}{\partial \theta} - \frac{v}{r} + \frac{\partial v}{\partial r} \right) & \frac{1}{2} \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial r} \right) \\ \frac{1}{2}\left( \frac{1}{r}\frac{\partial u}{\partial \theta} - \frac{v}{r} + \frac{\partial v}{\partial r} \right) & \frac{1}{r}\frac{\partial v}{\partial \theta} + \frac{u}{r} & \frac{1}{2} \left( \frac{1}{r}\frac{\partial w}{\partial \theta} + \frac{\partial v}{\partial z} \right) \\ \frac{1}{2} \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial r} \right) & \frac{1}{2} \left( \frac{1}{r}\frac{\partial w}{\partial \theta} + \frac{\partial v}{\partial z} \right) & \frac{\partial w}{\partial z} \end{pmatrix}\]

Then, using the last two definitions and the properties of the solenoid (axisymmetric and infinitely long), we can rewrite the equilibrium equation as:

\[\frac{\partial}{\partial r}\left( r \frac{\partial u}{\partial r} \right) - \frac{u}{r} = - \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r J_{\theta} B_z\]

4.1.2. Analytical solution

We want to find an analytical solution of the form :

\[u_{cyl}(r) = C_1 r + \frac{C_2}{r} + u_p(r)\]

where \(C_1\) and \(C_2\) are constants, and \(u_p(r)\) a particular solution of the equilibrium equation in cylindrical coordinates.

From Ampére’s theorem and considering that the soleinoid is axisymetrical and infinitely long, we deduce that \(B_r = 0\), and \(B_z\) depends only on \(r\), such that :

\[B_z(r_1) - B_z(r) = \frac{1}{\mu} \int_{r_1}^r J_{\theta} dr\]

with \(r_1\) the internal radius of the soleinoid.

For a uniform distribution of current in the solenoid (\(j_{\theta}\) constant), we deduce that \(B_z\) can be expressed as :

\[B_z(r) = B_z(r_1) - \frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1} - 1 \right)\]

Replacing \(b_z\) with his expression in the equilibrium equation, this gives :

\[ \frac{\partial}{\partial r}\left( r \frac{\partial u}{\partial r} \right) - \frac{u}{r} = - \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r J_{\theta} \left( B_z(r_1) - \frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1} - 1 \right) \right)\]

where \(r_2\) is the external radius, \(\alpha = \frac{r_2}{r_1}\) and \(\Delta b_z = B_z(r_1) - B_z(r_2)\).

A particular solution \(u_p(r)\) for this equation is given by:

\[u_p(r) = \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}r_1 J_{\theta} \left[ -\frac{r_1}{3}\left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1} \right) \left( \frac{r}{r_1}\right)^2 + \frac{r_1}{8}\frac{\Delta B_z}{\alpha - 1} \left( \frac{r}{r_1}\right)^3 \right]\]

The constants \(C_1\) and \(C_2\) are set by the boundary conditions, we consider here that there is no surface forces. That gives \(\bar{\bar{\sigma}}\cdot\mathbf{n} = 0\) on internal and external radius, that is \(\sigma_{rr}(r_1)=\sigma_{rr}(r_2)=0\).

Using the definition of \(u_{cyl}\) in

\[\sigma_{rr} = \frac{E}{(1+\nu)(1-2\nu)} \left[ (1-\nu)\frac{\partial u}{\partial r} + \nu \frac{u}{r} \right]\]

we can solve the system to find the constants:

\[\begin{align} C_1 = \frac{(1+\nu)(1-2\nu)}{E(1-\nu)}J_{\theta}r_1 &\left[ \left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1}\right) \left( \frac{2 - \nu}{3} \right) \left( 1 - \frac{r_2^2}{(r_2^2 - r_1^2)}\left( 1 - \frac{r_2}{r_1} \right) \right) \right. \\ &+ \left. \frac{\Delta B_z}{\alpha - 1}\left( \frac{2\nu - 3}{8} \right) \left( 1 - \frac{r_2^2}{(r_2^2 - r_1^2)} \left( 1 - \left( \frac{r_2}{r_1} \right)^2 \right) \right) \right] \end{align}\]
\[C_2 = \frac{-r_1^3 r_2^2 J_{\theta}}{(r_2^2 - r_1^2)}\frac{(1+\nu)}{E(1-\nu)} \left[ \left( B_z(r_1) + \frac{\Delta B_z}{\alpha - 1}\right) \left( \frac{2 - \nu}{3} \right)\left( 1 - \frac{r_2}{r_1} \right) +\frac{\Delta B_z}{\alpha - 1}\left( \frac{2\nu - 3}{8} \right)\left( 1 - \left( \frac{r_2}{r_1} \right)^2 \right) \right]\]

The final step is to translate this analytical solution \(u_{cyl}(r)\) into cartesian coordinates to obtain the analytical cartesian displacement \(\mathbf{u}_{cart}\):

\[\mathbf{u}_{cart}(x,y)= \begin{pmatrix} cos( \theta )u_{cyl}(\sqrt{x^2 + y^2})\\ sin( \theta )u_{cyl}(\sqrt{x^2 + y^2})\\ 0 \end{pmatrix} = \begin{pmatrix} \frac{x}{\sqrt{x^2 + y^2}}u_{cyl}(\sqrt{x^2 + y^2}) \\ \frac{y}{\sqrt{x^2 + y^2}}u_{cyl}(\sqrt{x^2 + y^2}) \\ 0 \end{pmatrix}\]

4.1.3. Geometry

We use a solenoïd of thickness one with \(r_1=1\) and \(r_2=2\) and with a length sufficiently important (\(l=10\,r_2\)) so that the influence of the top and of the bottom of the geometry, which are supposed not to exist, is close to zero.

4.1.4. Boundary conditions

The boundary conditions taken into account for the analytical solution have to be reproduced for the simulation. That means null pressure forces on internal and external radius, and displacement set to zero (Dirichlet) on the top and on the bottom to keep only the radial component.

We set:

  • \(\mathbf{u} = 0\) on \(\Gamma_{top}\cup\Gamma_{bottom}\)

  • \(\bar{\bar{\sigma}}\cdot \mathbf{n} = 0\) on \(\Gamma_{int}\cup\Gamma_{ext}\)

  • \(\mathbf{f} = \begin{pmatrix} \frac{x}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ \frac{y}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ 0 \end{pmatrix}\) in \(\Omega\)

4.2. Inputs

We use the following parameters:

Table 13. Inputs
Name Value

\(E\)

\(2.1e^6\)

\(\nu\)

\(0.33\)

\(\mathbf{f}\)

\(\begin{pmatrix} \frac{x}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ \frac{y}{\sqrt{x^2+y^2}}5000(10+20(\sqrt{x^2+y^2}-1))\\ 0 \end{pmatrix}\)

4.3. Output

We compare the radial component of the displacement on the segment \(z=l/2\), \(y=0\) and \(x\in \lbrack 1,2\rbrack \).

4.4. Results

Here are the analytical and the computed \(x\) component of the displacement. This has been obtain with a characteristic size of \(0.1\) and \(646 233\) dofs.

We can see that the errors grows as we approach the external radius. But the max of the error is \(5e^{-4}\) and it converges as the characteristic size decreases.

5. NAFEMS LE1 Benchmarck

This benchmark is extract from the Abaqus Benchmarks Manual.

5.1. Definition

We focus on the LE1 benchmarks in particular.

5.1.1. Geometry

The geometry is given here by :
geoLE10

5.1.2. Boundary conditions

We set:

  • \(u_y = 0\) on DC

  • \(u_x = 0\) on AB

  • \(\bar{\bar{\varepsilon}}\cdot\mathbf{n}=1e^7\) on BC.

5.2. Inputs

We have the following parameters:

Table 14. Inputs
Name Value

\(E\)

\(210\, GPa\)

\(\nu\)

\(0.3\)

\(\rho\)

\(7800\, kg/m^2\)

5.3. Outputs

We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(92.7\, MPa\).

5.4. Results

The value of \(\sigma_{yy}\) at the point D is \(94.09\, MPa\) for \(32 000\) dofs, which is \(1.49%\) higher than the target.

One possibility to get a more accurate output is to use a mixed formulation, where the stress tensor would also be an unknown.

6. NAFEMS LE10 Benchmarck

This benchmark is extract from the Abaqus Benchmarks Manual.

6.1. Definition

We focus on the LE10 benchmarks in particular.

6.1.1. Geometry

The geometry is given here by :

geoLE10

where the thickness is \(0.6\).

In addition, we define the point E which is the midpoint of CC' and E' the midpoint of BB'.

6.1.2. Boundary conditions

We set:

  • \(u_y = 0\) on DCD’C'

  • \(u_x = 0\) on ABA’B'

  • \(u_x = u_y = 0\) on BCB’C'

  • \(\bar{\bar{\varepsilon}}\cdot\mathbf{n}=-1e^6\) on the top surface.

6.2. Inputs

We have the following parameters:

Table 15. Inputs
Name Value

\(E\)

\(210\, GPa\)

\(\nu\)

\(0.3\)

\(\rho\)

\(7800\, kg/m^2\)

6.3. Outputs

We want to compare the value of \(\sigma_{yy}\) at the point D. The reference value is \(5.38\, MPa\).

6.4. Results

The value of \(\sigma_{yy}\) at the point D is \(5.53\, MPa\) for \(300 000\) dofs, which is \(2.78%\) higher than the target.

One possibility to get a more accurate output is to use a mixed formulation, where the stress tensor would also be an unknown.