true

This demo is a direct transposition of the transient
thermo-elasticity demo using a pure `FEniCS`

formulation.
We will show how to compute fully coupled thermo-mechanical problems
using `MFront`

, which can pave the way to more complex
thermo-mechanical behaviours including plasticity for instance.

Source files:

- Jupyter notebook: mgis_fenics_monolithic_transient_thermoelasticity.ipynb
- Python file: mgis_fenics_monolithic_transient_thermoelasticity.py
- MFront behaviour file: ThermoElasticity.mfront

The constitutive equations are derived from the following expression of the Gibbs free energy:

\[ \begin{aligned} \rho\,\Phi{\left(\boldsymbol{\varepsilon}^{\mathrm{to}},T\right)}&={{\displaystyle \frac{\displaystyle \lambda}{\displaystyle 2}}}\,{\left({\mathrm{tr}{\left(\boldsymbol{\varepsilon}^{\mathrm{to}}\right)}}-3\,\alpha\,{\left(T-T^{\mathrm{ref}}\right)}\right)}^2+ \mu\,{\left(\boldsymbol{\varepsilon}^{\mathrm{to}}-\alpha\,{\left(T-T^{\mathrm{ref}}\right)}\,\mathbf{I}\right)}\,\colon\,{\left(\boldsymbol{\varepsilon}^{\mathrm{to}}-\alpha\,{\left(T-T^{\mathrm{ref}}\right)}\,\mathbf{I}\right)}\\ &+{{\displaystyle \frac{\displaystyle \rho\,C_{\varepsilon}}{\displaystyle 2\,T^{\mathrm{ref}}}}}\,{\left(T-T^{\mathrm{ref}}\right)}^2+s_{0}\,{\left(T-T^{\mathrm{ref}}\right)} \end{aligned} \]

where:

- \(\lambda\) and \(\mu\) are the Lamé coefficients
- \(\rho\) is the mass density
- \(\alpha\) is mean linear thermal expansion coefficient
- \(C_{\varepsilon}\) is the specific heat at constant strain (per unit of mass).

This expression leads to the following expressions of the stress tensor \(\boldsymbol{\sigma}\) and entropy per unit of mass \(s\):

\[ \begin{aligned} \boldsymbol{\sigma}&=\rho \dfrac{\partial \Phi}{\partial \boldsymbol{\varepsilon}^{\mathrm{to}}}=\lambda\,{\mathrm{tr}{\left(\boldsymbol{\varepsilon}^{\mathrm{to}}\right)}}\,\mathbf{I}+2\,\mu\,\boldsymbol{\varepsilon}^{\mathrm{to}}-\kappa\,{\left(T-T^{\mathrm{ref}}\right)}\,\mathbf{I}\\ s&={\displaystyle \frac{\displaystyle \partial \Phi}{\displaystyle \partial T}}={{\displaystyle \frac{\displaystyle C_{\varepsilon}}{\displaystyle T^{\mathrm{ref}}}}}\,{\left(T-T^{\mathrm{ref}}\right)}+{{\displaystyle \frac{\displaystyle \kappa}{\displaystyle \rho}}}\,{\mathrm{tr}{\left(\boldsymbol{\varepsilon}^{\mathrm{to}}\right)}}\\ \end{aligned} \qquad{(1)}\]

where \(\kappa=\alpha\,{\left(3\,\lambda+2\,\mu\right)}\).

The heat flux \(\mathbf{j}\) is related to the temperature gradient \(\nabla\, T\) by the linear Fourier law:

\[ \mathbf{j}=-k\,\nabla\, T \qquad{(2)}\]

`MFront`

implementationThe constitutive equations (1) and (2) exhibit an explicit expression of the thermodynamic forces \({\left(\boldsymbol{\sigma}\, \mathbf{j}, s\right)}\) as a function of the gradients \({\left(\boldsymbol{\varepsilon}^{\mathrm{to}}, \nabla T, T\right)}\).

The most suitable domain specific language for this kind of behaviour
if the `DefaultGenericBehaviour`

.

`@DSL DefaultGenericBehaviour;`

The `@Behaviour`

keyword allows giving the name of the
behaviour:

`@Behaviour ThermoElasticity;`

The following lines add some metadata (authors of the implementation, date, description):

```
@Author Thomas Helfer, Jérémy Bleyer;
@Date 19/04/2020;
@Description {
This simple thermoelastic behaviour allows to perform-mechanical resolutions.
fully coupled thermo
://comet-fenics.readthedocs.io/ for details.
See https}
```

