Source files:

- Jupyter notebook: mgis_fenics_nonlinear_heat_transfer.ipynb
- Python file: mgis_fenics_nonlinear_heat_transfer.py
- MFront behaviour file: StationaryHeatTransfer.mfront

The thermal material is described by the following non linear Fourier Law:

\[ \mathbf{j}=-k\left(T\right)\,\mathbf{\nabla} T \]

where \(\mathbf{j}\) is the heat flux and \(\mathbf{\nabla} T\) is the temperature gradient.

The thermal conductivity is assumed to be given by:

\[ k\left(T\right)={\displaystyle \frac{\displaystyle 1}{\displaystyle A+B\,T}} \]

This expression accounts for the phononic contribution to the thermal conductivity.

As discussed below, the consistent linearisation of the heat transfer equilibrium requires to compute:

- the derivative \({\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial \mathbf{\nabla} T}}\) of the heat flux with respect to the temperature gradient. \({\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial \mathbf{\nabla} T}}\) is given by: \[ {\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial \mathbf{\nabla} T}}=-k\left(T\right)\,\matrix{I} \]
- the derivative \({\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial T}}\) of the heat flux with respect to the temperature. \({\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial T}}\) is given by: \[ {\displaystyle \frac{\displaystyle \partial \mathbf{j}}{\displaystyle \partial T}}=-{\displaystyle \frac{\displaystyle \partial k\left(T\right)}{\displaystyle \partial T}}\,\mathbf{\nabla} T=B\,k^{2}\,\mathbf{\nabla} T \]

`MFront`

’ implementationEvery `MFront`

file is handled by a domain specific
language (DSL), which aims at providing the most suitable abstraction
for a particular choice of behaviour and integration algorithm. See
`mfront mfront --list-dsl`

for a list of the available
DSLs.

The name of DSL’s handling generic behaviours ends with
`GenericBehaviour`

. The first part of a DSL’s name is related
to the integration algorithm used.

In the case of this non linear transfer behaviour, the heat flux is
explicitly computed from the temperature and the temperature gradient.
The `DefaultGenericBehaviour`

is the most suitable
choice:

`@DSL DefaultGenericBehaviour;`

The following lines define the name of the behaviour, the name of the author and the date of its writing:

```
@Behaviour StationaryHeatTransfer;
@Author Thomas Helfer;
@Date 15/02/2019;
```

Generic behaviours relate pairs of gradients and fluxes. Gradients and fluxes are declared independently but the first declared gradient is assumed to be conjugated with the first declared fluxes and so on…

The temperature gradient is declared as follows (note that Unicode characters are supported):

```
@Gradient TemperatureGradient ∇T;
.setGlossaryName("TemperatureGradient"); ∇T
```

Note that we associated to `∇T`

the glossary name
`TemperatureGradient`

. This is helpful for the calling
code.

After this declaration, the following variables will be defined:

- The temperature gradient
`∇T`

at the beginning of the time step. - The increment of the temperature gradient
`Δ∇T`

over the time step.

The heat flux is then declared as follows:

```
@Flux HeatFlux j;
.setGlossaryName("HeatFlux"); j
```

In the following code blocks, `j`

will be the heat flux at
the end of the time step.

By default, the derivatives of the gradients with respect to the
fluxes are declared. Thus the variable `∂j∕∂Δ∇T`

is
automatically declared.

However, as discussed in the next section, the consistent linearisation of the thermal equilibrium requires to return the derivate of the heat flux with respect to the increment of the temperature (or equivalently with respect to the temperature at the end of the time step).

`@AdditionalTangentOperatorBlock ∂j∕∂ΔT;`

The `A`

and `B`

coefficients that appears in
the definition of the thermal conductivity are declared as
parameters:

```
@Parameter real A = 0.0375;
@Parameter real B = 2.165e-4;
```

Parameters are stored globally and can be modified from the calling
solver or from `python`

in the case of the coupling with
`FEniCS`

discussed below.

