# Numerical simulation of heat extraction by a coaxial ground heat exchanger
under freezing conditions

G.P. Vasilyev  N.V. Peskov<EMAIL_ADDRESS>T.M. Lysak  OAO “INSOLAR-INVEST”,
Moscow, Russian Federation.  Faculty of Computational Mathematics and
Cybernetics,
Lomonosov Moscow State University, Moscow, Russian Federation.

###### Abstract

A new finite-difference model of heat transfer inside a shallow coaxial ground
heat exchanger and in the surrounding layered soil is presented, taking into
account the freezing of ground moisture. Three modes of heat exchanger
operation are numerically simulated: stationary mode, transient mode and
controlled mode. In the stationary mode, estimates of the sensitivity of the
heat carrier fluid outlet temperature to changes in the heat exchanger
parameters are calculated. In all modes, close attention is paid to
demonstrating the difference in the results at a negative temperature of the
fluid, calculated with and without taking into account the freezing of ground
moisture. It is shown that this difference, caused by the zero-curtain effect,
can range from 10% in the stationary mode to 35% in the control mode.

###### keywords:

coaxial heat exchanger, numerical simulations, freezing of ground moisture,
zero-curtain effect

##  1 Introduction

Ground source heat pump (GSHP) is a highly efficient renewable energy
technology which allows sufficient reduction in energy consumption for
building heating and cooling [1]. An important part of each GSHP system is a
ground heat exchanger (GHE) installed in a vertical borehole or horizontal
trench. Thermal energy is extracted from the ground or absorbed by the ground
due to the effective heat exchange between the fluid, which circulates in the
GHE, and the adjacent soil. The effective design of GHEs requires a good
understanding of the thermal processes in the system GHE-surrounding soil.

Three main approaches to the investigation of the system GHE-surrounding soil
can be distinguished: analytical, numerical and experimental. Analytical
approach allows obtaining analytical dependencies under a number of strong
simplifications about the thermal processes. More complicated 2D or 3D
geometry models require numerical simulation, but use the results of
analytical approach and experimental approaches for their testing. A lot of
information on the modeling of thermal processes for various types of GHE can
be found in reviews [2, 3]. In [3], a list of commonly used software packages
for GSHP design can be found as well.

While most models focus on normal GHE operation with above-zero inlet
temperature, the effect of phase change on GHE operation is of great
importance in cold regions, where the circulating fluid inlet temperature must
be below zero to extract the required amount of heat. Freezing-thawing
processes around the GHE change the thermal characteristics of the soil thus
affecting the thermal exchange between the pipe fluid and the surrounding
soil. In particular, a number of experiments show that freezing increases the
heat transfer rate and the coefficient of performance [4, 5, 6, 7, 8, 9, 10].
As a result, utilization of latent heat from the groundwater freezing allows,
for example, the reduction of vertical GHE length for the GSHP with the same
characteristics [5].

Another consequence is the possible damage of the pipe due to the ice lens
formation in the grout [11] and the impact of freeze-thaw cycles on GHE
hydraulic conductivity [12, 13, 14, 15]. In [13], the influence of freeze-thaw
cycles on aging effect of U-shaped GHEs were studied using the large-scale
test rig determining the hydraulic conductivity of GHE systems. In [12], a
pilot-scale experiment was built to test a 1-m section of a typical GHE under
freezing-thawing conditions. Freeze-thaw cycles seem to be a leading cause for
a potential violation of the system’s integrity [12, 14]. In [15], the thermal
processes in the saturated porous medium around the GHE were simulated using
FEFLOW 2D axisymmetric model with C++ plug-in for phase changes between solid
and liquid phases. Possible ways to avoid GHE pipes deformation and reduce
flow resistance under freezing conditions were studied via experimental tests
in [16] and using numerical simulation in [17].

Note that the damage of the pipes is mainly studied by experimental tests. On
the contrary, the role of numerical simulation in studying the efficiency of
thermal transfer processes in the system GHE-surrounding soil is very large.
Up to now a number of math models for this aim was developed, which take into
account freezing-thawing processes in the soil around the vertical and
horizontal GHEs. The major part of all math models is the heat transfer
equation for the surrounding soil. The freezing-thawing processes are
described in the model using a latent heat term in the heat transfer equation
and the heat transfer model [6] with effective heat capacity to describe the
three phases – solid, liquid and mushy. This method allows avoiding one of the
difficulties in numerical simulation for phase change: how to treat the
movable interface between the phases. The other important part of math model
is temperature distribution inside the GHE pipe. Some investigations do not
consider this temperature distribution, operating with the total heat
gain/loss in GHE [6, 18]. The others (see [19], for example ) use a dual
continuum approach [20]. Under such approach, the soil is described in 3D
geometry, while the two branches of U-shaped GHE are described in 1D geometry.
To simplify the description of the temperature distribution inside the
U-shaped pipe, the pipe is replaced by a mono-tube with equivalent diameter
[21].

Without a doubt, detailed 3D or even 2D calculation of thermal processes in
the GHE and surrounding soil can provide accurate results. But such
calculations are in most cases based on commercial software packages, such as
FEFLOW, COMSOL, ABAQUS, OpenGeoSys , CFD ANSYS FLUENT and so on. Commercial
software packages are also very useful at modeling the arrays of GHEs [7, 22,
23]. An obvious drawback of these packages is the impossibility of numerical
modeling of phase transfer processes in porous soil. To overcome this
drawback, some authors use C++ or FORTRAN plug-ins [15, 17] to commercial
software package for modeling phase transition between solid and liquid phases
in porous surrounding soil. Another disadvantage of detailed 3D simulation is
the large time of calculation and large computational resources necessary for
the long-term simulation of the system GHE-surrounding soil with phase change
processes.

To reduce the computational resources, simplified 1D and 2D model for vertical
GHEs were developed. For example in [24], an improved thermal resistance and
capacity model (RC) was proposed for the heat transfer modeling between
vertical single U-tube GHEs and the frozen soil. The model was verified
through experiment tests and numerical simulation on basis of 3D CFD ANSYS
model. In [6], a 2D model for temperature field modeling around the GHE was
developed. The model is based on heat transfer equation in cylindrical
geometry with a source term and Dirichlet boundary conditions to represent the
GHE. In [5], 1D radial numerical heat transfer model is developed to evaluate
heat transfer from the borehole wall to the ground with taking into account
three phases: ice-soil, water-soil and transition phase. For the boundary
conditions at the borehole wall and the outer boundary of computational area
the temperature, measured experimentally, was used. Another approach to reduce
the computational time was demonstrated in [18] for the modeling of a GHE
array applied to an existing historical building in Venice (Northern Italy).
In this paper, two different scale models were developed. The results of the
coarser large scale model were used as initial and boundary conditions for the
fully discretized small scale model, developed for the detailed description of
freezing-thawing processes in the close vicinity of the GHE.

The processes of freezing-thawing are closely related to the processes of
groundwater seepage and flow. The influence of these two factors – freezing-
thawing together with ground water seepage, are taken into account in [7, 10,
22] in modeling the temperature field in the surrounding soil for the array of
vertical U-shaped GHEs. The temperature field equation for the porous
surrounding soil with convection and ice/water phase and Brinkman equation for
water transport in porous media are the basis of the 3D model, which is solved
by the COMSOL software. The model does not describe the temperature inside the
pipes. Instead it uses the iteration processes to determine the inlet and
output temperature on basis of the soil model. The similar model was developed
in [23] for artificial freezing.

This paper proposes a combined mathematical model containing a 1D model of
heat transfer in a fluid circulating inside a vertical coaxial GHE and a 2D
model of heat conduction in a cylindrically symmetric surrounding soil. The
model takes into account temperature changes in the inner and annular pipes,
soil freezing processes, geothermal temperature gradient and horizontal soil
stratification. Numerical implementation of the model, based on the finite
difference method and the MatLab ode15s numerical solver, makes it possible to
investigate the effect of soil freezing on the outlet temperature of the
coaxial heat exchanger.

Coaxial GHE have a number of advantages in comparison with U-shaped one,
especially in the case of deep boreholes (1000 to 3000 meters depth). Among
them is a simple installation procedure, moderate temperature difference
between the secondary fluid and the surrounding ground, the possibility to use
water as a secondary fluid even in colder countries and to ignore its local
thermal resistance in comparison with U-shaped GHE [25, 26]. The comparison of
different types of GHE can be found in [27], where the equivalent ground
thermal conductivity was evaluated by Thermal Response Test. A review on
various approaches for investigation of coaxial GHEs can be found, for
example, in [28].

Geothermal gradient was taken into account in an analytical model [29] and
semi-analytical models [28, 30] proposed for shallow (100 to 200 m depth,
[28], [29]) and deep (1000 to 3000 meters depth, [30]) coaxial BHE.

The vertical geological structure of the adjacent ground was considered in
[31] at numerical modeling of a deep coaxial GHE.

However, to our knowledge, freezing-thawing processes have not yet been
considered for coaxial GHE. It should be noted that, in deep GHEs, only pure
water (without anti-freeze) with positive inlet temperature is allowed to use
to prevent pollution on deep groundwater resources [32]. Therefore, in our
paper, we focus on a shallow vertical coaxial GHE with negative inlet
temperature of circulating fluid.

##  2 Formulation of the problem

###  2.1 Scheme of coaxial ground heat exchanger

We imagine a heat exchanger in the form of two coaxial pipes located in a
vertical borehole and fixed in it with the help of a grout. The cross section
of the heat exchanger is shown schematically, not to scale, in the Figure 1.

Figure 1: Out of scale scheme of coaxial ground heat exchanger.

The cross section of the inner pipe is shown in the center of the figure,
$r_{11}$ and $r_{12}$ are the inner and the outer radius of this pipe,
respectively. Coaxial with the inner pipe is the outer pipe, its inner and
outer radius are designated $r_{21}$ and $r_{22}$, respectively. The radius of
the borehole is denoted by $r_{b}$, the space between the outer pipe and the
soil outside the borehole is filled with grout. In the heat exchanger scheme,
we take into account the thickness of the pipe walls and the thickness of the
grout layer to estimate the magnitude of heat flows between the pipes and
between the fluid and the soil.

In the process of heat exchange, the heat-carrying fluid moves along the inner
tube and in the annular space between the tubes in opposite directions. We
will consider a heat exchanger in which the fluid moves from top to bottom
along the annular area between the pipes and rises up the inner pipe. Such a
mode is usually used to extract heat from the ground. The temperature of the
fluid in the inner pipe will be denoted by $T_{1}$, in the space between the
pipes – by $T_{2}$, and through $T_{s}$ we will denote the temperature of the
soil. In the general case, the temperatures $T_{1}$, $T_{2}$, and $T_{s}$ are
unknown functions of time and spatial coordinates.

###  2.2 Mathematical statement of the problem

To calculate the temperature of the fluid and the soil surrounding the
borehole, we will solve the initial-boundary value problem in a circular
cylinder of radius $r_{d}$ and height $H$ with a heat exchanger located along
the vertical axis of the cylinder. The problem will be solved in a cylindrical
coordinate system with the origin located in the center of the upper end of
the inner pipe and the axis $z$ directed downward and coinciding with the axis
of the inner pipe. We assume that the temperature of the fluid $T_{1}$ and
$T_{2}$ can depend only on the time $t$ and the coordinate $z$, $T_{1}(t,z)$
and $T_{2}(t,z)$. We do not take into account the thermal conductivity of the
fluid, assuming that heat in the exchanger is distributed only due to the flow
of the fluid and heat transfer through the walls of pipes and borehole.

Regarding the soil, we assume that its physical properties may depend on the
coordinates $r$ and $z$, but do not depend on the azimuth angle,
$T_{s}=T_{s}(t,r,z)$. The side wall of the cylinder at $r=r_{d}$ is considered
to be thermally insulated, and a constant temperature is maintained on the
bases of the cylinder. Moreover, we will take into account the geothermal
gradient, which is essential for long heat exchangers. Thus, we will not take
into account seasonal changes in temperature on the surface of the earth, as
well as changes in the properties of the surface itself. It can be assumed
that the upper part of the heat exchanger is located at a depth of about 10 m
under the earth surface, where the soil temperature is practically constant
throughout the year and is equal to the average annual air temperature on the
surface.

Based on the assumptions made, the problem of calculating the temperature of
the fluid and the soil can be written in the following form.

$\displaystyle c_{f}\rho_{f}a_{1}\frac{\partial T_{1}(t,z)}{\partial
t}=c_{f}v_{f}\frac{\partial T_{1}(t,z)}{\partial z}-q_{1}(T_{1},T_{2}),$ (1)
$\displaystyle c_{f}\rho_{f}a_{2}\frac{\partial T_{2}(t,z)}{\partial
t}=-c_{f}v_{f}\frac{\partial T_{2}(t,z)}{\partial
z}+q_{1}(T_{1},T_{2})-q_{2}(T_{2},T_{sb}),$ (2) $\displaystyle
c_{s}\rho_{s}\frac{\partial T_{s}(t,r,z)}{\partial
t}=\frac{1}{r}\frac{\partial}{\partial r}\left[r\lambda_{s}\frac{\partial
T_{s}(t,r,z)}{\partial r}\right]+\frac{\partial}{\partial
z}\left[\lambda_{s}\frac{\partial T_{s}(t,r,z)}{\partial z}\right].$ (3)