The gradients are the strain \(\boldsymbol{\varepsilon}^{\mathrm{to}}\), the temperature gradient \(\nabla\,T\) and the temperature. The associated thermodynamic forces are respectively the stress \(\boldsymbol{\sigma}\), the heat flux \(\mathbf{j}\) and the entropy \(s\).

\(\boldsymbol{\varepsilon}^{\mathrm{to}}\), \(\nabla\,T\), \(\boldsymbol{\sigma}\) and \(\mathbf{j}\) are declared as follows:

```
@Gradient StrainStensor εᵗᵒ;
.setGlossaryName("Strain");
εᵗᵒ
@ThermodynamicForce StressStensor σ;
.setGlossaryName("Stress");
σ
@Gradient TemperatureGradient ∇T;
.setGlossaryName("TemperatureGradient");
∇T
@ThermodynamicForce HeatFlux j;
.setGlossaryName("HeatFlux"); j
```

The glossary names are the names seen from the calling solver. Glossary names are described on this page.

Due to a `MFront`

convention, the temperature is
automatically declared as an external state variable. For this reason,
the entropy is declared as a state variable:

```
@StateVariable real s;
.setEntryName("EntropyPerUnitOfMass"); s
```

In the current version of `MFront`

, there is no glossary
name associated with the entropy per unit of mass. In this case, the
`setEntryName`

is used to associate a name to this
variable.

By default, all the derivatives of the thermodynamic forces with respect to the increments of the gradients are declared as tangent operator blocks, i.e. derivatives that are meant to be used when building the stiffness matrix at the structural scale.

In this case, this is not appropriate as:

- some derivatives are known to be null, such as \({\displaystyle \frac{\displaystyle \partial \boldsymbol{\sigma}}{\displaystyle \partial \Delta\,\nabla\,T}}\) and \({\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial \Delta\,\boldsymbol{\varepsilon}^{\mathrm{to}}}}\).
- the derivative \({\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,\boldsymbol{\varepsilon}^{\mathrm{to}}}}\) of the entropy with respect to strain and the derivative \({\displaystyle \frac{\displaystyle \partial s}{\displaystyle \partial \Delta\,T}}\) of the entropy with respect to the temperature are also required.

The required tangent operator blocks are therefore explicitly requested:

`@TangentOperatorBlocks{∂σ∕∂Δεᵗᵒ, ∂σ∕∂ΔT, ∂s∕∂ΔT, ∂s∕∂Δεᵗᵒ, ∂j∕∂Δ∇T};`

The reference temperature is declared using the
`@StaticVariable`

keyword:

`@StaticVariable temperature Tʳᵉᶠ = 293.15;`

Internally `Tʳᵉᶠ`

is hold in an immutable static
variable.

The various material coefficients are now declared as parameters:

```
@Parameter stress E = 70e3;
.setGlossaryName("YoungModulus");
E@Parameter real ν = 0.3;
.setGlossaryName("PoissonRatio");
ν@Parameter massdensity ρ = 2700.;
.setGlossaryName("MassDensity");
ρ@Parameter thermalconductivity α = 2.31e-5 ;
.setGlossaryName("ThermalExpansion");
α@Parameter real Cₑ = 910e-6;
.setEntryName("SpecificHeatAtConstantStrainPerUnitOfMass");
Cₑ@Parameter thermalconductivity k = 237e-6;
.setGlossaryName("ThermalConductivity"); k
```

Parameters are global values that can be modified at runtime.

The computation of the thermodynamic forces and tangent operator
blocks is implemented in the `@Integrator`

code block:

`@Integrator{`

First, the Lamé coefficients are computed using the built-in
`computeLambda`

and `computeMu`

functions and then
we compute the \(\kappa\) factor:

```
const auto λ = computeLambda(E, ν);
const auto μ = computeMu(E, ν);
const auto κ = α ⋅ (2 ⋅ μ + 3 ⋅ λ);
```

For brevity, we compute the strain at the end of the time step as follows:

`const auto ε = εᵗᵒ + Δεᵗᵒ; `

The computation of the thermodynamic forces is then straightforward and closely looks like the constitutive equations (1) and (2):

```
= λ ⋅ trace(ε) ⋅ I₂ + 2 ⋅ μ ⋅ ε - κ ⋅ (T + ΔT - Tʳᵉᶠ) ⋅ I₂;
σ = Cₑ / Tʳᵉᶠ ⋅ (T + ΔT - Tʳᵉᶠ) + (κ / ρ) ⋅ trace(ε);
s = -k ⋅ (∇T + Δ∇T); j
```

The computation of the consistent tangent operator is only required
if the `computeTangentOperator_`

boolean value is true.
Again, their computations is straightforward [2]:

```
if (computeTangentOperator_) {
= λ ⋅ (I₂ ⊗ I₂) + 2 ⋅ μ ⋅ I₄;
∂σ∕∂Δεᵗᵒ = -κ ⋅ I₂;
∂σ∕∂ΔT = Cₑ / Tʳᵉᶠ;
∂s∕∂ΔT = κ ⋅ Cₑ / Tʳᵉᶠ ⋅ I₂;
∂s∕∂Δεᵗᵒ = -k ⋅ tmatrix<N, N, real>::Id();
∂j∕∂Δ∇T }
```

A final curly bracket then ends the `@Integrator`

code
block:

`}`

[1] We may also note that those blocks are third order tensors that
are not yet supported by `MFront`

.

[2] `N`

is the space dimension. `real`

is a
type alias to the numeric type used, which depends on the interface
used.

`FEniCS`

implementationThe problem consists of a quarter of a square plate perforated by a circular hole. A temperature increase of \(\Delta T=+10^{\circ}\text{C}\) will be applied on the hole boundary. Symmetry conditions are applied on the corresponding symmetry planes and stress and flux-free boundary conditions are adopted on the plate outer boundary. Similarly to the original demo, we will formulate the problem using the temperature variation as the main unknown.

We first import the relevant modules then define the mesh and some constants.

```
from dolfin import *
import mgis.fenics as mf
from mshr import Rectangle, Circle, generate_mesh
import numpy as np
import matplotlib.pyplot as plt
= 1.0
L = 0.1
R = 50 # mesh density
N
= Rectangle(Point(0.0, 0.0), Point(L, L)) - \
domain 0.0, 0.0), R, 100)
Circle(Point(= generate_mesh(domain, N)
mesh
= Constant(293.15)
Tref = Constant(10.0)
DThole = Constant(0) # time step dt
```

We now define the relevant FunctionSpace for the considered problem. Since we will adopt a monolithic approach i.e. in which both fields are coupled and solved at the same time, we will need to resort to a Mixed FunctionSpace for both the displacement \(\boldsymbol{u}\) and the temperature variation \(\Theta = T-T^\text{ref}\).

```
= VectorElement("CG", mesh.ufl_cell(), 2) # displacement finite element
Vue = FiniteElement("CG", mesh.ufl_cell(), 1) # temperature finite element
Vte = FunctionSpace(mesh, MixedElement([Vue, Vte]))
V
def inner_boundary(x, on_boundary):
return near(x[0]**2 + x[1]**2, R**2, 1e-3) and on_boundary
def bottom(x, on_boundary):
return near(x[1], 0) and on_boundary
def left(x, on_boundary):
return near(x[0], 0) and on_boundary
= [DirichletBC(V.sub(0).sub(1), Constant(0.0), bottom),
bcs 0).sub(0), Constant(0.0), left),
DirichletBC(V.sub(1), DThole, inner_boundary)] DirichletBC(V.sub(
```

The constitutive equations described earlier are completed by the quasi-static equilibrium equation:

\[ \text{div} \boldsymbol{\sigma}= 0 \]

and the transient heat equation (without source terms):

\[ \rho T^\text{ref} \dfrac{\partial s}{\partial t} + \text{div} \mathbf{j}= 0 \]

which can both be written in the following weak form:

\[ \begin{aligned} \int_{\Omega}\boldsymbol{\sigma}:\nabla^s\widehat{\boldsymbol{u}}\text{ d} \Omega &=\int_{\partial \Omega} (\boldsymbol{\sigma}\cdot\boldsymbol{n})\cdot\widehat{\boldsymbol{u}} dS \quad \forall \widehat{\boldsymbol{u}}\in V_U \\ \int_{\Omega}\rho T^\text{ref} \dfrac{\partial s}{\partial t}\widehat{T}d\Omega - \int_{\Omega} \mathbf{j}\cdot\nabla \widehat{T}d\Omega &= -\int_{\partial \Omega} \mathbf{j}\cdot\boldsymbol{n} \widehat{T} dS \quad \forall \widehat{T} \in V_T \end{aligned} \qquad{(3)}\]

with \(V_U\) and \(V_T\) being the displacement and temperature function spaces.

The time derivative in the heat equation is now replaced by an implicit Euler scheme, so that the previous weak form at the time increment \(n+1\) is now:

\[ \int_{\Omega}\rho T^\text{ref} \dfrac{s^{n+1}-s^n}{\Delta t}\widehat{T}d\Omega - \int_{\Omega} \mathbf{j}^{n+1}\cdot\nabla \widehat{T}d\Omega = -\int_{\partial \Omega} \mathbf{j}^{n+1}\cdot\boldsymbol{n} \widehat{T} dS \quad \forall \widehat{T} \in V_T \]