A local variable is accessible in each code blocks.

Here, we declare the thermal conductivity `k`

as a local
variable in order to be able to compute its value during the behaviour
integration and to reuse this value when computing the tangent
operator.

`@LocalVariable thermalconductivity k;`

The behaviour integration is straightforward: one starts to compute the temperature at the end of the time step, then we compute the thermal conductivity (at the end of the time step) and the heat flux using the temperature gradient (at the end of the time step).

```
@Integrator{
// temperature at the end of the time step
const auto T_ = T + ΔT;
// thermal conductivity
= 1 / (A + B ⋅ T_);
k // heat flux
= -k ⋅ (∇T + Δ∇T);
j } // end of @Integrator
```

The computation of the tangent operator blocks is equally simple:

```
@TangentOperator {
= -k ⋅ tmatrix<N, N, real>::Id();
∂j∕∂Δ∇T = B ⋅ k ⋅ k ⋅ (∇T + Δ∇T);
∂j∕∂ΔT } // end of @TangentOperator
```

`FEniCS`

implementationWe consider a rectanglar domain with imposed temperatures
`Tl`

(resp. `Tr`

) on the left (resp. right)
boundaries. We want to solve for the temperature field `T`

inside the domain using a \(P^1\)-interpolation. We initialize the
temperature at value `Tl`

throughout the domain.

```
%matplotlib notebook
import matplotlib.pyplot as plt
from dolfin import *
import mgis.fenics as mf
= 30e-3
length = 5.4e-3
width = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
mesh
= FunctionSpace(mesh, "CG", 1)
V = Function(V, name="Temperature")
T
def left(x, on_boundary):
return near(x[0], 0) and on_boundary
def right(x, on_boundary):
return near(x[0], length) and on_boundary
= 300
Tl = 800
Tr
T.interpolate(Constant(Tl))
= [DirichletBC(V, Constant(Tl), left),
bc DirichletBC(V, Constant(Tr), right)]
```

We use the `MFrontNonlinearMaterial`

class for describing
the material behaviour. The first argument corresponds to the path where
material librairies have been compiled, the second correspond to the
name of the behaviour (declared with `@Behaviour`

). Finally,
the modelling hypothesis is specified (default behaviour is
`"3d"`

).

```
= mf.MFrontNonlinearMaterial("./src/libBehaviour.so",
material "StationaryHeatTransfer",
="plane_strain") hypothesis
```

The `MFront`

behaviour declares the field
`"TemperatureGradient"`

as a Gradient variable, with its
associated Flux called `"HeatFlux"`

. We can check that the
`material`

object retrieves `MFront`

’s gradient
and flux names, as well as the different tangent operator blocks which
have been defined, namely `dj_ddgT`

and `dj_ddT`

in the present case:

```
print(material.get_gradient_names())
print(material.get_flux_names())
print(["d{}_d{}".format(*t) for t in material.get_tangent_block_names()])
```

```
['TemperatureGradient']
['HeatFlux']
['dHeatFlux_dTemperatureGradient', 'dHeatFlux_dTemperature']
```

When defining the non-linear problem, we will specify the boundary conditions and the requested quadrature degree which will control the number of quadrature points used in each cell to compute the non-linear constitutive law. Here, we specify a quadrature of degree 2 (i.e. 3 Gauss points for a triangular element).

`= mf.MFrontNonlinearProblem(T, material, quadrature_degree=2, bcs=bc) problem `

The `MFront`

behaviour implicitly declares the temperature
as an external state variable called `"Temperature"`

. We must
therefore associate this external state variable to a known mechanical
field. This can be achieved explicitly using the
`register_external_state_variable`

method. In the present
case, this can be done automatically since the name of the unknown
temperature field matches the TFEL Glossary name
`"Temperature"`

which is a predefined external state
variable. In this case, the following message will be printed:

`Automatic registration of 'Temperature' as an external state variable.`