Here are the constant parameters of the fluid: $c_{f}$, [J/(kg$\cdot$K)], is
the heat capacity; $\rho_{f}$, [kg/m3], is the density; and $v_{f}$, [kg/s],
is the mass velocity. The cross section area of the inner pipe and the annular
region between the pipes are denotes as $a_{1}=\pi r_{11}^{2}$ and $a_{2}=\pi
r_{21}^{2}-\pi r_{12}^{2}$. In contrast to the fluid parameters, the soil
parameters $c_{s}$ and $\rho_{s}$, as well as the thermal conductivity of the
soil, $\lambda_{s}$ [W/(m$\cdot$K)], in the general case, can depend on
spatial coordinates and on temperature $T_{s}$. However, in order not to
complicate the formulation of the problem and the numerical solution, we will
consider the smooth dependence of the soil parameters on the coordinates.

The terms $q_{1}$ and $q_{2}$, [W/m], describe heat transfer between fluids
flowing down and upwards, and between fluid and soil, respectively. As usual,
we calculate the values of these terms from the boundary value problem for the
stationary heat equation in the rings $r_{11}<r<r_{12}$ and $r_{21}<r<r_{b}$
with Robin boundary conditions at the pipe walls and the Dirichlet condition
at the borehole wall. From the solution of these problems,

$\displaystyle
q_{1}(T_{1},T_{2})=\frac{2\pi\lambda_{1}(T_{1}-T_{2})}{\ln\left(\frac{r_{12}}{r_{11}}\right)+\frac{\lambda_{1}}{\alpha_{11}r_{11}}+\frac{\lambda_{1}}{\alpha_{12}r_{12}}},$
(4) $\displaystyle
q_{2}(T_{2},T_{sb})=\frac{2\pi\lambda_{2}(T_{2}-T_{sb})}{\ln\left(\frac{r_{22}}{r_{21}}\right)+\frac{\lambda_{2}}{\lambda_{g}}\ln\left(\frac{r_{b}}{r_{22}}\right)+\frac{\lambda_{2}}{\alpha_{21}r_{21}}}.$
(5)

Here $\lambda_{1},\,\lambda_{2},\,\lambda_{g}$ are the thermal conductivity
coefficients of the walls of the inner and outer pipes and the grout,
respectively; $\alpha_{11},\,\alpha_{12},\,\alpha_{21}$ are the coefficients
of convective heat transfer between the fluid and pipe walls: $T_{sb}$ is the
temperature of the soil adjacent to the borehole wall.

Note that the equations, similar to equations (1)-(5), were written in [31]
for a deep coaxial GHE, without taking into account the grout and the
freezing-thawing processes.

The equations (1), (2) are defined in segment $0<z<H$ with the boundary
conditions

$T_{1}(t,H)=T_{2}(t,H);\;T_{2}(t,0)=T_{in},$ (6)

where $T_{in}$ is the fluid temperature at the heat exchanger inlet. In
addition, we will need the outlet temperature $T_{out}(t)=T_{1}(t,0)$, as well
as the heat extraction rate

$q_{ex}(t)=c_{f}v_{f}(T_{out}(t)-T_{in}).$ (7)

The equation (3) is defined in the rectangle $\\{r_{b}<r<r_{d},\,0<z<H\\}$. In
accordance with the assumptions made above, we set the boundary conditions for
Eq. (2) as follows

$\displaystyle\lambda_{s}r_{b}\frac{\partial T_{s}(t,r,z)}{\partial
r}\bigg{|}_{r=r_{b}}=q_{2},\;\frac{\partial T_{s}(t,r,z)}{\partial
r}\bigg{|}_{r=r_{d}}=0,$ (8) $\displaystyle
T_{s}(t,r,z)\big{|}_{z=0}=T_{top},\;T_{s}(t,r,z)\big{|}_{z=H}=T_{bot},$ (9)

where $T_{top}$ is the specified soil temperature at the top of the heat
exchanger (assumed to be equal to the average annual temperature of the
ambient air), and $T_{bot}$ is calculated by the equation

$T_{bot}=T_{top}+G_{t}H,$

where $G_{t}$ is the geothermal gradient

As the initial conditions for equations (1)-(3) we put

$T_{1}(0,z)=T_{2}(0,z)=T_{top}+G_{t}z,$ (10)

and also

$T_{s}(0,r,z)=T_{top}+G_{t}z,$ (11)

###  2.3 Numerical scheme

We will solve problem (1)-(11) numerically using the finite difference method.
For this purpose, we introduce a uniform grid along the $z$ coordinate,
$z_{i}=ih$, $i=0,1,\dots,N_{z}$, $h=H/N_{z}$; and an exponential grid along
the $r$ coordinate,

$r_{j}=r_{b}\exp\left[\frac{j-0.5}{N_{r}}\ln\left(\frac{r_{d}}{r_{b}}\right)\right],\,j=1,2,\dots,N_{r}.$

An exponential grid has been used in several papers cited above and is needed
to record the rapid change in soil temperature near borehole. Obviously, it is
inspired by the fundamental solution of the Laplace equation in polar
coordinates, $T(r)\sim\ln r$.

Passing to grid functions, $T^{1}_{i}(t)=T_{1}(t,z_{i})$,
$T^{2}_{i}(t)=T_{2}(t,z_{i})$, $T^{s}_{i,j}(t)=T_{s}(t,r_{j},z_{i})$; and
replacing spatial derivatives in equations (1)-(3) with finite differences, we
obtain a system of $N_{z}(N_{r}+2)$ ordinary differential equations.

$\displaystyle
c_{f}\rho_{f}a_{1}\frac{dT^{1}_{i}}{dt}=c_{f}v_{f}\frac{T^{1}_{i+1}-T^{1}_{i}}{h}-q_{1}(T^{1}_{i},T^{2}_{i}),$
(12) $\displaystyle i=0,1,\dots,N_{z}-1,\,T^{1}_{N_{z}}=T^{2}_{N_{z}};$
$\displaystyle
c_{f}\rho_{f}a_{2}\frac{dT^{2}_{i}}{dt}=-c_{f}v_{f}\frac{T^{2}_{i}-T^{2}_{i-1}}{h}+q_{1}(T^{1}_{i},T^{2}_{i})-q_{2}(T^{2}_{i},T^{s}_{i,1}),$
(13) $\displaystyle i=1,2,\dots,N_{z},\,T^{2}_{0}=T_{in};$
$\displaystyle(c_{s}\rho_{s})_{i,j}S_{j}\frac{dT^{s}_{i,j}}{dt}=R_{i,j}+Z_{i,j},\,i=1,2,...,N_{z},\,j=1,2,...,N_{r};$
(14)

where $S_{j}=\pi(r^{*2}_{j}-r^{*2}_{j-1})$. $r^{*}_{0}=r_{b}$,
$r^{*}_{j}=0.5(r_{j}+r_{j+1})$.

$\displaystyle R_{i,1}=q_{2}(T^{2}_{i},T^{s}_{i,1})-2\pi
r^{*}_{1}\lambda_{i,1}^{r}\frac{T^{s}_{i,1}-T^{s}_{i,2}}{r_{2}-r_{1}};$
$\displaystyle R_{i,j}=2\pi
r^{*}_{j-1}\lambda_{i,j-1}^{r}\frac{T^{s}_{i,j-1}-T^{s}_{i,j}}{r_{j}-r_{j-1}}-2\pi
r^{*}_{j}\lambda_{i,j}^{r}\frac{T^{s}_{i,j}-T^{s}_{i,j+1}}{r_{j+1}-r_{j}},$
(15) $\displaystyle j=2,3,...,N_{r}-1;$ $\displaystyle R_{i,N_{r}}=2\pi
r^{*}_{N_{r}-1}\lambda_{i,N_{r}-1}^{r}\frac{T^{s}_{i,N_{r}-1}-T^{s}_{i,N_{r}}}{r_{N_{r}}-r_{N_{r}-1}},$
$\displaystyle i=1,2,...,N_{z},$

where

$\lambda_{i,j}^{r}=\frac{2\lambda_{i,j}\lambda_{i,j+1}}{\lambda_{i,j}+\lambda_{i,j+1}}.$
$\displaystyle
Z_{i,j}=\frac{S_{j}}{h^{2}}\left[\lambda_{i-1,j}^{z}(T^{s}_{i-1,j}-T^{s}_{i,j})-\lambda_{i,j}^{z}(T^{s}_{i,j}-T^{s}_{i+1,j})\right],$
(16) $\displaystyle i=1,2,...,N_{z},\,j=1,2,...,N_{r};$

where

$\lambda_{i,j}^{z}=\frac{2\lambda_{i,j}\lambda_{i+1,j}}{\lambda_{i,j}+\lambda_{i+1,j}}.\;T^{s}_{0,j}=T_{top},\,T^{s}_{N_{z}+1,j}=T_{bot}.$

We have written finite-difference equations for the general case when soil
parameters depend on spatial coordinates. Obviously, for a homogeneous soil,
the equations are somewhat simplified.

##  3 Model parameters

For simulations that demonstrate the capabilities of the model, it is
necessary to set the values of the parameters of the heat exchanger, heat-
carrying fluid and ground. The parameters of shallow ground heat exchangers
can take values in fairly wide ranges [33]. In this work, we do not consider
any special heat exchanger. Therefore, for numerical simulation, we have
chosen the typical values of the parameters that are found in the literature.
These values are presented in Table 1.

Table 1: Heat exchanger parameters | Description | Value
---|---|---
$H$ | Pipes length | 100 m
$r_{11}$ | Inner radius of inner pipe | 0.04 m
$r_{12}$ | Outer radius of inner pipe | 0.043 m
$r_{21}$ | Inner radius of outer pipe | 0.074 m
$r_{22}$ | Outer radius of outer pipe | 0.08 m
$r_{b}$ | Radius borehole | 0.1 m
$\lambda_{1}$ | Thermal conductivity of inner pipe | 0.4 W/(m$\cdot$K)
$\lambda_{2}$ | Thermal conductivity of outer pipe | 40 W/(m$\cdot$K)
$\lambda_{g}$ | Thermal conductivity of grout | 1.5 W/(m$\cdot$K)
$\alpha_{11}$ | Convective coefficient at $r_{11}$ | 1000 W/(m${}^{2}\cdot$K)
$\alpha_{12}$ | Convective coefficient at $r_{12}$ | 500 W/(m${}^{2}\cdot$K)
$\alpha_{21}$ | Convective coefficient at $r_{21}$ | 500 W/(m${}^{2}\cdot$K)
$c_{f}$ | Fluid heat capacity | 4500 J/(kg$\cdot$K)
$\rho_{f}$ | Fluid density | 1000 kg/m3
$v_{f}$ | Fluid mass velocity | 0.3 kg/s

Difference equations (12)-(16) make it possible to calculate the temperature
distribution and heat fluxes in inhomogeneous soil. To demonstrate this
possibility in this work, we have chosen an imaginary soil composed of several
horizontal layers of minerals with different physical properties. The
composition of the soil and the location of minerals in depth are presented in
Table 2. The properties of minerals are described by the parameters,
$\rho_{s}^{0}$, $c_{s}^{0}$, $\lambda_{s}^{0}$ – density, heat capacity and
thermal conductivity of dry matter, and volumetric moisture content in the
substance of the layer, $m$. The bottom part of the table shows the
thermophysical parameters for water and ice. In the next section, the
parameters of water are denoted by the subscript $w$ ($\rho_{w}$, $c_{w}$,
$\lambda_{w}$), and of ice – by the subscript $ice$ ($\rho_{ice}$, $c_{ice}$,
$\lambda_{ice}$).

Table 2: Strati-graphic column of borehole Soil type | $z_{t}-z_{b}$ | $\rho_{s}^{0}$ | $c_{s}^{0}$ | $\lambda_{s}^{0}$ | $m$
---|---|---|---|---|---
| m | kg/m3 | J/(m$\cdot$K) | W/(m$\cdot$K) | m3/m3
clay | 0 - 5 | 1700 | 920 | 1.1 | 0.2
limestone | 5 - 60 | 2500 | 840 | 1.0 | 0.15
mudstone | 60 - 90 | 2600 | 800 | 1.8 | 0.1
granite | 90 - 100 | 2700 | 790 | 1.1 | 0.05
Average | 0 - 100 | 2510 | 827 | 1.25 | 0.13
water |  | 997 | 4200 | 0.57 |
ice |  | 919 | 2108 | 2.25 |

##  4 Simulation results

Equations (12)-(14) can be written in matrix form as

$\frac{d{\bm{T}}}{dt}={\bf A}{\bm{T}}+{\bm{R}},$ (17)

where $\bm{T}$ is a column vector containing $N_{e}=(N_{r}+2)N_{z}$
components.

$T_{i}=T^{1}_{i},\;T_{N_{z}+i}=T^{2}_{i},\;i=1,\dots,N_{z};$

and

$\\{T_{(j+1)N_{z}+i}=T^{s}_{i,j},\;i=1,\dots,N_{z}\\},\;j=1,\dots,N_{r}.$

The matrix $\bf A$ depends on the fluid and soil parameters, while the vector
$\bm{R}$ takes into account the boundary conditions of the problem.The
elements of the matrix $\bf A$ and the column vector $\bm{R}$ are easily
determined from equations (12)-(14).

The matrix $\bf A$ contains $N_{e}^{2}$ elements, but it is essentially
5-diagonal. Namely, only $11N_{z}+(N_{r}-1)(5N_{z}-2)-5$ of them are nonzero.
With the values of $N_{z}$ and $N_{r}$ taken in our calculations, the ratio of
the number of non-zero elements to the total number of elements is of the
order of $10^{-4}$. Thus, the matrix $\bf A$ is very strongly sparse. In this
case, when solving system (17), it is very efficient to use the technique of
sparse matrices, which makes it possible to radically reduce the size of RAM
and the computation time. In this work, the solution of system (17) was
carried out in MatLab using the ode15s function, which effectively uses the
sparseness of the matrix and overcomes the stiffness of the problem.