where \(s^{n+1},\mathbf{j}^{n+1}\)
correspond to the *unknown* entropy and heat flux at time \(t_{n+1}\).

Since both the entropy and the stress tensor depend on the temperature and the total strain, we obtain a fully coupled problem at \(t=t_{n+1}\) for \((\boldsymbol{u}_{n+1},T_{n+1})\in V_U\times V_T\). With the retained boundary conditions both right-hand sides in (3).

We now load the material behaviour and define the corresponding
`MFrontNonlinearProblem`

. One notable specificity of the
present example is that the unknown field `v`

belongs to a
mixed function space. Therefore, we cannot rely on automatic
registration for the strain and temperature gradient. We will have to
specify explicitly their UFL expression with respect to the displacement
`u`

and temperature variation `Theta`

sub-functions of the mixed unknown `v`

. We also register the
`"Temperature"`

external state variable with respect to
`Theta`

.

```
= mf.MFrontNonlinearMaterial("./src/libBehaviour.so",
material "ThermoElasticity",
="plane_strain")
hypothesis= Constant(material.get_parameter("MassDensity"))
rho = Function(V)
v = split(v)
(u, Theta) = mf.MFrontNonlinearProblem(v, material, quadrature_degree=2, bcs=bcs)
problem "Strain", sym(grad(u)))
problem.register_gradient("TemperatureGradient", grad(Theta))
problem.register_gradient("Temperature", Theta + Tref) problem.register_external_state_variable(
```

Similarly to the Transient
heat equation with phase change demo, we need to specify explicitly
the coupled thermo-mechanical residual expression using the stress, heat
flux and entropy variables. For the implicit Euler scheme, we will need
to define the entropy at the previous time step. For the mechanical
residual, note that the stress variable `sig`

is represented
in the form of its vector of components. The computation of \(\boldsymbol{\sigma}:\nabla^s\widehat{\boldsymbol{u}}\)
therefore requires to express \(\widehat{\boldsymbol{\varepsilon}}=\nabla^s\widehat{\boldsymbol{u}}\)
in the same way. For this purpose, we could use the
`mgis.fenics.utils.symmetric_tensor_to_vector`

on the
tensorial UFL expression `sym(grad(u))`

. Another possibility
is to get the corresponding `"Strain"`

gradient object
(expressed in vectorial form) and get his variation with respect to
`v_`

.

```
= problem.get_flux("Stress")
sig = problem.get_flux("HeatFlux")
j = problem.get_state_variable("EntropyPerUnitOfMass")
s
problem.initialize()
= s.copy(deepcopy=True)
s_old = TestFunction(V)
v_ = split(v_) # Displacement and temperature test functions
u_, T_ = problem.gradients["Strain"].variation(v_)
eps_ = dot(sig, eps_)*problem.dx
mech_residual = (rho*Tref*(s - s_old)/dt*T_ - dot(j, grad(T_)))*problem.dx
thermal_residual = mech_residual + thermal_residual
problem.residual problem.compute_tangent_form()
```

The problem is now solved by looping over time increments. Because of
the typical exponential time variation of temperature evolution of the
heat equation, time steps are discretized on a non-uniform (logarithmic)
scale. \(\Delta t\) is therefore
updated at each time step. The previous entropy field `s_old`

is updated at the end of each step.

```
= 10
Nincr = np.logspace(1, 4, Nincr + 1)
t = 100
Nx = np.linspace(R, L, Nx)
x = np.zeros((Nx, Nincr + 1))
T_res for (i, dti) in enumerate(np.diff(t)):
print("Increment " + str(i + 1))
dt.assign(dti)
problem.solve(v.vector())
s_old.assign(s)
+ 1] = [v(xi, 0.0)[2] for xi in x] T_res[:, i
```

```
Increment 1
Increment 2
Increment 3
Increment 4
Increment 5
Increment 6
Increment 7
Increment 8
Increment 9
Increment 10
```

At each time increment, the variation of the temperature increase
\(\Delta T\) along a line \((x, y=0)\) is saved in the
`T_res`

array. This evolution is plotted below. As expected,
the temperature gradually increases over time, reaching eventually a
uniform value of \(+10^{\circ}\text{C}\) over infinitely long
waiting time. We check that we obtain the same solution as the pure
`FEniCS`

demo.

```
%matplotlib notebook
plt.figure()1::Nincr // 10])
plt.plot(x, T_res[:, "$x$-coordinate along $y=0$")
plt.xlabel("Temperature variation $\Delta T$")
plt.ylabel("$t={:.0f}$".format(ti) for ti in t[1::Nincr // 10]], ncol=2)
plt.legend([ plt.show()
```

` <IPython.core.display.Javascript object>`