For problems in which the temperature only acts as a parameter (no
jacobian blocks with respect to the temperature), the temperature can be
automatically registered as a constant value (\(293.15 \text{ K}\) by default) or to any
other (`dolfin.Constant`

, `float`

or
`dolfin.Function`

) value using the
`register_external_state_variable`

method.

In the `FEniCS`

interface, we instantiate the main
mechanical unknown, here the temperature field `T`

which has
to be named `"Temperature"`

in order to match
`MFront`

’s predefined name. Using another name than this will
later result in an error saying:

`ValueError: 'Temperature' could not be associated with a registered gradient or a known state variable.`

Finally, we need to associate to `MFront`

gradient object
the corresponding UFL expression as a function of the unknown field
`T`

. To do so, we use the `register_gradient`

method linking `MFront`

`"TemperatureGradient"`

object to the UFL expression `grad(T)`

.

`"TemperatureGradient", grad(T)) problem.register_gradient(`

Similarly to the case of external state variables, common gradient
expressions for some TFEL Glossary names
have been already predefined which avoid calling explicitly the
`register_gradient`

method. Predefined expressions can be
obtained from:

` mf.list_predefined_gradients()`

```
'TFEL gradient name' (Available hypotheses)
---------------------------------------------
'Strain' ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')
'TemperatureGradient' ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')
'DeformationGradient' ('Tridimensional', 'PlaneStrain', 'Axisymmetrical')
```

We can see that the name `"Temperature Gradient"`

is in
fact a predefined gradient. Omitting calling the
`register_gradient`

method will in this case print the
following message upon calling `solve`

:

`Automatic registration of 'TemperatureGradient' as grad(Temperature).`

meaning that a predefined gradient name has been found and registered as the UFL expression \(\nabla T\). > Note that automatic registration is not supported when using mixed function spaces.

No external loading has been specified so that the residual form is automatically defined as: \[\begin{equation} F(\widehat{T}) = \int_\Omega \mathbf{j}\cdot \nabla \widehat{T} \text{dx} = 0 \quad \forall \widehat{T} \end{equation}\]

From the two tangent operator blocks `dj_ddgT`

and
`dj_ddT`

, it will automatically be deduced that the heat flux
\(\mathbf{j}\) is a function of both
the temperature gradient \(\mathbf{g}=\nabla
T\) and the temperature itself i.e. \(\mathbf{j}=\mathbf{j}(\mathbf{g}, T)\). The
following tangent bilinear form will therefore be used when solving the
above non-linear problem:

\[\begin{equation} a_\text{tangent}(\widehat{T},T^*) = \int_{\Omega} \nabla \widehat{T}\cdot\left(\dfrac{\partial \mathbf{j}}{\partial \mathbf{g}}\cdot \nabla T^*+\dfrac{\partial \mathbf{j}}{\partial T}\cdot T^*\right) \text{dx} \end{equation}\]

We finally solve the non-linear problem using a default Newton
non-linear solver. The `solve`

method returns the number of
Newton iterations (4 in the present case) and converged status .

` problem.solve(T.vector())`

```
Automatic registration of 'Temperature' as an external state variable.
(4, True)
```

We finally check that the thermal conductivity coefficient \(k\), computed from the ratio between the
horizontal heat flux and temperature gradient matches the
temperature-dependent expressions implemented in the `MFront`

behaviour.

```
= problem.fluxes["HeatFlux"].function
j = problem.gradients["TemperatureGradient"].function
g = -j.vector().get_local()[::2]/g.vector().get_local()[::2]
k_gauss = problem.state_variables["external"]["Temperature"].function.vector().get_local()
T_gauss = material.get_parameter("A");
A = material.get_parameter("B");
B = 1/(A + B*T_gauss)
k_ref 'o', label="FE")
plt.plot(T_gauss, k_gauss, '.', label="ref")
plt.plot(T_gauss, k_ref, r"Temperature $T\: (K)$")
plt.xlabel(r"Thermal conductivity $k\: (W.m^{-1}.K^{-1})$")
plt.ylabel(
plt.legend() plt.show()
```

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