During the operation of the heat exchanger, the temperature of the fluid and
the surrounding soil can drop below the freezing temperature of the ground
moisture, $T_{f}$ (for definiteness, we assume $T_{f}=0^{\circ}$C). Antifreeze
can be used as a heat carried fluid, but ground moisture will freeze at low
temperatures, and then, when the temperature rises, melt. Thus, in ground
moisture, a water-ice phase transition will occur. A change in the state of
aggregation of moisture leads to a change in the physical properties of the
soil and is accompanied by the release or absorption of the latent heat of the
phase transition, $L$ = 334 kJ/kg. Therefore, adequate modeling of heat
exchangers operating at the phase transition temperature must take this
transition into account. However, due to computational difficulties, it is
often ignored. In this work, for comparison, we present the results of
calculations obtained both with and without taking into account the freezing-
thawing of ground moisture.

In calculations that do not take into account the freezing of moisture, we
assume that moisture behaves like antifreeze. In the following, a model with
this assumption will be referred to as model A. In model A, at any
temperature, the ground parameters in any stratum are calculated by the
formula

$p_{s}=(1-m)p^{0}_{s}+mp_{w},$ (18)

where $p$ denotes $\rho$, $c$ or $\lambda$.

To model the freezing-thawing of ground moisture, we use the method of
apparent heat capacity [34, 35]. In this method, it is assumed that the water-
ice transition occurs in a fairly narrow temperature range,
$T_{f}-\Delta_{f}<T<T_{f}+\Delta_{f}$. When the temperature $T$ lies in this
interval, then it is assumed that water and ice are simultaneously present in
the soil. Moreover, the proportion of water, $w$, in the mixture is calculated
by the formula

$w=\frac{1}{2}\left(\frac{T}{\Delta_{f}}+1\right).$ (19)

And the ground parameters are determined by the equations

$\displaystyle\rho_{s}=(1-m)\rho^{0}_{s}+m[w\rho_{w}+(1-w)\rho_{ice}],$
$\displaystyle c_{s}=(1-m)c^{0}_{s}+m[wc_{w}+(1-w)c_{ice}+L/(2\Delta_{f})],$
(20)
$\displaystyle\lambda_{s}=(1-m)\lambda^{0}_{s}+m[w\lambda_{w}+(1-w)\lambda_{ice}].$

At $T>T_{f}+\Delta_{f}$ the soil parameters are calculated by Eq. 18. While at
$T<T_{f}-\Delta_{f}$, when all the moisture in the soil is frozen, the soil
parameters are calculated by the equation

$p_{s}=(1-m)p^{0}_{s}+mp_{ice},$ (21)

Allowance for moisture freezing makes system (17) non-linear and significantly
increases the computation time. When solving (17) numerically, the matrix
elements are recalculated only in the vicinity of the phase transition band.
As $\Delta_{f}$ decreases, the number of elements to be updated decreases too,
but the stiffness of the problem increases and the computation time increases.
Therefore, the value of $\Delta_{f}$ is chosen so as to balance these opposite
tendencies. In the calculations presented below, it is assumed that
$\Delta_{f}=0.1^{\circ}$C. In the following, the model that takes into account
the freezing of soil moisture will be referred to as the model F.

###  4.1 Stationary mode

Under unchanged boundary conditions, the solution ${\bm{T}}(t,r,z)$ of system
(17) at $t\to\infty$ asymptotically approaches the stationary solution
${\bm{T}}_{st}(r,z)$, which satisfies the system of algebraic equations

${\bf A}{\bm{T}}_{st}+{\bm{R}}=0.$ (22)

Figure 2 shows graphs of the stationary temperature of the fluid and the soil
adjacent to the borehole wall, calculated both with allowance for moisture
freezing (solid lines) and without freezing (dotted lines).

Figure 2: Temperature of incoming, $T_{2}(z)$, outgoing, $T_{1}(z)$, fluid and
adjacent soil, $T_{s,b}(z)$, in stationary mode at $T_{in}=-10^{o}$C, $v_{f}$
= 0.3 kg/s. Solid lines – model F, dotted lines – model A.

Freezing of ground moisture in our model leads to an increase in the thermal
conductivity of the soil, and hence to an increase in the heat flux from the
external soil to the heat exchanger. Therefore, in Fig. 2 temperature graphs
for the model with freezing are shifted to the right relative to the graphs
for the model without freezing. In particular, under given conditions, the
temperature of the fluid at the outlet of the heat exchanger
$T_{out}=-8.2^{\circ}$C in the model without freezing and $-8.0$ with
freezing. Under considered conditions, a difference of 0.2 degrees in
temperature means a difference of 972 kJ/h in heat extraction rate or 10
percent of the heat extraction rate in model A.

On the graphs of the temperature of the adjacent soil, a protrusion is
distinguished at a depth of 60 to 90 m. This protrusion corresponds to a layer
in which the thermal conductivity coefficient is significantly higher than in
neighboring layers. Therefore, the heat flux from the outer soil in this layer
is greater than in the neighboring ones, and, accordingly, the soil
temperature is higher.

In our model and in many other models of heat exchangers, it is assumed that
the heat flow between the coolant and the adjacent soil is proportional to the
temperature difference between the soil and the coolant. Figure 3 shows the
graphs of the temperature difference between the incoming liquid and the
adjacent soil for model A and for model F.

Figure 3: The temperature difference between the adjacent soil and the
incoming fluid, $T_{s,b}(z)-T_{2}(z)$, at $T_{in}=-10^{o}$C, $v_{f}$ = 0.3
kg/s.

The color map of the temperature field in the vicinity of the borehole is
shown in Fig. 4, on the left – for model A and on the right – for model F. The
numbers near the temperature level lines show the temperature value on the
corresponding line. It can be seen that the soil temperature in the model with
freezing (F) is lower than the soil temperature in the model without freezing.
This is understandable, since it follows from Figures 2 and 3 that the amount
of heat extracted from the soil per unit time in model F is greater than in
model A.

Figure 4: Color map of soil temperature in stationary mode at $T_{in}$ =
-10∘C, $v_{f}$ = 0.3 kg/s. The number at the level line indicates the
temperature in ∘C. Left – model A, right – model F.

However, it also follows from Fig. 2 that the temperature of the soil adjacent
to the borehole in model F is higher than in model A. This apparent
contradiction is clarified by Fig. 5, which shows plots of soil temperature
depending on the radius for $z$ = 50 m (the middle of the exchanger). The
temperature graphs intersect at a distance of approximately 10 cm from the
borehole wall. To the left of the intersection point, the temperature in model
F is higher than in model A, and to the right, it is vice versa.

Figure 5: Soil temperature versus radius at a depth of 50 m in models A and F.
$T_{in}$ = -10∘C, $v_{f}$ = 0.3 kg/s.

In stationary mode, it is convenient to evaluate the sensitivity of various
characteristics of the heat exchanger to changes in its parameters. In
particular, let us estimate the sensitivity of the stationary heat extraction
rate, $q_{st}$. In a linear approximation, the change in $q_{st}$ with a
change in the parameter $p$ is calculated by the formula

$\Delta q_{st}\approx\frac{\partial q_{st}}{\partial p}\Delta p,$

where $p$ is one of the GHE parameters listed in Table 1. We estimate the
partial derivatives numerically using the central finite difference by
calculating $q_{st}$ at $p=(1\pm 0.01)p_{0}$, where $p_{0}$ is the value of
parameter $p$ specified in the table. The results of calculation are presented
in Table 3.

Table 3: Parameters sensitivity |  $T_{in}=+2^{\circ}$C |  $T_{in}=-10^{\circ}$C, A |  $T_{in}=-10^{\circ}$C, F
---|---|---|---
$r_{11}$ |  $-$1.732e+04 |  $-$3.758e+04 |  $-$4.582e+04
$r_{12}$ | 1.568e+04 | 3.401e+04 | 4.148e+04
$r_{21}$ | 7.572e+01 | 1.672e+02 | 2.825e+02
$r_{22}$ | 1.868e+03 | 4.125e+03 | 5.444e+03
$r_{b}$ | 2.985e+02 | 7.039e+02 | 1.517e+01
$\lambda_{1}$ |  $-$1.231e+02 |  $-$3.563e+01 |  $-$3.293e+02
$\lambda_{2}$ | 1.092e$-$02 | 2.412e$-$02 | 2.833e$-$02
$\lambda_{g}$ | 2.223e+01 | 4.909e+01 | 6.518e+01
$\alpha_{11}$ |  $-$6.808e$-$03 |  $-$1.477e$-$02 |  $-$1.786e$-$02
$\alpha_{12}$ |  $-$2.533e$-$02 |  $-$5.495e$-$02 |  $-$6.507e$-$02
$\alpha_{21}$ | 1.121e$-$02 | 2.475e$-$02 | 4.180e$-$02
$c_{f}$ | 3.964e$-$02 | 8.725e$-$02 | 1.094e$-$01
$\rho_{f}$ | 1.137e$-$12 |  $-$1.439e$-$11 |  $-$1.105e$-$11
$v_{f}$ | 5.946e+02 | 1.309e+03 | 1.642e+03
$T_{in}$ |  $-$1.124e+02 |  $-$1.242e+02 |  $-$1.382e+02

Table 3 shows estimates of partial derivatives for both positive,
$T_{in}=2^{\circ}$C, and negative, $T_{in}=-10^{\circ}$C, inlet temperatures.
Moreover, for a negative temperature, the derivatives were calculated in both
models A and F. In all calculations for the heat extraction rate $q_{st}$,
$r_{11}$ and $r_{12}$ turned out to be the most sensitive parameters. The
derivatives with respect to them have different signs. The derivative with
respect to $r_{11}$ is negative. An increase in $r_{11}$ leads to a decrease
in the thermal resistance of the inner pipe wall and to a decrease in the
linear velocity of the fluid in the inner pipe. Therefore, the fluid cools
down more strongly when moving through the inner pipe, which leads to a
decrease in $q_{st}$. The derivative with respect to $r_{12}$ is positive. An
increase in $r_{12}$ leads to an increase in thermal resistance between the
fluids in the outer and inner pipes, so the outgoing fluid cools less and
$q_{st}$ increases. Similar conclusions about the effect of inlet temperature
and inner and outer pipes diameters on the thermal performance of a coaxial
GHE were reported in [31] for a deep GHE.

Note that under given conditions, all derivatives of $q_{st}$ with respect to
the heat exchanger parameters have the same signs and almost all of them are
of the same order of magnitude. The exceptions are the derivatives with
respect to $r_{b}$ and $\lambda_{1}$. The derivative with respect to $r_{b}$
in model A is 46 times greater than in model F, and the derivative with
respect to $\lambda_{1}$ in model A is approximately 10 times smaller in
absolute value than in model F.

The stationary solution is easy to calculate and provides useful information
about the heat extraction process. However, with a significant difference
between the average soil temperature (11.5∘C in our example) and the inlet
temperature $T_{in}$, the stationary solution turns out to be practically
unattainable due to the large relaxation time. Figure 6 gives an idea of the
relaxation time for the assumed conditions.

Figure 6: Estimation of the duration of the relaxation time of soil and fluid
depending on the inlet temperature.

The relaxation time increases rapidly as the difference between the natural
ground temperature and the inlet fluid temperature increases. Moreover, the
temperature of the fluid in the heat exchanger approaches the equilibrium
temperature much faster than the temperature of the soil. We estimate the time
of approaching the soil temperature to equilibrium by the time $t_{*}$ at
which the inequality $\max_{r,z}|T(t,r,z)-T_{st}|<\delta$ is satisfied for the
first time, where $T$ is the solution of (17) and $\delta=1^{\circ}$C. For the
fluid, we estimate the time $t_{\circ}$ at which the inequality
$T_{out}(t)-T_{out}^{st}<\delta$ begins to hold for $\delta=0.1^{\circ}$C. The
graphs of $t_{*}$ and $t_{\circ}$ in dependence on $T_{in}$ are shown in
Figure 6.

###  4.2 Transient mode

Let us now return to system (17) and present the time-dependent solution of
the system immediately after the start of working. Here and below, all
calculations are performed with the initial conditions (10), (11).

Figure 7 shows graphs of fluid temperature, $T_{out}(t)$, at the outlet of the
heat exchanger for the first five days of operation at $T_{in}=-10^{\circ}$C
for both models A and F.

Figure 7: Outlet temperature $T_{out}(t)$ in model A and F.
$T_{in}=-10^{\circ}$C, $v_{f}$ = 0.3 kg/s.

During the first day of operation, the temperature $T_{out}(t)$ drops rapidly.
By the end of the fifth day, the temperature drop noticeably slows down,
although it is still very far from the equilibrium temperature,
$-8.2^{\circ}$C for model A and $-8.0^{\circ}$C for model F (see Fig. 2). The
convergence of the solution of both models to the stationary solution is very
slow. From Figure 6, we can conclude that the relaxation time for
$T_{in}=-10^{\circ}$C must be very large and thermal equilibrium in this
example is practically unattainable.

By the end of the fifth day the outlet temperature in model A equals to
$-5.85^{\circ}$C and in model F it equals $-5.10^{\circ}$C. The heat
extraction rate $q_{ex}$ is proportional to the difference $T_{out}-T_{in}$.
Thus, under given conditions, by the end of the fifth day of work, $q_{ex}$
calculated using model F is 18% more than $q_{ex}$ calculated using model A.

The evolution of soil temperature near the borehole in the beginning five days
of operation is shown in Figures 8 and 9. These figures show temperature
versus time graphs at nine spatial grid nodes with coordinates $z$ = 0.25 m
and $r$ = 0.106, 0.120, 0.136, 0.154, 0.175, 0.198, 0.224, 0.254, 0.288 m,
from bottom to top (borehole radius $r_{d}$ = 0.1 m). These nodes are located
in horizontal plane just near the top of the exchanger ($z$ = 0).

