2. Principle of the maximum#

2.1. Statement of principle for the continuous case#

One gives here one of the possible statements of the principle of maximum for the heat operator (in the absence of source terms, and in isotropic homogeneous linear thermics) [bib2].

Let \(\Omega\) be an open bounded by \({ℝ}^{n}\) with a border of \(\Gamma\), whose grip is rated \(\stackrel{ˉ}{\Omega }\).

Let \(u(\text{x,t})\) be such that:

\(\frac{\mathrm{\partial }u}{\mathrm{\partial }t}\mathrm{-}\Delta u\mathrm{\le }0\text{sur}\Omega \mathrm{\times }\mathrm{]}\mathrm{0,}T\mathrm{[},(T>0)\)

Of class \({C}^{2}\) compared to \(x\) and \(u\) of class \({C}^{1}\) compared to \(t\) out of \(\Omega \text{'}\mathrm{\times }\mathrm{]}\mathrm{0,}T\mathrm{[}\)

So \(\underset{\overline{\Omega }\mathrm{\times }\left[\mathrm{0,}T\right]}{\text{Max}}u\mathrm{=}\underset{P}{\text{Max}}u\), where \(P\mathrm{=}(\overline{\Omega }\mathrm{\times }\left\{0\right\})\mathrm{\cup }(\Gamma \mathrm{\times }\left[\mathrm{0,}T\right])\) is the border of cylinder \(\Omega \mathrm{\times }\mathrm{]}0\text{,T}\mathrm{[}\).

This result therefore ensures that the maximum of \(u\) is necessarily reached either during the initial conditions or on one edge of the domain during the transition.

2.2. Respect for the principle of maximum discretion#

We consider the heat equation (thermal conduction):

\(\text{div}(\lambda \text{.}\nabla T)+s(x,t)={\mathrm{\rho C}}_{p}\frac{\partial T}{\partial t}\) + Boundary conditions + initial condition \(T({t}_{0}\text{,x})={T}^{0}(x)\)

with

\(T\)

temperature

\(s\)

heat per unit volume (internal sources)

\(t\)

time variable

\(x\)

space variable

\(\lambda\)

thermal conductivity coefficient

\({\mathrm{\rho C}}_{p}\)

volume heat at constant pressure

Boundary condition types (linear) :

  • Imposed temperature: Dirichlet limit condition

\(T(x,t)={T}_{\text{imp}}(x,t)\text{sur}{\Gamma }_{\text{imp}}\)

  • Imposed normal flow: Neumann condition defining the flow entering the domain

\(-q(x,t)\text{.}n=f(x,t)\text{sur}{\Gamma }_{\text{flux}}\)

  • Exchange: Fourier limit condition modeling convective exchanges at the edges of the domain

\(-q(x,t)\text{.}n=h(x,t)({T}_{\text{ext}}(x,t)-T(x,t))\text{sur}{\Gamma }_{\text{échange}}\)

The variational formulation of the problem is as follows: [bib3]

\(\underset{\Omega }{\int }{\mathrm{\rho C}}_{p}\frac{\partial T}{\partial t}\text{.}\mathrm{\nu d}\Omega +\underset{\Omega }{\int }\lambda \nabla T\text{.}\nabla \nu d\Omega +\underset{{\Gamma }_{\text{Žéchange}}}{\int }\mathrm{hT}\text{.}\nu d\Gamma =\underset{\Omega }{\int }s\text{.}\nu d\Omega +\underset{{\Gamma }_{\text{flux}}}{\int }f\text{.}\nu d\Gamma +\underset{{\Gamma }_{\text{Žéchange}}}{\int }{\mathrm{hT}}_{\text{ext}}\text{.}\nu d\Gamma\)

\(\forall v\) checking \(v={T}_{\text{imp}}(x,t)\text{sur}{\Gamma }_{\text{imp}}\)

After discretizing this equation in space, we obtain the system:

\(M\left\{\frac{\partial T}{\partial t}(t)\right\}+\text{KT}(t)=F(t)\).

with

\(T(t)\): vector of nodal temperatures

\(M\): thermal mass matrix

\(M=\sum _{e}\underset{{\Omega }_{e}}{\int }\mathrm{\rho C}N{N}^{T}\text{dV}\)

\(K\): thermal stiffness matrix

\(K=\sum _{e}(\underset{{\Omega }_{e}}{\int }\lambda \nabla N\text{.}\nabla {N}^{T}\text{dV}+\underset{{\Gamma }_{{\text{échange}}_{e}}}{\int }hN\text{.}{N}^{T}d\Gamma )\)

\(F\): vector of the second member

\(F=\sum _{e}(\underset{{\Omega }_{e}}{\int }sNd\Omega +\underset{{\Gamma }_{{\text{flux}}_{e}}}{\int }fNd\Gamma +\underset{{\Gamma }_{{\text{échange}}_{e}}}{\int }{\mathrm{hT}}_{\text{ext}}Nd\Gamma )\)

\(N\): (shape functions)

For time discretization, we apply a \(\theta\) -\((\theta \in [\mathrm{0,1}])\) method, which leads to:

\((M+\theta \Delta tK){T}^{n+1}=(M+(\theta -1)\Delta tK){T}^{n}+\theta \langle {F}^{n+1}\rangle +(1-\theta )\langle {F}^{n}\rangle\)

where \({T}^{n},{T}^{n+1}\) are the vectors of the nodal temperatures at the times \({t}_{n}{\text{,t}}_{\text{n+}1}\).

2.3. Sufficient conditions for compliance with the principle of maximum at the discrete level#

One of the characteristics of non-compliance with the principle of maximum is the appearance of oscillations (temporal or spatial): if we observe the variation in temperature at a node over time, we note that the solution oscillates and exceeds the minimum and maximum values determined by the minimum and maximum values determined by the initial conditions and the limit conditions. Or again, at a given moment, spatial oscillations are observed.

We are therefore looking for sufficient conditions on \(\Delta t\), \(K\) and \(M\) so that the solution does not oscillate over time ([bib1], [bib4], [bib5]). In fact, we do not know how to obtain necessary and sufficient conditions. We therefore seek conditions under which the solution does not oscillate over time. If these are verified, it will be verified that the spatial oscillations have also disappeared, and then the respect of the maximum principle is ensured.

Assumptions:

To be able to express these sufficient conditions of non-oscillation, two hypotheses must be added:

  • we are going to the elementary level. Compliance with the properties at the elementary level is sufficient for the non-oscillation conditions to be verified for the assembled matrices.

  • it is considered that the stiffness matrix \(K\) is only formed from the volume term \({K}_{V}=\underset{{\Omega }_{e}}{\int }\lambda \nabla N\text{.}\nabla {N}^{T}\text{dV}\)

This assumption is not valid for all boundary conditions (see paragraph :ref:`3.1.3 <RefHeading__34144941>`) . *

Sufficient non-oscillation conditions amount to expressing conditions on the time step and on the diagonal and extra-diagonal terms of \(M\) and \(K\) for certain properties of these matrices to be verified (based on the monotony of the matrices) [bib1]:

\({M}_{\text{ij}}+\theta \text{.}\Delta {\text{tK}}_{\text{ij}}\le 0\text{}i\ne j\) eq 2.3-1

\({M}_{\text{ij}}+(\theta -1)\Delta {\text{tK}}_{\text{ij}}\ge 0\text{}i\ne j\) eq 2.3-2

\({M}_{\text{ii}}+(\theta -1)\Delta {\text{tK}}_{\text{ii}}\ge 0\text{}\forall i\) eq 2.3-3

In the general case, extra-diagonal terms can be of any sign. A quick study makes it possible to determine the conditions on t according to their signs so that the previous equations are verified:

\({K}_{\text{ij}}\ge 0\)

\({K}_{\text{ij}}\le 0\)

\({M}_{\text{ij}}\ge 0\)

\({M}_{\text{ij}}+\theta \text{.}\Delta {\text{tK}}_{\text{ij}}\le 0\text{}i\ne j\)

[éq 2.3-1] unconditionally false except \({M}_{\text{ij}}={K}_{\text{ij}}=0\)

\(\underset{i\ne j}{\text{max}}(\frac{{M}_{\mathrm{ij}}}{-\theta {K}_{\text{ij}}})\le \Delta t\le \underset{i}{\text{min}}(\frac{{M}_{\text{ii}}}{(1-\theta ){K}_{\text{ii}}})\)

\({M}_{\text{ij}}\le 0\)

\({M}_{\text{ij}}+(\theta -1)\Delta {\text{tK}}_{\text{ij}}\ge 0\text{}i\ne j\)

[éq 2.3-2] unconditionally false except \({M}_{\text{ij}}={K}_{\text{ij}}=0\)

\(\underset{i\ne j}{\text{max}}(\frac{{M}_{\text{ij}}}{(1-\theta ){K}_{\text{ij}}})\le \Delta t\le \underset{i}{\text{min}}(\frac{{M}_{\text{ii}}}{(1-\theta ){K}_{\text{ii}}})\)

Regardless of \(\Delta t\) and the shape of \(M\), there is risk of oscillations.

Interval to be respected on \(\Delta t\). The diagonalization of \(M\) allows to remove the lower bound.

The conditions sufficient to avoid oscillations are then:

\({K}_{\text{ij}}\le 0\text{}i\ne j\)

\(\Delta {t}_{\text{min}}\le \Delta t\le \Delta {t}_{\text{max}}\)

with:

  • \(\Delta {t}_{\text{max}}=\underset{i}{\text{min}}(\frac{{M}_{\text{ii}}}{(1-\theta ){K}_{\text{ii}}})\) and

  • \(\underset{}{\Delta {t}_{\text{min}}=\underset{i\ne j}{\text{max}}}(\frac{{M}_{\text{ij}}}{(1-\theta ){K}_{\text{ij}}},\frac{{M}_{\text{ij}}}{-\theta {K}_{\text{ij}}})\),

Consequently, the elementary matrices must first verify \({K}_{\text{ij}}\le 0\) (this is the case of the linear finite elements studied later).

As for the interval on the time step:

If the oscillations are due to a time step that is too big \((\Delta t>\Delta {t}_{\text{max}})\), we can advise:

  • or to choose an implicit time integration scheme \((\theta =1)\), to eliminate the upper bound of the interval.

  • or to reduce \(\Delta t\). (In practice, it is difficult to know an order of magnitude of \(\Delta {t}_{\text{max}}\)).

Quite often, the problem of oscillations arises for small time steps \((\Delta t<\Delta {t}_{\text{min}})\); in fact, to take into account the variations of the solution (for example, during a thermal shock), we are led to choose a discretization that is fine in time. In this case, to avoid oscillations, we can suggest:

  • or to increase the value of the time steps. In practice, this is not always possible because \(\mathrm{Dt}\) can be imposed by the nature of the problem (rapid load variation). Plus, it’s hard to get an order of magnitude of \(\Delta {t}_{\text{min}}\).

  • or to reduce the size of the meshes and therefore increase the number of elements. In fact, the value of \(\Delta {t}_{\text{min}}\) depends directly on spatial discretization:

The expressions of the elementary matrices are in fact:

\({M}_{\text{ij}}=\underset{{\Omega }_{e}}{\int }\rho {\text{CN}}_{i}{N}_{j}\text{dV}\)

\({K}_{\text{ij}}=\underset{{\Omega }_{e}}{\int }\lambda \nabla {N}_{i}\nabla {N}_{j}\text{dV}\)

For 2D elements, the terms in \(M\) are therefore in the form \(\mathrm{\rho C}\times \mathrm{Surface}\) while those in \(K\) are only a function of \(\lambda\). This solution remains the best if one is not limited by the calculation cost, because the thermal and especially mechanical solution will be all the more accurate.

  • Or to diagonalize the matrix \(M\), which removes the lower bound of the interval.

That is the solution proposed here.

In the rest of the study, we are only interested in the problem of oscillations that occur for too small time steps: \(\Delta t<\Delta {t}_{\text{min}}\). More precisely, the method of diagonalization of the matrix \(M\) chosen is presented, and the different types of elements to which it applies.