Figure 8 shows plots of soil temperature in model A, obtained without taking
into account moisture freezing. Here the temperature changes smoothly,

Figure 8: Soil temperature in model A at 9 grid nodes (see text). $T_{in}$ =
-10∘C, $v_{f}$ = 0.3 kg/s. Figure 9: Soil temperature in model F at the same 9
grid nodes as in Fig. 8. $T_{in}$ = -10∘C, $v_{f}$ = 0.3 kg/s.

decreasing rapidly on the first day and slowing down noticeably at the end of
the fifth day. Note that the temperature curves are almost equidistant. This
fact is a consequence of the exponential grid in the radial direction.

A different view have the soil temperature curves at the same grid nodes in
Figure 9, calculated using the model F. When the temperature in the node and
in the corresponding cell of the spatial grid decreases to $T_{f}-\Delta_{f}$,
moisture begins to freeze in the cell. During freezing, latent heat $L$ is
released, which maintains the temperature in the interval
($T_{f}-\Delta_{f},\,T_{f}+\Delta_{f})$ in the grid cell and slows down the
temperature change in neighboring cells in the radial direction. This
phenomenon is known as the zero-curtain effect [36].

The farther the node is from the borehole, the larger the volume of the cell
and the more latent heat is released. Therefore, with the distance from the
borehole, the freezing time of moisture in the cell increases. Temperature
stabilization in the cell for the time of moisture freezing slows down the
temperature change not only in the cells with a negative temperature located
closer to the borehole wall, but also in several adjacent cells with a
positive temperature, which lie farther from the borehole in the same
horizontal plane. The release of the latent heat of freezing increases the
heat flux to the heat exchanger and increases the temperature of the fluid in
the outer pipe, which is confirmed by the graphs in Figure 7.

Obviously, the stepped form of the curves is a consequence of the
discretization of the problem. With refinement of the spatial grid, the number
of steps will increase, their height will decrease, and the curves will be
smoothed out.

Figure 10 shows color maps of soil temperature near the borehole – in the area
of the most rapid temperature change with distance from the borehole. The
general slope of the level lines to the right is due to the geothermal
gradient. Soil temperature at $z$ = 100 (13∘C) is three degrees higher than at
$z$ = 0. The breaks in the temperature level lines on the left map (model A)
correspond to the boundary between the soil layers with different physical
properties. These kinks are also visible on the right map (model F), but here
one can also see much more kinks associated with temperature stabilization
when moisture freezes. Level lines corresponding to the same temperature in
model A are located farther from the borehole than in model F. Therefore, at
the end of the fifth day, the amount of heat in the soil calculated by model A
is less than in model F. There is no contradiction here with Fig. 4, because
in Fig. 4 shows the stationary (final) temperature distribution in the soil,
and Fig. 10 shows the unsteady transitional distribution.

Figure 10: Color map of soil temperature near the well at the end of the fifth
day of heat exchanger operation. Left – model A, right – model F. $T_{in}$ =
-10oC, $v_{f}$ = 0.3 kg/s.

###  4.3 Controlled mode

And, finally, we give an example of a controlled mode of operation of the heat
exchanger. Suppose that during the heating season (October 1 - April 31) every
day the heat exchanger must eject from the soil the amount of heat $Q_{d}$,
proportional to the difference between the average daily ambient temperature
per day $d$, $T_{a}^{d}$, and room temperature $T_{r}=20^{\circ}$C,

$Q_{d}=k_{q}(T_{a}^{d}-T_{r}),$

where $k_{q}$ [J/K] is a constant coefficient.

For definiteness, let’s take the Moscow region. The graph of the average daily
air temperature in Moscow is shown in Fig. 11.

Figure 11: Average daily temperature in Moscow. First date – October 1st.

Let the temperature of the fluid at the inlet to the heat exchanger $T_{in}$
during the $d$-th day be constant and equal to $T_{in}^{d}$, then the heat
exchanger takes away on the $d$-th day the amount of heat $Q_{ex}^{d}$ equal
to

$Q_{ex}^{d}(T_{in}^{d})=c_{f}v_{f}\int_{d}{(T_{out}(t)-T_{in}^{d})\,dt},$

where $T_{out}(t)$ is obtained from the solution of Eq. (17) with
$T_{in}=T_{in}^{d}$ and the temperature field at the end of the previous day
as the initial condition. The integral is taken over $d$-th day. To obtain the
required amount of heat on each day of the heating season, it is necessary to
solve the following chain of equations for $T_{in}^{d}$,

$Q_{ex}^{d}(T_{in}^{d})=Q_{d},\;d=1,2,...$ (23)

Figure 12 shows the results of calculations of the heating period, performed
according to the F and A models for two values of the coefficient $k_{q}$, 104
and $2\times 10^{4}$. The top panel shows plots of extracted heat $Q_{ex}(t)$.
On the scale of the Figure, the graphs of $Q_{ex}(t)$ and the required heat
$Q_{d}(t)$ coincide. The lower panel shows graphs of the inlet temperature of
the fluid $T_{in}$ corresponding to the required amount of heat $Q_{d}$,
calculated by model A and model F.

Figure 12: Extracted heat (top) and inlet temperature (bottom) in controlled
mode.

From Fig. 12 it follows that for both considered values of $k_{q}$, negative
fluid temperatures are necessary to satisfy the required amount of heat
extracted from the soil. Note that if at $k_{q}=10^{4}$ the difference in the
inlet temperatures of the fluid in the heat exchanger, calculated according to
models A and F, is small, then at $k_{q}=2\times 10^{4}$ the maximum
difference reaches 37 percent: $-23.9^{\circ}$C in model A and
$-17.4^{\circ}$C in model F.

The numerical solution of equation (23) for the $d$-th day was obtained by a
method similar to the method of chords for solving nonlinear equations. First,
the amount of heat $Q_{1}$ extracted at the inlet temperature
$T_{1}=T_{in}^{d-1}$ is calculated. Then the amount of heat $Q_{2}$ extracted
from the soil at the inlet temperature $T_{2}=T_{1}+\delta$, if $Q_{1}>Q_{d}$,
or at $T_{2}=T_{1}-\delta$, if $Q_{1}<Q_{d}$ is calculated. Then

$T_{3}=T_{1}+\frac{Q_{d}-Q_{1}}{Q_{2}-Q_{1}}(T_{2}-T_{1})$

was assigned and the amount of heat $Q_{3}$ extracted at the inlet temperature
$T_{3}$ was calculated. The described iterations stopped when the condition

$|Q_{3}-Q_{d}|/Q_{d}<0.001.$

Note that in all our calculations, this condition was already satisfied at the
first iteration.

The temperature field in the soil in the vicinity of the borehole at the end
of the heating period is shown in Fig. 13. The temperature is calculated
according to the F model, i.e. taking into account the freezing/thawing of
moisture, on the left – for $k_{q}=10^{4}$, on the right – for $k_{q}=2\times
10^{4}$. At $k_{q}=10^{4}$, the ice has completely melted by the end of the
heating period, and the ground temperature is everywhere positive. At
$k_{q}=2\times 10^{4}$, a significant part of the moisture in the soil remains
frozen. In the upper part of the borehole, the maximum ice boundary
approximately runs at a distance of 0.78 m from the borehole axis.

Figure 13: Soil temperature at the end of heating season (model F). Left –
$k_{q}=10^{4}$ J/K, right – $k_{q}=2\times 10^{4}$ J/K.

Thus, at the end of the heating season, in both considered cases, the soil
temperature near the borehole is significantly lower than the undisturbed soil
temperature. Therefore, the question arises of restoring the soil temperature
by the beginning of the next heating season. In this paper, we do not consider
this issue, but note that the natural temperature recovery due to the
influence of boundary conditions cannot be correctly described by our model.
Our model does not take into account the convective heat transfer in fluids in
the vertical direction. When fluid moves in a working heat exchanger, this
transfer can be ignored. However, in a stagnant fluid in the summer, it can
have a significant effect on the restoration of the soil temperature near the
borehole, in particular, on the melting of adjacent ice.

##  5 Conclusion

The paper proposes a new mathematical model of heat transfer in a shallow
vertical coaxial ground heat exchanger and in the surrounding ground. The
model takes into account soil heterogeneity and geothermal gradient. At a
negative temperature of the heat-carrying fluid, the process of freezing-
thawing of ground moisture is taken into account. The water-ice phase
transition is modeled using the apparent heat capacity method. The numerical
implementation of the model was performed using the finite difference method
on a non-homogeneous spatial grid under the assumption of axial symmetry of
the problem. The computer implementation of the model is made in the Matlab
environment on a personal computer, which makes the model accessible to a wide
range of specialists.

In the considered examples, three modes of operation of the heat exchanger are
modeled: stationary, transitional and controlled. In the calculations, the
main attention was paid to demonstrating the differences in the results
obtained with and without taking into account the water-ice phase transition
in ground moisture. From the above examples, it follows that the greatest
difference is noted in the controlled extraction of heat from the soil during
the heating season.

The proposed model will be useful in the design of heat pump installations for
heating and cooling buildings if it is necessary to work with negative
temperatures of the coolant, taking into account the structural features and
composition of the soil. To simulate the release of heat into the soil, it is
enough to change the direction of fluid movement in the model. The model can
also be used to optimize the parameters of the heat exchanger, taking into
account the specific conditions of its operation.

## Disclosure of Potential Conflicts of Interest

The Authors declare that there is no conflict of interest.

## Acknowledgment

## References

  * [1] D. Menegazzo, G. Lombardo, S. Bobbo, M. De Carli and L. Fedele. State of the art, perspective and obstacles of ground-source heat pump technology in the European building sector: A review. Energies 15 (2022) 2685. DOI: 10.3390/en15072685
  * [2] Z. Zhao, Y.-F. Lin, A. Stumpf, X. Wang. Assessing impacts of groundwater on geothermal heat exchangers: A review of methodology and modeling. Renewable Energy 190 (2022) 121-147. DOI: 10.1016/j.renene.2022.03.089
  * [3] S. Cho and S. Mirianhosseinabadi. Simulation modeling of ground source heat pump systems for the performance analysis of residential buildings. Proceedings of the 13th Conference of International Building Performance Simulation Association, Chambéry, France, 2013, p. 1960-1967.
  * [4] J. Cao, M. Bottarelli, M. Bortoloni, G. Pei. Small-scale lab analysis of the ground freezing effect on the thermal performance of a Flat-Panel ground heat exchanger. Geothermics 74 (2018) 247-254. DOI: 10.1016/j.geothermics.2018.03.013
  * [5] P. Eslami-nejad, M. Bernier. Freezing of geothermal borehole surroundings: A numerical and experimental assessment with applications. Applied Energy 98 (2012) 333–345. DOI: 10.1016/j.apenergy.2012.03.047
  * [6] W. Yang, L. Kong, Y. Chen. Numerical evaluation on the effects of soil freezing on underground temperature variations of soil around ground heat exchangers. Applied Thermal Engineering 75 (2015) 259–269. DOI: 10.1016/j.applthermaleng.2014.09.049
  * [7] X. Meng, Z. Han, H. Hu, H. Zhang, X. Li. Studies on the performance of ground source heat pump affected by soil freezing under groundwater seepage. Journal of Building Engineering 33 (2021) 101632. DOI: 10.1016/j.jobe.2020.101632
  * [8] S. Tu, X. Zhang, X. Zhou. A revised thermal resistance and capacity model for the ground heat exchanger under freezing soil conditions and thermal performance analysis. Procedia Engineering 205 (2017) 19-26. DOI: 10.1016/j.proeng.2017.09.929
  * [9] J. Wu, D. Wu, W. Wang, H. Zhao. Research on thermal storage and release capacity of freezing sandy soil with phase change, Applied Thermal Engineering 166 (2020) 114638. DOI: 10.1016/j.applthermaleng.2019.114638
  * [10] H. Zhang, Z. Han, X. Li, M. Ji, X. Zhang, G. Li, L. Yang. Study on the influence of borehole spacing considering groundwater flow and freezing factors on the annual performance of the ground source heat pump, Applied Thermal Engineering 182 (2021) 116042. DOI: 10.1016/j.applthermaleng.2020.116042
  * [11] S. Erol, Experimental and analytical evaluations of freezing effects in borehole heat exchangers. Conference Proceedings, Second EAGE Workshop on Geomechanics and Energy, Oct 2015, cp-466-00023. DOI: 10.3997/2214-4609.201414309
  * [12] J.‑H. Kupfernagel, J.C. Hesse, M. Schedel, B. Welsch, H. Anbergen, L. Müller and I. Sass. Impact of operational temperature changes and freeze–thaw cycles on the hydraulic conductivity of borehole heat exchangers. Geotherm Energy (2021) 9:24. DOI: 10.1186/s40517-021-00206-y
  * [13] A. Kirschbaum, J.M. Kuckelkorn, K. Hagel. Vertical hydraulic conductivity of borehole heat exchanger systems before and after freeze-thaw cycle stress. Proceedings of the IGSHPA Research Track 2018. DOI: 10.22488/okstate.18.000032
  * [14] H. Anbergen, J. Frank, L. Müller, I. Sass. Freeze–thaw-cycles on borehole heat exchanger grouts: impact on the hydraulic properties. Geotech Test Journal 37 (2014) 639–651. DOI: 10.1520/ GTJ20 130072
  * [15] H. Anbergen, W. Rühaak, J. Frank, I. Sass/ Numerical simulation of a freeze–thaw testing procedure for borehole heat exchanger grouts. Canadian Geotech. Journal 52 (2015) 1087–1100. DOI: 10.1139/cgj-2014-0177
  * [16] Y. Wang, Q. Gao, X. Zhu, M. Yu, X. Zhao. Experimental study on interaction between soil and ground heat exchange pipe at low temperature, Applied Thermal Engineering 60 (2013) 137–144. DOI: 10.1016/j.applthermaleng.2013.06.053
  * [17] Y. Wang, Y. Liu, Y. Cui, W. Guo, J. Lv, Numerical simulation of soil freezing and associated pipe deformation in ground heat exchangers. Geothermics 74 (2018) 112-120. DOI: 10.1016/j.geothermics.2018.02.010
  * [18] G. Dalla Santa, Z. Farina, H. Anbergen, W. Rühaak, A. Galgaro. Relevance of computing freeze-thaw effects for borehole heat exchanger modelling: A comparative case study, Geothermics 79 (2019) 164-175. DOI: 10.1016/j.geothermics.2019.02.001
  * [19] T. Zheng, H. Shao, S. Schelenz, P. Hein, T. Vienken, Z. Pang, O. Kolditz, T. Nagel. Efficiency and economic analysis of utilizing latent heat from groundwater freezing in the context of borehole heat exchanger coupled ground source heat pump systems, Applied Thermal Engineering 105 (2016) 314-326. DOI: 10.1016/j.applthermaleng.2016.05.158
  * [20] H.-J. Diersch, D. Bauer, W. Heidemann, W. Rühaak, P. Schätzl. Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals. Computers & Geosciences 37 (2011) 1122–1135. DOI: 10.1016/j.cageo.2010.08.003
  * [21] T. Yang, X. Zhang, B. Zhou, M. Zheng. Simulation and experimental validation of soil cool storage with seasonal natural energy. Energy and Buildings 63 (2013) 98–107. DOI: 10.1016/j.enbuild.2013.03.019
  * [22] H. Zhang, Z. Han, G. Li, M. Ji, X. Cheng, X. Li, L. Yang. Study on the influence of pipe spacing on the annual performance of ground source heat pumps considering the factors of heat and moisture transfer, seepage and freezing. Renewable Energy 163 (2021) 262-275. DOI: 10.1016/j.renene.2020.08.149
  * [23] Z. Li, J. Chen, M. Sugimoto, H. Ge. Numerical simulation model of artificial ground freezing for tunneling under seepage flow conditions. Tunneling and Underground Space Technology 92 (2019) 103035. DOI: 10.1016/j.tust.2019.103035
  * [24] S. Tu, X. Yang, X. Zhou, M. Luo, X. Zhang. Experimenting and modeling thermal performance of ground heat exchanger under freezing soil conditions. Sustainability 11 (2019) 5738. DOI: 10.3390/su11205738
  * [25] J. Acuña, B. Palm. A novel coaxial borehole heat exchanger: description and first distributed thermal response test measurements. In Proceedings of the World Geothermal Congress, 2010, April (p. 7).
https://www.geothermal-energy.org/pdf/IGAstandard/WGC/2010/2953.pdf

  * [26] P. Li, P. Guan, J. Zheng, B. Dou, H. Tian, X. Duan, H. Liu. Field test and numerical simulation on heat transfer performance of coaxial borehole heat exchanger. Energies 13 (2020) 5471. DOI: 10.3390/en13205471
  * [27] A. Zarrella, G. Emmi, S. Graci, M. De Carli, M. Cultrera, G.D. Santa, A. Galgaro, D. Bertermann, J. Müller, L. Pockelé, G. Mezzasalma, D. Righini, M. Psyk, A. Bernardi. Thermal response testing results of different types of borehole heat exchangers: An analysis and comparison of interpretation methods. Energies 10 (2017) 801. DOI: 10.3390/en10060801
  * [28] C. Wang, Y. Lu, L. Chen, Z. Huang, H. Fang. A semi-analytical model for heat transfer in coaxial borehole heat exchangers. Geothermics 89 (2021)10195. DOI: j.geothermics.2020.101952
  * [29] Y. Luo, H. Guo, F. Meggers, L. Zhang. Deep coaxial borehole heat exchanger: Analytical modeling and thermal analysis. Energy 185 (2019), 1298-1313. DOI: 10.1016/j.energy.2019.05.228
  * [30] C. Wang, H. Fang, X. Wang, J. Lu, Y. Sun. Study on the influence of borehole heat capacity on deep coaxial borehole heat exchanger. Sustainability 14 (2022) 2043. DOI: 10.3390/su14042043
  * [31] J. Liu, F. Wang, W. Cai, Z. Wang, Q. Wei, J. Deng. Numerical study on the effects of design parameters on the heat transfer performance of coaxial deep borehole heat exchanger. Int J Energy Res. (2019) 1–16. DOI: 10.1002/er.4357
  * [32] W. Cai, F. Wang, J. Jiang, Z. Wang, J. Liu, C. Chen. Long-term performance evaluation and economic analysis for deep borehole heat exchanger heating system in Weihe basin. Front. Earth Sci. 10 (2022) 806416. DOI: 110.3389/feart.2022.806416
  * [33] A. Shah, M. Krarti, J. Huang. Energy performance evaluation of shallow ground source heat pumps for residential buildings. Energies 15 (2022)1025. DOI:10.3390/en15031025
  * [34] H.T. Hashemi and C.M. Sliepcevich. A numerical method for solving two-dimensional problems of heat conduction with change of phase. Chem. Eng. Prog. Symp. Series 63(3) (1967) 34–41.
  * [35] H. Hu, S. A Argyropoulos. Mathematical modelling of solidification and melting: a review. Modeling Simul. Mater. Sci. Eng. 4 (1996) 371–396. DOI: 10.1088/0965-0393/4/4/004
  * [36] S.I. Outcalt, K.M. Hinkel. The zero-curtain effect – heat and mass transfer across an isothermal region in freezing soil. Water Resources Research 26 (1990) 1509-1516. DOI: 10.1029/wr026i007p01509


# A Locking-free modified conforming FEM for planar elasticity††thanks: School
of Mathematics and Statistics, University of New South Wales, Sydney,
Australia††thanks: This work was supported by the Australian Research Council
grant DP220101811.

K. Mustapha, W. McLean, J. Dick and Q. T. Le Gia

###### Abstract

Due to the divergence-instability, low-order conforming finite element methods
(FEMs) for nearly incompressible elasticity equations suffer from the so-
called locking phenomenon as the Lamé parameter $\lambda\to\infty$ and
consequently the material becomes more and more incompressible. For the
piecewise linear case, the error in the $L^{2}$-norm of the standard Galerkin
conforming FEM is bounded by $C_{\lambda}h^{2}$. However,
$C_{\lambda}\to\infty$ as $\lambda\to\infty$, resulting in poor accuracy for
practical values of $h$ if $\lambda$ is sufficiently large. In this short
paper, we show that for 2D problems the locking phenomenon can be controlled
by replacing $\lambda$ with $\lambda^{\alpha}$ in the stiffness matrix, for a
certain choice of $\alpha=\alpha_{*}(h,\lambda)$ in the range $0<\alpha\leq
1$. We prove that for this optimal choice of $\alpha$, the error in the
$L^{2}$-norm is bounded by $Ch$ where $C$ does not depend on $\lambda$.
Numerical experiments confirm the expected convergence behaviour and show
that, for practical meshes, our locking-free method is more accurate than the
standard method if the material is nearly incompressible. Our analysis also
shows that the error in the $H^{1}$-norm is bounded by $Ch^{1/2}$, but our
numerical experiments suggest that this bound is not sharp.

##  1 Introduction

Consider the following linear elasticity equation

$-\nabla\cdot\boldsymbol{\sigma}({\bf u}({\boldsymbol{x}}))={\bf
f}({\boldsymbol{x}})\quad\text{for
${\boldsymbol{x}}\in\Omega\subseteq\mathbb{R}^{2}$,}$ (1.1)

subject to homogeneous Dirichlet boundary conditions ${\bf u}={\bf 0}$ on the
polygonal boundary $\Gamma:=\partial\Omega$. The Cauchy stress tensor
$\boldsymbol{\sigma}\in[L^{2}(\Omega)]^{2\times 2}$ is defined as

$\boldsymbol{\sigma}({\bf u})=\lambda\nabla\cdot{\bf u}\,{\bf
I}+2\mu{\boldsymbol{\varepsilon}}({\bf u})\quad\text{on $\Omega$,}$

with ${\bf u}$ being the displacement vector field and the symmetric strain
tensor ${\boldsymbol{\varepsilon}}({\bf u}):=\frac{1}{2}(\nabla{\bf
u}+(\nabla{\bf u})^{T})$. Here, ${\bf f}$ is the body force per unit volume,
and ${\bf I}$ is the identity tensor. The gradient ($\nabla$) and the
divergence ($\nabla\cdot$) are understood to be with respect to the physical
variable ${\boldsymbol{x}}\in\Omega$. The constant Lamé parameters $\mu$ and
$\lambda$ are strictly positive, and we are interested in the case when the
Poisson ratio $\nu$ of the elastic material approaches $1/2$. Thus, $\lambda$
is large and we are dealing with a nearly incompressible material. It is known
in this case that locking is a major issue when conforming FEMs are employed
to approximate the solution of (1.1). It was shown that locking is absent for
polynomials of degree $\geq 4$ on a triangular mesh [28]. However, on
quadrilateral meshes, it was shown that locking cannot be avoided for any
polynomial of degree $\geq 1$ [4]. In a related result, Vogelius [31] proved
that the $p$-version of the conforming FEM on smooth domains is unaffected by
locking.

The locking phenomenon [1] in conforming FEMs has motivated researchers to
explore and develop various locking-free numerical methods ($h$-, $p$-, and
$hp$-versions, low-order, and high-order) for solving the nearly
incompressible elasticity equation (1.1). These methods include nonconforming
FEMs [3, 7, 9, 16, 23, 26, 34], conforming-nonconforming mixed FEMs [5, 17,
19, 35], the weak Galerkin FEM [10, 20, 22, 24, 32], the discontinuous
Petrov–Galerkin method [8], the discontinuous Galerkin method (with
hybridization) [11, 14, 18, 29, 33], and the virtual element method [15, 21,
30, 36, 37]. It is worth mentioning that most of these methods are more
complicated both computationally and in their analysis than the conforming
FEM. Furthermore, these methods become more complicated or not applicable for
the case of an inhomogeneous material, where the Lamé parameters are functions
of ${\boldsymbol{x}}$.

The aim of this work is to demonstrate that the locking phenomenon in the
popular piecewise linear conforming FEM can be avoided with a slight
modification of the scheme. More precisely, we replace the Lamé parameter
$\lambda$ with $\lambda^{\alpha}$ in the standard conforming FEM, where
$0<\alpha\leq 1$. The mechanism for choosing the optimal value
$\alpha=\alpha_{*}(h,\lambda)$ depends on the finite element mesh size $h$,
the Lamé parameter $\lambda$ and on diameter $d_{\Omega}$ of $\Omega$.
Although most of our analysis is independent of the dimension, we need to
assume that $\Omega\subseteq\mathbb{R}^{2}$ so that we can apply certain
$H^{2}$-regularity results of Brenner and Sung [7].

In the next section, we succinctly discuss the weak formulation of Problem
(1.1) and the well-posedness of the modified conforming FEM. In Section 3, we
demonstrate that the scheme is $O(h)$ accurate, uniformly in $\lambda$.
Additionally, we include a remark on the applicability of these results in the
case of mixed Dirichlet-traction boundary conditions. Furthermore, we include
another remark on how to extend the proposed finite element scheme and the
convergence analysis for the case of an inhomogeneous body, where the Lamé
parameters $\lambda$ and $\mu$ depend on the space variable
${\boldsymbol{x}}\in\Omega$. To support our theoretical findings, we present
some numerical results in Section 4 for two different examples, including
Cook’s benchmark problem.

##  2 Modified conforming FEM

In this section, we briefly discuss the existence and uniqueness of the weak
solution of (1.1), and then define our piecewise linear modified conforming
Galerkin FEM. We investigate the error estimates in a separate section.

In the weak formulation of equation (1.1), we seek ${\bf u}\in{\bf V}$
satisfying

$\mathcal{B}({\bf u},{\bf v})=\ell({\bf v}),\quad\text{for all~{}~{}${\bf
v}\in{\bf V}$,}$ (2.1)

where the linear functional $\ell$ and the bilinear form $\mathcal{B}$ are
defined by

$\ell({\bf v})=\int_{\Omega}{\bf f}\cdot{\bf
v}\,\mathrm{d}{\boldsymbol{x}}\quad\text{and}\quad\mathcal{B}({\bf u},{\bf
v})=\int_{\Omega}[2\mu\,{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf v})+\lambda(\nabla\cdot{\bf
u})(\nabla\cdot{\bf v})]\,\mathrm{d}{\boldsymbol{x}}.$

The colon operator is the inner product between tensors. Let $\|\cdot\|$
denote the norm in ${\bf L}^{2}(\Omega)=[L^{2}(\Omega)]^{2}$, and let
$\|\cdot\|_{\bf V}$ and $\|\cdot\|_{\bf H}$ denote the norms in ${\bf
V}=[H^{1}_{0}(\Omega)]^{2}$ and ${\bf H}=[H^{2}(\Omega)]^{2}$, respectively.
Here, $H^{1}_{0}(\Omega)$ and $H^{2}(\Omega)$ are the usual Sobolev spaces of
scalar-valued functions on $\Omega$. Finally, ${\bf V}^{*}$ is the dual space
of ${\bf V}$ with respect to the inner product in ${\bf L}^{2}(\Omega)$.

Since $\mathcal{B}({\bf v},{\bf w})\leq C_{\lambda,\mu}\|\nabla{\bf
v}\|\,\|\nabla{\bf w}\|\leq C_{\lambda,\mu}\|{\bf v}\|_{{\bf V}}\,\|{\bf
w}\|_{{\bf V}}$ for ${\bf v},{\bf w}\in{\bf V}$, the bilinear form
$\mathcal{B}(\cdot,\cdot)$ is bounded over ${\bf V}\times{\bf V}$. The
coercivity of $\mathcal{B}(\cdot,\cdot)$ on ${\bf V}$ follows from Korn’s
inequality,

$\|{\bf v}\|_{\mathcal{B}}^{2}:=\mathcal{B}({\bf v},{\bf v})\geq
2\mu\|{\boldsymbol{\varepsilon}}({\bf v})\|^{2}\geq c\mu\|{\bf v}\|_{{\bf
V}}^{2},\quad{\bf v}\in{\bf V}.$ (2.2)

Owing to these two properties, and since $|\ell({\bf v})|\leq\|{\bf f}\|_{{\bf
V}^{*}}\|{\bf v}\|_{{\bf V}}$, an application of the Lax–Milgram theorem
completes the proof of the next theorem.

######  Theorem 2.1.

For every $f\in{\bf V}^{*}$, the problem (2.1) has a unique solution ${\bf
u}\in{\bf V}$.

For the convergence analysis, we require the following regularity properties
[7, Lemma 2.2].

######  Theorem 2.2.

For every $f\in{\bf V}^{*}$, the weak solution ${\bf u}$ of (2.1) satisfies

$\|{\bf u}\|_{{\bf V}}+\lambda\|\nabla\cdot{\bf u}\|\leq C\,\|{\bf f}\|_{{\bf
V}^{*}}.$ (2.3)

Furthermore, if the domain $\Omega$ is convex and if ${\bf f}\in{\bf
L}^{2}(\Omega)$, then ${\bf u}\in{\bf V}\cap{\bf H}$ and $\|{\bf u}\|_{\bf
H}\leq C\,\|{\bf f}\|.$ In both estimates, the constant $C$ depends on
$\Omega$ and $\mu$, but is independent of $\lambda$.

######  Remark 2.1.

Our need for (2.3) is the only reason that Theorem 3.1 is restricted to two
dimensional problems.

For the finite element scheme, we introduce as usual a family of regular
quasi-uniform triangulations $\mathcal{T}_{h}$ of the domain $\Omega$ and let
$h=\max_{K\in\mathcal{T}_{h}}(h_{K})$, where $h_{K}$ denotes the diameter of
the element $K$. The space of polynomials with degree at most $1$ is denoted
by $P_{1}$, and the conforming finite element space is then defined by

${\bf V}_{h}=\\{\,{\bf w}_{h}\in{\bf V}:\text{${\bf w}_{h}|_{K}\in[P_{1}]^{2}$
for all $K\in\mathcal{T}_{h}$}\,\\}.$

In the standard conforming FEM we seek ${\bf u}_{h}\in{\bf V}_{h}$ such that
$\mathcal{B}({\bf u}_{h},{\bf v}_{h})=\ell({\bf v}_{h})$ for all ${\bf
v}_{h}\in{\bf V}_{h}$. To avoid locking when $\lambda$ is large, more
precisely, when $\lambda\gg 1/h$, we modify this scheme as follows: find ${\bf
u}_{h}\in{\bf V}_{h}$ such that

$\mathcal{B}^{\alpha}({\bf u}_{h},{\bf v}_{h})=\ell({\bf
v}_{h}),\quad\text{for all ${\bf v}_{h}\in{\bf V}_{h}$},$ (2.4)

where the bilinear operator $\mathcal{B}^{\alpha}$ is defined in the same way
as $\mathcal{B}$ in (2), except with $\lambda^{\alpha}$ in place of $\lambda$,
that is,

$\mathcal{B}^{\alpha}({\bf u},{\bf
v})=\int_{\Omega}[2\mu\,{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf v})+\lambda^{\alpha}(\nabla\cdot{\bf
u})(\nabla\cdot{\bf v})]\,\mathrm{d}{\boldsymbol{x}},$

for an appropriately chosen constant $\alpha$ that depends on $h$, $\lambda$
and $d_{\Omega}$, and that satisfies $0<\alpha\leq 1$; see Theorem 3.1 below.
The existence and uniqueness of the solution ${\bf u}_{h}$ of (2.4) is
guaranteed because the operator $\mathcal{B}^{\alpha}$ is bounded and
positive-definite on ${\bf V}_{h}$.

##  3 Convergence analysis

For convenience, we introduce the following notations: let

$\phi(\lambda)=\lambda^{\alpha},\quad\psi(\lambda)=(\lambda-\lambda^{\alpha})^{2}\lambda^{-\alpha}\quad{\rm
and}\quad\chi(\lambda)=(\lambda-\lambda^{\alpha})^{2}\lambda^{-\alpha-2}.$

For the error analysis, we use a projection operator $\Pi_{h}:{\bf V}\to{\bf
V}_{h}$ with the following error estimate properties:

$\|{\bf v}-\Pi_{h}{\bf v}\|+h\|\nabla({\bf v}-\Pi_{h}{\bf v})\|\leq
Ch^{2}\|{\bf v}\|_{\bf H}\,.$ (3.1)

The regularity results in Theorem 2.2 imply that

$\|{\bf u}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}\leq
C\|{\boldsymbol{\varepsilon}}({\bf u}-\Pi_{h}{\bf
u})\|^{2}+\|\sqrt{\lambda^{\alpha}}\,\nabla\cdot({\bf u}-\Pi_{h}{\bf
u})\|^{2}\leq C(1+\phi(\lambda))h^{2}\|{\bf u}\|_{\bf H}^{2}\leq
C\phi(\lambda)h^{2}\|{\bf f}\|^{2}.$ (3.2)

######  Theorem 3.1.

Assume that $\Omega\subseteq\mathbb{R}^{2}$ is convex with diameter
$d_{\Omega}$, and that ${\bf f}\in{\bf L}^{2}(\Omega)$. Let ${\bf u}$ and
${\bf u}_{h}$ be the solutions of problems (1.1) and (2.4), respectively, and
choose

$\alpha=\alpha_{*}(h,\lambda)\equiv\min\biggl{(}1,\frac{\log(d_{\Omega}/h)}{\log(\lambda)}\biggr{)}\quad\text{so
that}\quad\lambda^{\alpha}=\begin{cases}d_{\Omega}/h&\text{for
$\lambda>d_{\Omega}/h$,}\\\ \lambda&\text{otherwise,}\end{cases}$

with $0<\alpha\leq 1$. Then, there is a constant $C$ depending on $\Omega$ and
$\mu$, but independent of $\lambda$, such that

$\|{\bf u}_{h}-{\bf u}\|_{\mathcal{B}^{\alpha}}\leq Ch^{1/2}\|{\bf f}\|.$

Moreover, if $\mathcal{L}:{\bf L}^{2}(\Omega)\to\mathbb{R}$ is a bounded
linear functional (that is, $|\mathcal{L}({\bf w})|\leq\|\mathcal{L}\|\,\|{\bf
w}\|$), then

$|\mathcal{L}({\bf u}_{h})-\mathcal{L}({\bf u})|\leq C\,h\|{\bf
f}\|\|\mathcal{L}\|.$ (3.3)

###### Proof.

From the weak formulation in (2.1) and the numerical scheme in (2.4), we have
the property

$\displaystyle\mathcal{B}^{\alpha}({\bf u}_{h}-{\bf u},{\bf v}_{h})=$
$\displaystyle\mathcal{B}^{\alpha}({\bf u}_{h},{\bf
v}_{h})-\mathcal{B}^{\alpha}({\bf u},{\bf v}_{h})-\mathcal{B}({\bf u},{\bf
v}_{h})+\mathcal{B}({\bf u},{\bf v}_{h})$ $\displaystyle=$
$\displaystyle\mathcal{B}({\bf u},{\bf v}_{h})-\mathcal{B}^{\alpha}({\bf
u},{\bf v}_{h})$ $\displaystyle=$
$\displaystyle\int_{\Omega}(\lambda-\lambda^{\alpha})\nabla\cdot{\bf
u}\,\nabla\cdot{\bf v}_{h}\,\mathrm{d}{\boldsymbol{x}}\quad\text{for all ${\bf
v}_{h}\in{\bf V}_{h}$.}$ (3.4)

Using this formula, the fact that $\mathcal{B}^{\alpha}(\cdot,\cdot)$ is a
bounded bilinear operator, and applying the Cauchy–Schwarz inequality, we
obtain,

$\displaystyle\|{\bf u}_{h}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}$
$\displaystyle=\mathcal{B}^{\alpha}({\bf u}-\Pi_{h}{\bf u},{\bf
u}_{h}-\Pi_{h}{\bf
u})+\int_{\Omega}(\lambda-\lambda^{\alpha})\,\nabla\cdot{\bf
u}\,\nabla\cdot({\bf u}_{h}-\Pi_{h}{\bf u})\,\mathrm{d}{\boldsymbol{x}}$
$\displaystyle\leq\frac{1}{4}\|{\bf u}_{h}-\Pi_{h}{\bf
u}\|_{\mathcal{B}^{\alpha}}^{2}+\|{\bf u}-\Pi_{h}{\bf
u}\|_{\mathcal{B}^{\alpha}}^{2}+C\psi(\lambda)\|\nabla\cdot{\bf
u}\|^{2}+\frac{1}{2}\|\sqrt{\lambda^{\alpha}}\,\nabla\cdot({\bf
u}_{h}-\Pi_{h}{\bf u})\|^{2}\,.$

Using $\|\sqrt{\lambda^{\alpha}}\,\nabla\cdot({\bf u}_{h}-\Pi_{h}{\bf
u})\|^{2}\leq\|{\bf u}_{h}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2},$ then
cancelling the similar terms, we get

$\|{\bf u}_{h}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}\leq\|{\bf
u}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}+C\psi(\lambda)\|\nabla\cdot{\bf
u}\|^{2}.$

Hence, using (3.2) and the regularity results in Theorem 2.2, we have

$\|{\bf u}_{h}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}\leq
C\phi(\lambda)h^{2}\|{\bf f}\|^{2}+C\chi(\lambda)\|{\bf f}\|^{2}.$

For $\lambda>d_{\Omega}/h$, we choose $\alpha$ to balance the quantities
$\phi(\lambda)h^{2}$ and $\chi(\lambda)$. More precisely, we require that
$\lambda-\lambda^{\alpha}=\lambda^{1+\alpha}h/d_{\Omega}$. Dividing through by
$\lambda^{1+\alpha}$ yields $\lambda^{-\alpha}=1/\lambda+h/d_{\Omega}$, and
then taking the natural logarithm of both sides, we obtain
$\alpha=-\log(h/d_{\Omega}+1/\lambda)/\log(\lambda)\approx\log(d_{\Omega}/h)/\log(\lambda)=\alpha_{*}(h,\lambda)$.
If this condition is satisfied, then

$\|{\bf u}_{h}-\Pi_{h}{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}\leq Ch\|{\bf
f}\|^{2}.$ (3.5)

Using this together with (3.2), shows

$\|{\bf u}-{\bf u}_{h}\|_{\mathcal{B}^{\alpha}}^{2}\leq Ch\|{\bf f}\|^{2}.$
(3.6)

Now, if $\lambda\leq d_{\Omega}/h$, we simply repeat the above steps with
$\alpha=1$ (that is, standard conforming FEM) and arrive at the same bound.

To show the second estimate, we replace $\ell$ with $\mathcal{L}$ in (2.1),
and then, by Theorem 2.1, there exists a unique ${\bf u}_{\mathcal{L}}\in{\bf
V}$ satisfying

$\mathcal{B}({\bf u}_{\mathcal{L}},{\bf v})=\mathcal{L}({\bf v})\quad\text{for
all ${\bf v}\in{\bf V}$.}$ (3.7)

Furthermore, using Theorem 2.2 (with ${\bf u}_{\mathcal{L}}$ in place of ${\bf
u}$), we deduce that

$\|{\bf u}_{\mathcal{L}}\|_{{\bf V}}+\lambda\|\nabla\cdot{\bf
u}_{\mathcal{L}}\|\leq C\|\mathcal{L}\|_{{\bf V}^{*}}\quad{\rm and}\quad\|{\bf
u}_{\mathcal{L}}\|_{\bf H}\leq C\|\mathcal{L}\|.$ (3.8)

The modified conforming finite element approximation to ${\bf
u}_{\mathcal{L}}$ is the function ${\bf\Theta}_{h}\in{\bf V}_{h}$ such that

$\mathcal{B}^{\alpha}({\bf\Theta}_{h},{\bf v}_{h})=\mathcal{L}({\bf
v}_{h})\quad\text{for all ${\bf v}_{h}\in{\bf V}_{h}$.}$ (3.9)

By repeating the arguments that led to (3.6) and using (3.8), we obtain

$\|{\bf u}_{\mathcal{L}}-{\bf\Theta}_{h}\|_{\mathcal{B}^{\alpha}}^{2}\leq
Ch\|{\bf u}_{\mathcal{L}}\|_{\bf H}\|\mathcal{L}\|\leq Ch\|\mathcal{L}\|^{2}.$
(3.10)

Notice that

$\mathcal{B}({\bf w},{\bf v})-\mathcal{B}^{\alpha}({\bf w},{\bf
v})=\int_{\Omega}(\lambda-\lambda^{\alpha})(\nabla\cdot{\bf
w})(\nabla\cdot{\bf v})\,\mathrm{d}{\boldsymbol{x}}\quad\text{for all ${\bf
w}$, ${\bf v}\in{\bf V}$,}$

so by (3.7) and (3.9) it follows that

$\displaystyle\mathcal{B}^{\alpha}({\bf\Theta}_{h}-{\bf u}_{\mathcal{L}},{\bf
v}_{h})$ $\displaystyle=\mathcal{L}({\bf v}_{h})-\mathcal{B}({\bf
u}_{\mathcal{L}},{\bf v}_{h})+\mathcal{B}({\bf u}_{\mathcal{L}},{\bf
v}_{h})-\mathcal{B}^{\alpha}({\bf u}_{\mathcal{L}},{\bf v}_{h})$ (3.11)
$\displaystyle=\int_{\Omega}(\lambda-\lambda^{\alpha})(\nabla\cdot{\bf
u}_{\mathcal{L}})(\nabla\cdot{\bf
v}_{h})\,\mathrm{d}{\boldsymbol{x}}\quad\text{for all ${\bf v}_{h}\in{\bf
V}_{h}$.}$

By (2.4) and (3.9), we have $\mathcal{L}({\bf
u}_{h})=\mathcal{B}^{\alpha}({\bf\Theta}_{h},{\bf
u}_{h})=\mathcal{B}^{\alpha}({\bf
u}_{h},{\bf\Theta}_{h})=\ell({\bf\Theta}_{h})=\mathcal{B}({\bf
u},{\bf\Theta}_{h})$, whereas (3.7) and (2.1) imply that $\mathcal{L}({\bf
u})=\mathcal{B}({\bf u}_{\mathcal{L}},{\bf u})=\mathcal{B}({\bf u},{\bf
u}_{\mathcal{L}})$, so

$\mathcal{L}({\bf u}_{h})-\mathcal{L}({\bf u})=\mathcal{B}({\bf
u},{\bf\Theta}_{h}-{\bf u}_{\mathcal{L}})=\mathcal{B}^{\alpha}({\bf
u},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})+\int_{\Omega}(\lambda-\lambda^{\alpha})(\nabla\cdot{\bf
u})(\nabla\cdot({\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})\,\mathrm{d}{\boldsymbol{x}}.$

Since $\mathcal{B}^{\alpha}({\bf u}_{h},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})=\mathcal{B}^{\alpha}({\bf\Theta}_{h},{\bf
u}_{h})-\mathcal{B}^{\alpha}({\bf u}_{\mathcal{L}},{\bf
u}_{h})=\mathcal{L}({\bf u}_{h})-\mathcal{B}^{\alpha}({\bf
u}_{\mathcal{L}},{\bf u}_{h})=\mathcal{B}({\bf u}_{\mathcal{L}},{\bf
u}_{h})-\mathcal{B}^{\alpha}({\bf u}_{\mathcal{L}},{\bf u}_{h})$,

$\mathcal{B}^{\alpha}({\bf u},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})=\mathcal{B}^{\alpha}({\bf u}-{\bf
u}_{h},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})+\int_{\Omega}(\lambda-\lambda^{\alpha})(\nabla\cdot{\bf
u}_{\mathcal{L}})(\nabla\cdot{\bf u}_{h})\,\mathrm{d}{\boldsymbol{x}}$

and hence

$\mathcal{L}({\bf u}_{h})-\mathcal{L}({\bf u})=\mathcal{B}^{\alpha}({\bf
u}-{\bf u}_{h},{\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})+\int_{\Omega}(\lambda-\lambda^{\alpha})(\nabla\cdot{\bf
u})\nabla\cdot({\bf\Theta}_{h}-{\bf
u}_{\mathcal{L}})\,\mathrm{d}{\boldsymbol{x}}\\\
+\int_{\Omega}(\lambda-\lambda^{\alpha})\bigl{[}(\nabla\cdot{\bf
u}_{\mathcal{L}})\nabla\cdot({\bf u}_{h}-{\bf u})+(\nabla\cdot{\bf
u}_{\mathcal{L}})(\nabla\cdot{\bf u})\bigr{]}\,\mathrm{d}{\boldsymbol{x}}.$
(3.12)

By Theorem 2.2 and (3.8), $\|\nabla\cdot{\bf u}\|\leq C\lambda^{-1}\|{\bf
f}\|$ and $\|\nabla\cdot{\bf u}_{\mathcal{L}}\|\leq
C\lambda^{-1}\|\mathcal{L}\|$. Therefore, from the achieved estimates in (3.6)
and (3.10), and using $\|\nabla\cdot{\bf v}\|\leq C\lambda^{-\alpha/2}\|{\bf
v}\|_{\mathcal{B}^{\alpha}}$

$\displaystyle|\mathcal{L}({\bf u})-\mathcal{L}({\bf u}_{h})|$
$\displaystyle\leq\|{\bf u}-{\bf u}_{h}\|_{\mathcal{B}^{\alpha}}\|{\bf
u}_{\mathcal{L}}-\Theta_{h}\|_{\mathcal{B}^{\alpha}}+C\lambda\|\nabla\cdot{\bf
u}\|\|\nabla\cdot({\bf\Theta}_{h}-{\bf u}_{\mathcal{L}})\|$
$\displaystyle\qquad{}+C\lambda\|\nabla\cdot{\bf
u}_{\mathcal{L}}\|\|\nabla\cdot({\bf u}_{h}-{\bf
u})\|+C\lambda\|\nabla\cdot{\bf u}_{\mathcal{L}}\|\|\nabla\cdot{\bf u}\|$
$\displaystyle\leq
C\bigr{(}h+\lambda^{-\alpha/2}h^{1/2}+\lambda^{-1}\bigr{)}\|\mathcal{L}\|\|{\bf
f}\|.$

If $\lambda>d_{\Omega}/h$, then $0<\alpha\leq 1$ and
$\lambda^{\alpha}=d_{\Omega}/h$, so $\lambda^{-\alpha/2}\leq Ch^{1/2}$ and
$\lambda^{-1}\leq Ch^{1/\alpha}\leq Ch$, implying the estimate (3.3). If
$\lambda\leq d_{\Omega}/h$, then $\alpha=1$ so $\lambda-\lambda^{\alpha}=0$
and (3.3) again holds. ∎

######  Corollary 3.2.

There is a constant $C$, independent of $\lambda$, such that

$\|{\bf u}-{\bf u}_{h}\|\leq Ch\|{\bf f}\|\quad\text{and}\quad\|{\bf
u}_{h}-{\bf u}\|_{{\bf V}}\leq Ch^{1/2}\|{\bf f}\|.$

###### Proof.

By choosing $\mathcal{L}({\bf v})=\int_{\Omega}({\bf u}_{h}-{\bf u})\cdot{\bf
v}\,\mathrm{d}{\boldsymbol{x}}$ in (3.3), we have

$\|{\bf u}_{h}-{\bf u}\|^{2}=\mathcal{L}({\bf u}_{h}-{\bf u})=\mathcal{L}({\bf
u}_{h})-\mathcal{L}({\bf u})\leq Ch\|{\bf f}\|\|{\bf u}_{h}-{\bf u}\|,$

implying that $\|{\bf u}_{h}-{\bf u}\|\leq Ch\|{\bf f}\|$. By (2.2) and (3.6),

$\|{\bf u}_{h}-{\bf u}\|_{{\bf V}}^{2}\leq C\|\varepsilon({\bf u}_{h}-{\bf
u})\|^{2}\leq C\|{\bf u}_{h}-{\bf u}\|_{\mathcal{B}^{\alpha}}^{2}\leq Ch\|{\bf
f}\|^{2},$

implying the second estimate. ∎

######  Remark 3.1.

To discuss mixed boundary conditions, let $\boldsymbol{n}$ denote the outward
unit normal vector with respect to $\Omega$. Recall that $\Omega$ is a convex
polygonal domain, and let the closed line segment $\Gamma_{i}$ denote the
$i$th edge of the boundary $\Gamma=\partial\Omega$, ordered counterclockwise.
Thus, $\Gamma=\bigcup_{i=1}^{m}\Gamma_{i}$ for some finite $m$. We assume that
$u={\bf 0}$ on $\Gamma_{D}=\bigcup_{i=1}^{m_{0}}\Gamma_{i}$ and that
$\boldsymbol{\sigma}\boldsymbol{n}={\bf g}$ on
$\Gamma_{N}=\bigcup_{m_{0}+1}^{m}\Gamma_{i}$. In other words, homogeneous
Dirichlet boundary conditions apply on $\Gamma_{D}$, whereas traction boundary
conditions apply on $\Gamma_{N}$. The measure of $\Gamma_{D}$ must be
positive, and the Neumann data ${\bf g}\in{\bf L}^{2}(\Gamma_{2})$.

For $\ell({\bf v})=\int_{\Omega}({\bf f}\cdot{\bf
v})\,\mathrm{d}{\boldsymbol{x}}+\int_{\Gamma_{2}}({\bf g}\cdot{\bf
v})\,\mathrm{d}s$, the existence and uniqueness theorem (Theorem 2.1) remains
valid with $|\ell({\bf v})|\leq C(\|{\bf f}\|_{{\bf V}^{*}}+\|{\bf g}\|_{{\bf
L}^{2}(\Gamma_{2})})\|{\bf v}\|_{{\bf V}}$ and ${\bf V}=\\{\,{\bf
w}\in[H^{1}(\Omega)]^{2}:\text{${\bf w}=0$ on $\Gamma_{D}$}\,\\}$. The
conforming finite element scheme (2.4) is also well-defined if we put

${\bf V}_{h}=\\{\,{\bf w}_{h}\in{\bf V}:\text{${\bf w}_{h}|_{K}\in[P_{1}]^{2}$
for all $K\in\mathcal{T}_{h}$}\,\\}.$

The convergence results in Theorem 3.1 hold true provided the regularity
results in Theorem 2.2 remain valid. In these two theorems, we anticipate that
$\|{\bf f}\|_{{\bf V}^{*}}$ and $\|{\bf f}\|$ should be replaced with $\|{\bf
f}\|_{{\bf V}^{*}}+\|{\bf g}\|_{{\bf L}^{2}(\Gamma_{N})}$ and $\|{\bf
f}\|+\|{\bf g}\|_{[H^{1/2}(\Gamma_{N})]^{2}}$, respectively. For more details,
see for example Arnold and Falk [2], Brenner and Sung [7, Lemma 2.3], and
Brenner and Scott [6, Page 319].

######  Remark 3.2.

In this remark, we discuss the case of variable Lamé parameters $\lambda$ and
$\mu$. More precisely, when $\mu=\mu({\boldsymbol{x}})$ and
$\lambda=\Lambda\widehat{\lambda}({\boldsymbol{x}})$, where the constant
factor $\Lambda$ is expected to be very large (and approaches infinity for
materials which are increasingly incompressible). We assume that
$\mu_{\min}\leq\mu\leq\mu_{\max}$ and
$\widehat{\lambda}_{\min}\leq\widehat{\lambda}\leq\widehat{\lambda}_{\max}$ on
$\Omega$ for some positive constants $\mu_{\min}$, $\mu_{\max}$,
$\widehat{\lambda}_{\min},$ and $\widehat{\lambda}_{\max}.$ The regularity
results in Theorem 2.1 remains valid in this case with $\Lambda$ in place of
$\lambda$, where the constant $C$ depends on $\Omega$, $\mu$ and
$\widehat{\lambda}$ , but is independent of $\Lambda$; see Dick et al. [13,
Theorem 3.2].

We define the modified conforming FEM as in (2.4) with

$\mathcal{B}^{\alpha}({\bf u},{\bf
v})=\int_{\Omega}[2\mu\,{\boldsymbol{\varepsilon}}({\bf
u}):{\boldsymbol{\varepsilon}}({\bf
v})+\Lambda^{\alpha}\widehat{\lambda}(\nabla\cdot{\bf u})(\nabla\cdot{\bf
v})]\,\mathrm{d}{\boldsymbol{x}}.$

For the convergence results, we replace $\lambda$ with $\Lambda$ in the
statements of Theorem 3.1 and Corollary 3.2. In the proof of Theorem 3.1, we
replace $\lambda$ with $\Lambda$ in: (1) the functions $\phi,\psi,$ and
$\chi,$ (2) the text before (3.5) and after (3.6), and (3) the material after
(3.12).

##  4 Numerical experiments

In this section, we illustrate the theoretical convergence results. Based on
Corollary 3.2, we expect to observe $O(h)$ and $O(h^{1/2})$ convergence rates
in the ${\bf L}^{2}(\Omega)$-norm for ${\bf u}$ and $\nabla{\bf u}$,
respectively, uniformly in $\lambda$, provided that all the imposed
assumptions are satisfied. The Julia code used in these experiments is
available online [25].

### Example 1.

We chose $\Omega=(0,\pi)^{2}$ and chose the body force so that the exact
solution is

${\bf u}(\boldsymbol{x})=\begin{bmatrix}u_{1}(x_{1},x_{2})\\\
u_{2}(x_{1},x_{2})\end{bmatrix}=\begin{bmatrix}\bigl{(}\cos(2x_{1})-1\bigr{)}\sin(2x_{2})\\\
\bigl{(}1-\cos(2x_{2})\bigr{)}\sin(2x_{1})\end{bmatrix}+\frac{\sin(x_{1})\sin(x_{2})}{\lambda}\begin{bmatrix}1\\\
1\end{bmatrix},\quad{\boldsymbol{x}}=(x_{1},x_{2}).$

Notice ${\bf u}=\mathbf{0}$ on $\partial\Omega$. We used a family of
unstructured triangulations obtained by uniform refinement of an initial
coarse mesh. Table 1 compares the errors and convergence rates in the ${\bf
L^{2}}(\Omega)$-norm of the standard conforming, piecewise-linear Galerkin FEM
($\alpha=1$) with our conforming scheme ($\alpha=\alpha_{*}(h,\lambda)$ from
Theorem 3.1). The number of unknowns in the linear system for each mesh is
denoted by $N_{h}$. As $\lambda$ increases from $1{,}000$ to $10{,}000$ to
$100{,}000$ the material becomes more nearly incompressible and, as expected,
the standard method exhibits poor convergence. By contrast, our new method
achieves $O(h)$ convergence uniformly in $\lambda$.

Table 2 compares the errors in $\nabla{\bf u}$ for the two methods. Once
again, our new method is more accurate. In fact, it seems to achieve $O(h)$
convergence, implying that the second error bound in Corollary 3.2 is not
sharp.

Table 1: Errors and convergence rates for ${\bf u}$ in Example 1. $\lambda=1{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
---|---|---
$N_{h}$ | $h$ | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 2.96e+00 | 0.283 | 3.67e-01 | 1.209 | 0.423
2466 | 0.119 | 1.76e+00 | 0.748 | 1.68e-01 | 1.129 | 0.524
10114 | 0.060 | 7.02e-01 | 1.326 | 8.01e-02 | 1.067 | 0.624
40962 | 0.030 | 2.15e-01 | 1.706 | 3.90e-02 | 1.039 | 0.724
164866 | 0.015 | 5.86e-02 | 1.877 | 1.90e-02 | 1.034 | 0.825
$\lambda=10{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
$N_{h}$ | $h$ | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 3.73e+00 | 0.034 | 3.68e-01 | 1.208 | 0.317
2466 | 0.119 | 3.42e+00 | 0.127 | 1.68e-01 | 1.126 | 0.393
10114 | 0.060 | 2.57e+00 | 0.409 | 8.07e-02 | 1.061 | 0.468
40962 | 0.030 | 1.33e+00 | 0.949 | 3.96e-02 | 1.026 | 0.543
164866 | 0.015 | 4.77e-01 | 1.484 | 1.97e-02 | 1.012 | 0.618
$\lambda=100{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
$N_{h}$ | $h$ | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\|{\bf u}_{h}-{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 3.84e+00 | 0.003 | 3.68e-01 | 1.208 | 0.254
2466 | 0.119 | 3.80e+00 | 0.014 | 1.68e-01 | 1.125 | 0.314
10114 | 0.060 | 3.66e+00 | 0.054 | 8.08e-02 | 1.060 | 0.374
40962 | 0.030 | 3.20e+00 | 0.194 | 3.97e-02 | 1.025 | 0.435
164866 | 0.015 | 2.16e+00 | 0.569 | 1.97e-02 | 1.009 | 0.495
Table 2: Errors and convergence rates for $\nabla{\bf u}$ in Example 1. $\lambda=1{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
---|---|---
$N_{h}$ | $h$ | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 6.83e+00 | 0.281 | 1.31e+00 | 0.989 | 0.423
2466 | 0.119 | 4.12e+00 | 0.730 | 6.67e-01 | 0.972 | 0.524
10114 | 0.060 | 1.73e+00 | 1.253 | 3.44e-01 | 0.955 | 0.624
40962 | 0.030 | 5.96e-01 | 1.536 | 1.79e-01 | 0.945 | 0.724
164866 | 0.015 | 1.96e-01 | 1.605 | 9.29e-02 | 0.944 | 0.825
$\lambda=10{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
$N_{h}$ | $h$ | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 8.62e+00 | 0.034 | 1.31e+00 | 0.989 | 0.317
2466 | 0.119 | 7.89e+00 | 0.127 | 6.67e-01 | 0.972 | 0.393
10114 | 0.060 | 5.95e+00 | 0.406 | 3.45e-01 | 0.954 | 0.468
40962 | 0.030 | 3.14e+00 | 0.921 | 1.79e-01 | 0.944 | 0.543
164866 | 0.015 | 1.20e+00 | 1.389 | 9.32e-02 | 0.942 | 0.618
$\lambda=100{,}000$ | $\alpha=1$ | $\alpha=\alpha_{*}(h,\lambda)$
$N_{h}$ | $h$ | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\|\nabla{\bf u}_{h}-\nabla{\bf u}\|$ | Rate | $\alpha$
586 | 0.239 | 8.86e+00 | 0.003 | 1.31e+00 | 0.989 | 0.254
2466 | 0.119 | 8.77e+00 | 0.014 | 6.67e-01 | 0.972 | 0.314
10114 | 0.060 | 8.45e+00 | 0.054 | 3.45e-01 | 0.954 | 0.374
40962 | 0.030 | 7.39e+00 | 0.194 | 1.79e-01 | 0.944 | 0.435
164866 | 0.015 | 5.01e+00 | 0.561 | 9.32e-02 | 0.942 | 0.495
$44$$60$$48$$x_{1}$$x_{2}$$A$ Figure 1: The domain $\Omega$ for Cook’s
membrane problem. Figure 2: The deformed mesh obtained when each node is
displaced by ${\bf u}_{h}$. Left: using the standard conforming method
($\alpha=1$). Right: using our new method ($\alpha=\alpha_{*}(h,\lambda))$.
Figure 3: Convergence of the computed value of the vertical displacement
$u_{2}(A)$

### Example 2.

(Cook’s membrane problem [12]) This is a benchmark problem for linear
elasticity which combines bending and shearing. As illustrated in Figure 1,
$\Omega$ is the convex domain formed by connecting the vertices $(0,0)$,
$(48,44)$, $(48,60),$ and $(0,44)$. The left-hand side of $\Omega$ is clamped
(that is, ${\bf u}={\bf 0}$), and a constant shear load in the vertical
direction is applied on the right-hand side, indicated by the red arrow:
$\boldsymbol{\sigma}{\bf n}=(0,g)$ where $g$ is a positive constant. The
remaining part of the boundary is traction-free, that is,
$\boldsymbol{\sigma}{\bf n}={\bf 0}$, and the body force ${\bf f}={\bf 0}$.
The elastic material has Young’s modulus $E=1.12499998125$ MPa and Poisson
ratio $\nu=0.499999975$, so the Lamé parameters are

$\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)}=7.5\times
10^{6}\quad\text{and}\quad\mu=\frac{E}{2(1+\nu)}=0.375.$

The shear load gives rise to the functional $\ell({\bf
v})=\int_{44}^{60}gv_{2}\,\mathrm{d}x_{2}$, where ${\bf v}=(v_{1},v_{2})^{T}$,
${\boldsymbol{x}}=(x_{1},x_{2})$ and $g=1/16$. Starting from a common initial
mesh with $h=2.57881$, we solved for ${\bf u}_{h}$ using $\alpha=1$ and
$\alpha=\alpha_{*}(h,\lambda)=0.214$, and then displaced each node by the
respective solution ${\bf u}_{h}$ to obtain the deformed meshes shown in
Figure 2. Even with this relatively coarse mesh, the locking effect is
apparent.

We also compared the computed values for the vertical displacement at the
midpoint $A=(48,50)$ of the right-hand edge with the benchmark value
$u_{2}(A)=16.442$ [24, p. 3491]. Figure 3 shows that using $\alpha=1$ gives a
much less accurate result than using $\alpha=\alpha_{*}(h,\lambda)$.

## References

  * [1] M. Ainsworth and C. Parker, Unlocking the secrets of locking, Comput. Methods Appl. Mech. Engrg., 395, 115034 (2022).
  * [2] D. N. Arnold and R. S. Falk, Well-posedness of the fundamental boundary value problems for constrained anisotropic elastic materials, Arch. Rat. Mech. Anal., 98, 143–167 (1987).
  * [3] D. N. Arnold, G. Awanou and R. Winther, Nonconforming tetrahedral mixed finite elements for elasticity, Math. Models Methods Appl. Sci., 24, 783–796 (2014).
  * [4] I. Babuška and M. Suri, Locking effects in the finite element approximation of elasticity problems, Numer. Math., 62, 439–463 (1992).
  * [5] D. Boffi, F. Brezzi and M. Fortin, Mixed finite element methods and applications. Springer, Berlin (2013).
  * [6] S. C. Brenner and L. R. Scott, The mathematical Theory of Finite Element Methods, Texts in Applied Mathematics, third edition, Springer, New York, 2008.
  * [7] S. C. Brenner and L.-Y. Sung, Linear finite element methods for planar linear elasticity, Math. Comput., 59, 321–338 (1992).
  * [8] J. Bramwell, L. Demkowicz, J. Gopalakrishnan and W. Qiu, A locking-free hp DPG method for linear elasticity with symmetric stresses, Numer. Math., 122, 671–707 (2012).
  * [9] S. Chen, G. Ren and S. Mao, Second-order locking-free nonconforming elements for planar linear elasticity, J. Comput. Appl. Math., 233, 2534–2548 (2010).
  * [10] G. Chen and X. Xie, A robust weak Galerkin finite element method for linear elasticity with strong symmetric stresses, Comput. Methods Appl. Math., 16, 389–408 (2016).
  * [11] B. Cockburn, D. Schötzau and J. Wang, Discontinuous Galerkin methods for incompressible elastic materials, Comput. Methods Appl. Mech. Eng., 195, 3184–3204 (2006).
  * [12] R. D. Cook, Improved two-dimensional finite element, J. Structural. Division., 100, 1851–1863 (1974).
  * [13] J. Dick, Q. T. Le Gia, W. McLean, K. Mustapha, T. Tran, High-order QMC nonconforming FEMs for nearly incompressible planar stochastic elasticity equations, arXiv preprint arXiv:2402.11545
  * [14] D. A. Di Pietro and S. Nicaise, A locking-free discontinuous Galerkin method for linear elasticity in locally nearly incompressible heterogeneous media, Appl. Numer. Math., 63, 105–116 (2013).
  * [15] A. Edoardo, M. Stefano, L. Carlo and P. Luca, A dual hybrid virtual element method for plane elasticity problems, ESAIM: Math. Model. Numer. Anal., 54, 1725–1750 (2020).
  * [16] R. Falk, Nonconforming finite element Methods for the equations of linear elasticity, Math. Comput., 57, 529–550 (1991).
  * [17] Gopalakrishnan, J., Guzmán, J.: Symmetric nonconforming mixed finite elements for linear elasticity. SIAM J. Numer. Anal. 49, 1504–1520 (2011).
  * [18] P. Hansbo and M. G. Larson, Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’s method, Comput. Methods Appl. Mech. Eng., 191, 1895–1908 (2002).
  * [19] J. Hu and Z. C. Shi, Lower order rectangular nonconforming mixed finite elements for plane elasticity, SIAM J. Numer. Anal., 46, 88–102 (2008).
  * [20] F. Huo, R. Wang, Y. Wang and R. Zhang, A locking-free weak Galerkin finite element method for linear elasticity problems, Computers Math. Appl., 160, 181–190 (2024).
  * [21] J.G. Huang, S. Lin, Y. Yu, A novel locking-free virtual element method for linear elasticity problems, arXiv:2112.13848v3, 2022.
  * [22] F. Huo, R. Wang, Y. Wang and R. Zhang, An arbitrary order locking-free weak Galerkin method for linear elasticity problems based on a reconstruction operator, (2023).
  * [23] C.-O. Lee, J. Lee and D. Sheen, A locking-free nonconforming finite element method for planar linear elasticity, Adv. Comput. Math., 19, 277–291 (2003).
  * [24] Y. Liu and J. Wang, A locking-free $P_{0}$ finite element method for linear elasticity equations on polytopal partitions, IMA J. Numer. Anal., 42, 3464–3498, (2022).
  * [25] W. McLean, https://github.com/billmclean/ControlledLocking.jl.git, (2024).
  * [26] S. Mao and S. Chen, A quadrilateral nonconforming finite element for linear elasticity problem, Adv. Comput. Math., 28, 81–100 (2008).
  * [27] B. D. Reddy and D. van Huyssteen, A virtual element method for transversely isotropic elasticity, Comput. Mech., 64, 971–988 (2019).
  * [28] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, RAIRO Math. Modeling Num. Anal. 19, 111–143 (1985).
  * [29] S.-C. Soon, B. Cockburn and H. K. Stolarski, A hybridizable discontinuous Galerkin method for linear elasticity, Int. J. Numer. Meth. Eng., 80, 1058–1092 (2009).
  * [30] L.B. Da Veiga, F. Brezzi and L. D. Marini, Virtual elements for linear elasticity problems, SIAM J. Numer. Anal., 51, 794–812 (2013).
  * [31] M. Vogelius, An analysis of the p-version of the finite element method for nearly incompressible materials, Numer. Math., 41, 39–53 (1983).
  * [32] C. Wang, J. Wang, R. Wang and R. Zhang, A locking-free weak Galerkin finite element method for elasticity problems in the primal formulation, J. Comput. Appl. Math., 307, 346–366 (2016).
  * [33] T. P. Wihler, Locking-free adaptive discontinuous Galerkin FEM for linear elasticity problems, Math. Comp., 75, 1087–1102 (2006).
  * [34] Y. Yang and S. Chen, A locking-free nonconforming triangular element for planar elasticity with pure traction boundary condition, J. Comput. Appl. Math., 233, 2703–2710 (2010).
  * [35] B. Zhang, J. Zhao, S. Chen and Y. Yang, A locking-free stabilized mixed finite element method for linear elasticity: the high order case, CALCOLO, 55, 1–17 (2018).
  * [36] B. Zhang, J. Zhao, Y. Yang and S. Chen, The nonconforming virtual element method for elasticity problems, J. Comput. Phys., 378, 394–410 (2019).
  * [37] J. Zhao, T. Wang and B. Zhang, The stabilized nonconforming virtual element method for linear elasticity problem, J. Sci. Comput., 92, 68, (2022).