Finite Element Simulation of Thermo-Mechanical Model with Phase Change

: In this work, we consider a mathematical model and ﬁnite element implementation of heat transfer and mechanics of soils with phase change. We present the construction of the simpliﬁed mathematical model based on the deﬁnition of water and ice fraction volumes as functions of temperature. In the presented mathematical model, the soil deformations occur due to the porosity growth followed by the difference between ice and water density. We consider a ﬁnite element discretization of the presented thermoelastic model with implicit time approximation. Implementation of the presented basic mathematical model is performed using FEniCS ﬁnite element library and openly available to download. The results of the numerical investigation are presented for the two-dimensional and three-dimensional model problems for two test cases in three different geometries. We consider algorithms with linearization from the previous time layer (one Picard iteration) and the Picard iterative method. Computational time is presented with the total number of nonlinear iterations. A numerical investigation with results of the convergence of the nonlinear iteration is presented for different time step sizes, where we calculate relative errors for temperature and displacements between current solution and reference solution with the largest number of the time layers. Numerical results illustrate the inﬂuence of the porosity change due to the phase-change of pore water into ice on the deformation of the soils. We observed a good numerical convergence of the presented implementation with the small number of nonlinear iterations, that depends on time step size. previous time layer (one Picard iteration) and the Picard iterative method. We present computational time with the total number of nonlinear iterations for the two-dimensional and three-dimensional model problems for two test cases in three different geometries. Figures of the relative errors in L 2 norm for temperature and displacements between the given solution and the reference solution using a smaller time step size illustrate the convergence of the presented ﬁnite element model for different numbers of different time step sizes and linearization algorithms. Displacements and temperature ﬁelds with phase-change interface illustrate an inﬂuence of the porosity change due to phase-change of pore water into ice on the deformation of the soils. From the presented numerical results for the basic thermo-elasticity model with phase-change of porous water to ice, we observe good convergence with a few numbers of the nonlinear iterations. We saw that inﬂuence of the time step size is larger for the displacements ﬁeld than for the temperature ﬁeld. Furthermore, we observe that more nonlinear iterations are needed for a larger time step size that converges to the one iteration per time layer. The algorithms with one Picard iteration provide worse results than the algorithm with nonlinear iterations, where errors are larger for less number of time layers. We observed a good numerical convergence of the presented implementation with the small number of nonlinear iterations that depends on time step size for two- and three-dimensional model problems.


Introduction
Permafrost covers almost a quarter of the land area in the Northern Hemisphere. Many people around the world live in permafrost and seasonal frozen areas in Alaska, Canada, and Russia (Na and Sun [1], Yu et al. [2], Tounsi et al. [3], Marchenko et al. [4], Knoblauch et al. [5]). The ground freezing and thawing process can change the shape of the land surface and damage constructions and buildings in the permafrost area (Zhou and Meschke [6], Sweidan et al. [7], Xu et al. [8]). Moreover, buildings constructed in the permafrost area can create heat that thaws the underlying grounds. Thawed soil becomes soft and uneven and can destroy houses and buildings. In order to prevent ground warming, buildings are constructed on pile foundations (Shang et al. [9], Nixon [10], Foriero and Ladanyi [11]). Construction of roads, bridges, railways, and other types of transport infrastructure is also a major problem in the permafrost zone as soil thawing damages them and requires constant repairs to keep them safe. Moreover, the overall rise in annual temperatures and releasing greenhouse gases into the atmosphere lead to further climate change on a global level. The global climate is getting warmer and making the frozen ground thaw (Nixon [10], Buteau et al. [12], Delisle [13]). Permafrost contains deposits of carbon matter, and global climate warming can release carbon dioxide and methane the freezing and thawing of soils. The mathematical model describes a frost heaving using the porosity rate function, which simulates the growth of ice lenses as average growth in porosity. The model was implemented in the commercial finite-element system ABAQUS.
Simulations of the applied problems require the construction of complex domains. In order to numerically solve a mathematical model, unstructured computational meshes should be constructed to the accurate approximation of the problem in geometrically complex domains. The common approach used for approximation of multiphysical problem on an unstructured grid is the finite element method (Brenner and Scott [31], Hughes [32]). In our previous works (Pavlova et al. [33], Vabishchevich et al. [34]), we considered the heat transfer problem with phase change and presented finite element approximation on unstructured grids for complex geometries (Pavlova et al. [33], Vabishchevich et al. [34]). The results of simulations for several applied problems in three-dimensional formulations illustrated the applicability of the method. For numerical implementation, we used an open-source finite-element library FEniCS (Logg et al. [35]). In FEniCS, we used a PETSc based solver to solve a large system of equations on high-performance computing systems (Balay et al. [36]). The PETSc library is a leading software framework for parallel scientific computing. In (Kolesov et al. [37]), we presented an implementation of the poroelasticity model using the PETSc library, where splitting schemes were constructed and investigated for multidimensional test problems on high-performance systems with distributed memory. Recently, we considered the construction of the mixed dimensional mathematical model for artificial freezing problems in permafrost area in (Vasilyeva et al. [38]), where we presented a multiscale method with additional multiscale basis functions for the solution of the Stefan problem in heterogeneous media for two and three-dimensional formulations. In this work, we continue the development of mathematical models with finite element approximation for problems in permafrost.
In this work, we consider a basic thermoelastic model that takes into account the deformations of porous media due to the phase change of pore water to ice. We start with the description of a thermo-mechanical model of soils with phase-change. The mathematical model is based on a definition of water fraction volume as a function of temperature using the relationship for unfrozen water content. In (Tice et al. [39], Michalowski [40]), the estimation of unfrozen water content function based on experimental tests is given, where parameters of the model depend on soil type and, in general, are different for freezing and thawing processes. Using mass conservation of frozen and thawed water, we describe ice volume fraction. As the sum of ice and unfrozen water fraction equals porosity, we can define a relationship for porosity as a function of temperature. Soil deformations occur due to porosity growth caused by the difference in ice and water density. The presented mathematical model is described by a nonlinear system of equations for temperature and displacements. The novelty of the paper is related to (1) derivation of the simplified mathematical model that takes into account porous water/ice phase-change and its influence on the mechanical deformations of soils; (2) construction and implementation of finite element approximation in order to numerically investigate model with an implicit scheme for approximation by time and Picard iterations for nonlinear coefficients for two and three-dimensional formulations. Implementation of the presented basic mathematical model is performed using FEniCS finite element software (Logg et al. [35]). The code is written using C++ Programming Language to perform fast calculations of multidimensional problems. We present numerical results for model problems in two-and three-dimensional formulations for two test cases in three different geometries. A numerical investigation is presented for different time step sizes, where we calculate relative errors for temperature and displacements. We consider algorithms with linearization from the previous time layer (one Picard iteration) and the Picard iterative method. Numerical results are presented with the number of nonlinear iterations. Finally, we discuss the computational time of the presented implementation with linearization from the previous time layer and the Picard iterative method. For the numerical solution of the arisen linear system of equations, we use a direct solver for both temperature and displacements fields. In future works, we will consider the numerical solution of the complex thermoelasticity model with phase-change on high-performance computing systems. The more complex models will be considered in future works by adding elastoplastic effects, the effect of the fluid flow, and ice lens formations.
The paper is organized as follows. In Section 2, we describe the basic concepts used in the construction of the mathematical model. In Section 3, we describe heat transfer processes in porous media and consider the mechanics of soils due to phase-change of the pore water into ice. For the numerical solution of the thermo-mechanical model, we present a finite element approximation in Section 4. Numerical results for model problems are presented in Section 5 for two-and three-dimensional formulations. The paper ends with conclusions.

Preliminaries
Let V w , V i , and V s be water, ice, and solid volumes of porous media, and where V is a total volume, and V v is a volume of void space that is filled with water or/and ice. The volume representation is shown in Figure 1 (Zhang [30]).
Let Θ w , Θ i , and Θ s be water, ice, and solid fraction volumes where superscripts s, w, and i denote solid, unfrozen water and ice, respectively. The porosity is defined as follows We define unfrozen water content as follows where ρ w and ρ s are density of water and solid, respectively. Note that the unfrozen water content function parameters depend on soil type and, in general, are different for freezing and thawing processes. The following model is used in (Tice et al. [39], Michalowski [40]) where T is a given temperature, T f is a freezing temperature and a is a model parameter [ • C −1 ], w and w * are maximum and minimum water contents. Therefore, using w, we have In order to find a relationship that determines porosity changes due to water freezing or thawing processes, mass conservation of water in frozen and thawed phases is used where ρ i is an ice density and V w is a water volume for T ≥ T f . Therefore Next, we can find a relationship for porosity An obtained relationship for porosity, water, and ice fraction volumes will be used in the thermo-mechanical model.

Mathematical Model
In this section, we present a thermo-mechanical mathematical model for the simulation of heat transfer and mechanics in soils. We start with the equation for heat transfer in porous media. Then, we describe the mechanical model. Finally, we present the thermomechanical problem formulation.

Heat Transfer in Porous Media
The heat transfer in porous media is based on the energy conservation principle that can be written as follows (Michalowski and Zhu [28], Michalowski [40]) where T is a temperature, L is the latent heat of fusion of water per unit mass, k is a heat conductivity, and c is a specific heat capacity, and ρ is a density. By mixture theory for saturated soils, we have where c j and ρ j are heat capacities (per unit mass) and densities for water, ice, and solid (j = w, i, s).
Heat conductivity can be calculated using the following model (Michalowski and Zhu [28]) where k w , k i and k s are water, ice and solid heat conductivities. For water, ice, and solid fraction volumes and porosity, we have the following relationships for the given relationship for water content with w * = 0 We rewrite the heat transfer Equation (2) in the following form for the water content relationship given in Equation (9).
We note that more complex models can be used to simulate phase change processes, for example, taking into account the ice lenses formation process (Michalowski and Zhu [28], Zhang [30]). In this work, we focus on basic ingredients for porosity change that relate to the ratio of water density and ice density (ρ w /ρ i ≈ 1.09).

Mechanics of Soils Due to Phase Change
Let T t = T(t) and T t+dt = T(t + dt) be the temperature at times t and t + dt, where dt is the time step. We define porosity at time t and t + dt as follows (see Figure 2) where dV is a change of volume. The volumetric strain due to porosity change is given as follows (Zhang [30]) where φ t+dt = φ is a current porosity and dφ = (φ t+dt − φ t ) is porosity change.
For the linear elasticity model, we have div(dσ) = 0, dσ = C : dε, where C is a stiffness tensor, σ and ε are total stress and strain tensors, respectively. The total strain is induced by loading and by a thermal process where dε e is an elastic mechanical increment, dε f is a thermal porosity growth increment, I is the identity matrix and u is displacements. For the two-dimensional isotropic media case, we have where u = (u 1 , u 2 ), λ and µ are Lame parameters, and β is a thermal expansion due to porosity change.

Thermo-Mechanical Model
Finally, we obtain the following thermo-mechanical model in the computational domain Ω The system of Equations (11a) and (11b) is an uncoupled system of equations because the presented model neglects the influence of the displacement field on the temperature distribution. In general, the mathematical model can be constructed in a coupled way with an additional part that takes into account the influence of the deformation on the temperature field. In this work, we investigate and implement a simplified model, where the temperature does not depend on displacements and can be calculated first. After using the calculated temperature field, we can find the value of porosity change dφ to calculate a correspondent displacement field of soils.
We consider the thermo-mechanical model with the following initial conditions and boundary conditions where n is the outward unit normal vector to ∂Ω and Γ 1 ∪ Γ 2 = Γ 3 ∪ Γ 4 = ∂Ω. Note that in (12), (13a) and (13b) the general formulation of the initial and boundary conditions is presented. In Figure 3, we present one possible variant of the domain Ω and the boundaries for which we can set such conditions. In this variant, we set Γ 2 ⊂ Γ top (green color), and Γ 1 = ∂Ω \ Γ 2 (purple color), and Γ 3 = Γ top (blue color), and Γ 4 = ∂Ω \ Γ 3 (red color).
The presented basic mathematical model (11a) and (11b), is the simplification of the model presented in (Zhang [30]). We use the presented model to illustrate the influence of the water/ice phase-change in porous media on deformations of the soils. Next, we will present a finite element approximation of the thermoelastic model with porous water/ice phase-change. In the numerical implementation, we don't use a proprietary "black box" software and present an openly available implementation based on the FEniCS library using a C++ programming language for fast simulations of multidimensional problems.

Finite Element Approximation
For the numerical solution of the mathematical model, we use a finite element method with implicit time approximation (Brenner and Scott [31], Kolesov et al. [37], Vasilyeva et al. [38]). We start with variational formulation of the problem (11a), (11b), (13a) and (13b) with initial conditions (12). Then, the formulation of the discrete problem is presented in matrix form.

Variational Formulation
and Q = H 1 (Ω) be the functional spaces for displacements and temperature, T ∈ Q and Multiplying (11a) and (11b) by test functions r ∈ Q and v ∈ V, respectively, and integrating by parts with boundary conditions given by (13a) and (13b), we obtain the following variational formulation of the thermo -mechanical problem (Brenner and Scott [31]): find T ∈ Q and u ∈ V such that where operator: is the inner product between tensors. For approximation by time, the Backward Euler scheme is used ∂T ∂t where τ is the size of the time step, and m is a time step number, m = 0, 1, ..., N t . Therefore, the variational formulation is the following: find T m+1 ∈ Q and u m+1 ∈ V such that Because, in our model, the heat equation does not depend on deformations, we can split the heat equation from the mechanical model.
To solve the nonlinear problem for temperature, the Picard iterative process is used (Logg et al. [35]). At each time layer, we solve the following linear system of equations: where l is the number of nonlinear iteration (l = 0, N nl ) and initial guess T m+1,0 = T m . In order to stop iterations, we calculate the relative difference in L 2 norm and use the following criteria where relative tolerance is given in % and ||v|| L 2 = Ω v 2 dx. We get the following algorithm: • Find T m+1,l+1 ∈ Q such that: where with convergence criteria (15) and where and porosity is the function of temperature Here σ(u) = 2 µ ε + λ tr (ε)I and ε = (∇u + ∇u T )/2.

Finite Element Discretization
Let T h be a grid partition of the computational domain Ω into finite elements K i where N c is a number of grid cells. Relative to T h , we define finite-dimensional spaces V h ⊂ V and Q h ⊂ Q. We suppose that Q h and V h are composed of piecewise linear basis functions on grid cells (Brenner and Scott [31], Logg et al. [35]) where ψ j and ϕ j are linear basis functions defined on T h , Thus, we obtain the following discrete systems in matrix form: In the presented algorithm, we first solve a temperature equation using the Picard iterative method. Then, we find displacements for the given temperature change.
For the implementation of the presented basic mathematical model with finite element approximation, we use the FEniCS finite element software (Logg et al. [35]) (version 2019.1.0+). The FEniCS library is a tool for automated scientific computing, with a particular focus on the solution of differential equations by the finite element method. The code for the presented thermoelastic model is written using C++ Programming Language to perform fast calculations of multidimensional problems. We use the GMSH program to construct computational geometries, and unstructured grids (Geuzaine and Remacle [41]). Paraview program is used for visualization of the results (Ahrens et al. [42]). The code and geometry files are openly available on Bitbucket at https://bitbucket.org/vmasha/thermoelastic/.

Numerical Results
In this section, we present the results of numerical simulations for two-and threedimensional model problems. In particular, we numerically examine the presented ap-proximation for different time step sizes and their influences on the number of nonlinear iterations. We present errors for temperature and displacements between the given solution and reference solution using a smaller number of time steps. Finally, we discuss the computational time of the presented implementation.
Phase-change parameters For the unfrozen water content model (1), we set Illustration of water fraction volume Θ w , ice fraction volume Θ i and porosity φ is presented in Figure 4.

Two-Dimensional Problem
In this section, we present numerical results for the two-dimensional domain Ω. Three geometries of the computational domain are considered: • Geometry 1 (see Figure 5a).
with L x = 6 m and L y = 3 m. On the boundary Γ H (green line), we set temperature T 1 = T H and T 1 = T C on the top boundary Geometry 2 (see Figure 5b).

Rectangular domain with two perforations (pipes)
with L x = 6 m and L y = 3 m and Ω p is the two circle domains. Two subdomains have different properties (Ω = Ω 1 ∪ Ω 2 ), where Ω 1 is the subdomain of soil (blue color in Figure 5b) and Ω 2 is the subdomain of pipe walls (red color). On the boundary of pipes Γ H = ∂Ω ∩ ∂Ω 2 (green line), we set temperature T 1 = T H and T 1 = T C on the top boundary Γ C = Γ top .
Computational domain with two subdomains Ω = Ω 1 ∪ Ω 2 , where Ω 1 is the subdomain of soil (blue color in Figure 5c), and Ω 2 is the subdomain of the pillar (red color). On the green boundary Γ H , we set temperature T 1 = T H and T 1 = T C on other top boundary In second subdomain Ω 2 , we simulate stiff material E 2 = 10 11 Pa, ν 2 = 0.3 with low heat conductivity Two-dimensional domains with unstructured grids are represented in Figure 5. The first geometry (Figure 5a) shows a thawing process caused by positive temperature on the boundary. This example illustrates the impact of a building/construction on the soil when poor floor insulation causes thawing during the winter period. The second geometry in Figure 5b illustrates the impact of underground pipes filled with warm water. This example shows how warm temperature can change the temperature field around and lead to thawing if poorly insulated. The third geometry (Figure 5c) illustrates deformation and heat transfer under the building/construction, where subdomain Ω 2 is the construction pillar. The example represents temperature changes under the construction and affects the pillar.
Unstructured meshes are constructed with local refinement near boundaries Γ C and Γ H . Computational mesh sizes and the number of degrees of freedom are the following: We simulate for t max = 10 6 s with two initial conditions for temperature: • Test 1. Frozen initial condition with T 0 = −2 • C. • Test 2. Unfrozen initial condition with T 0 = 2 • C.
As boundary conditions, we set Robin boundary conditions on Γ C and Γ H On other boundaries ∂Ω/(Γ C ∪ Γ H ), we set zero Neumann boundary conditions. For displacements, we set σ 1 = 0 Pa u 2 = 0 m on the bottom boundary, and u 1 = 0 m, σ 2 = 0 Pa on the left and right boundary. As initial condition, we set u 0 = 0 m. Numerical results for Test 1 are presented in Figure 6 and for Test 2 in Figure 7. Results have been performed with 300 time steps using Picard iterations with = 1%. On the left figure, we presented a temperature field distribution at the final time. The corresponding porosity field is shown in the middle figure. On the right figure, we depict a magnitude of displacements field on the moved mesh, where displacements are exaggerated by a factor of 25. The temperature field and porosity are depicted with the phase-change interface (white color line).
In Figure 6, we depicted results for Test 1 with frozen initial condition for temperature, T 0 = −2 • C. On the left figures, we observe the thawing process near the boundary Γ H , where we set Robin boundary conditions with temperature, T 1 = 10 • C. On the right figures, we see nearly vertical displacements caused by soil thawing. The results in Figure 6 illustrates the impact of a building/construction or pipes on the soil when poor floor insulation causes thawing during the winter period. We observe how warm temperature can change the temperature field around and lead to thawing if poorly insulated. The thawing has changed the porosity field as pore water has a larger density, and it occupies a smaller volume and leads to soil deformations.
The results of the numerical simulations for Test 2 are presented in Figure 7. In this example, we consider a case with unfrozen initial condition for temperature, T 0 = 2 • C. By setting freezing boundary conditions on the top boundary Γ H , we simulate the influence of the different constructions on the temperature and displacements fields. From this example, we observe how different boundary conditions on Γ C and Γ H affect the phasechange interface. From the middle figures of the porosity field related to the temperature at the final time, we observe porosity change due to freezing of the pore water. The right figures show displacements that occur for the given porosity change.
Next, we numerically investigate the solution of the problem using Picard iterations with a different number of time steps N t = 50, 100, 200. To compare the difference between solutions with different time steps, we calculate relative L 2 errors in % where T re f and u re f are reference solutions. Because we cannot calculate the analytical solution for the considered geometries in two and three-dimensional formulations, as the reference solutions in order to investigate numerical convergence of the presented algorithm with a different number of time steps and calculate errors, we use numerical results for N t = 300 with Picard iterations. In Figures 8 and 9, we present relative errors in % for 10 time layers using Picard iterations for Test 1 and Test 2, respectively. In the first and second columns, relative errors for temperature and displacements are depicted. In the third column, the number of Picard iterations are presented for each time step. We observe the influence of the time step size on the solutions and have the following conclusions from the presented results: • From the first and second columns in Figures 8 and 9, we observe larger errors for Test 2 than for Test 1 for all computational domains at initial time layers. For example, in Geometry 1 using N t = 50 time iterations, we have (6.3%, 3 . We obtain similar behavior for Geometry 2 and 3 and see that the Picard iterative method needs more iterations to converge for a small number of time steps N t . In Figure 10, we present the difference between solution using linearization from the previous time layer (one Picard iteration) and solution using Picard iterations using the same numbers for time steps N t = 50, 100, 200, and 300. We obtain the following conclusions: • For Test 1, we have a small difference between solution using linearization from the previous time layer and solution using Picard iterations. For example, in Test 1 (Geometry 1), we have near or less than one percent of errors for temperature and displacements for the algorithm with linearization from the previous time layer for N t = 50, 100, 200, 300. This happens because, for this test case, phase change interface movement occurs in a smaller domain and leads to a less nonlinear process for the temperature distribution with a small number of nonlinear iterations. In Test 2, we obtain a more nonlinear process with a larger thawing domain. This leads to a larger number of nonlinear iterations at initial time layers to converge, and therefore larger errors occur between solution using linear algorithm and algorithm with Picard iterations. From the presented results, we see how the computational domain and initial conditions affect the solution for a different number of time steps. We also observe the influence on the solution of the linearization method.
In Table 1, we present the computational time with the sum of Picard iterations over time (∑ iter). Note that ∑ iter = N t for linearization from the previous time layer. We have the following conclusions: • For the Picard iteration method, we observe that the sum of iterations is larger for Test 2 than for Test 1 for all geometries. For Geometry 1, we have in total 72 nonlinear iterations in Test 1 and 105 nonlinear iterations in Test 2 for N t = 50 time iterations.
We have more nonlinear iterations for Test 2, because, as we see before, the temperature field for initial and boundary conditions of Test 2 we obtain faster phase change interface moving and therefore we have more nonlinear iterations to converge. Solution time also depends on the size of the system that is solved at each time/nonlinear iteration. The size of the systems was presented above with the mesh size for all geometries. Solution time for Geometry 1 is slightly smaller than for Geometry 2 and 3 because we have a smaller number of degrees of freedom for temperature and displacements.
Solution time is presented in seconds and does not include the time of saving solution, where we have approximately 15 s for saving solution ten times. Computations were performed on 2.3 GHz Intel Core i7 with 32 GB of memory. For a solution, we used a direct solver for both temperature and displacements fields.

Three-Dimensional Problem
Finally, we present the numerical results for a three-dimensional domain based on the cuboid Ω m = [0, L x ] × [0, L y ] × [0, L z ] with L x = L y = 6 m and L z = 3 m.
Similarly to previous results, we consider three geometries: • Geometry 1 (see Figure 11a).
We set Ω = Ω m is the soil domain. On the boundary Γ H (elliptical surface on the top boundary), we set temperature T 1 = T H and T 1 = T C on the top boundary Γ C = Γ top /Γ H .
Computational domain with two cylindrical perforations (pipes), Ω = Ω m ∪ Ω p . We have two subdomains with different properties Ω = Ω 1 ∪ Ω 2 , where Ω 1 is the subdomain of soil (blue color), and Ω 2 is the thin subdomain of pipe walls (red color). On the boundary of pipes Γ H = ∂Ω ∩ ∂Ω 2 , we set temperature T 1 = T H and T 1 = T C on the top boundary Γ C = Γ top .
Computational domain with two subdomains Ω = Ω 1 ∪ Ω 2 , where Ω 1 is the subdomain of soil (blue color), and Ω 2 is the subdomain of the pillar (red color). On the boundary Γ H of the construction, we set temperature T 1 = T H and T 1 = T C on other top boundary In the second subdomain Ω 2 , we simulate material with similar parameters as in the previous two-dimensional case.
Numerical calculations were performed with 300 time steps using Picard iterations with = 1%. On the Figures 12 and 13, we present results at the final time for Test 1 and 2, respectively. On the left figures, we presented a temperature field distribution. On the middle figures, we depicted a phase-change interface with domain clip for x 2 = L y /2 m. The magnitude of the displacement field on the moved mesh is presented on the right, where displacements are exaggerated by a factor of 25.
From the presented results for Geometries 1, 2, and 3 in Test 1, we observed similar behavior as in the previous two-dimensional examples. Results for Geometry 1 on Figure 12a illustrate the influence of the construction with temperature T H = 10 • C and coefficient γ H = 20 W · m −2 · • C −1 without good insulation of the floor. We saw how the thawed zone corresponded to the soil deformations under construction. Results in Figure 12b for Geometry 2 are the present influence of the pipes laying under the ground, and we observe the corresponding phase-change interface and deformations on the thawed zone around pipes. The geometry 3 results in Figure 12c illustrate temperature, phase-change interface, and deformation under construction (warehouse structure) and near the cylindrical pillar.
In Figure 13, the results of the numerical simulations for Test 2 are presented, where an unfrozen initial condition is given for temperature, T 0 = 2 • C. In results for Geometries 1, 2, and 3, we observe the freezing from the top of the boundary with T C = −15 • C. In the presented geometries, we have different boundary Γ H with T H = 10 • C. We observe how a temperature field with a different shape of phase-change interface leads to the different deformations.
We present relative L 2 errors between solutions using Picard iterations with N t = 300 and solutions with different number of time steps N t = 50, 100, 200. In Figures 14 and 15, we present relative errors in % for 10 time layers using Picard iterations for Test 1 and Test 2, respectively. In the first and second columns, relative errors for temperature and displacements are depicted. We observe larger errors for temperature in Test 2 than for Test 1 for all computational domains. Similar to the results for the two-dimensional problem, the influence of the time step size is larger for displacements than for temperature for both tests. Temperature and displacements errors are decreasing over time. In the third column, the numbers of Picard iterations are presented for each time step. We observe that the Picard iterative method needs more iterations to converge for a small number of time steps N t .
Finally, we present the computational time with the sum of Picard iterations over time in Table 2. For linearization from the previous time layer, we have ∑ iter = N t . We observe that the sum of nonlinear iterations is slightly larger for Test 2 than for Test 1. (c). Geometry 3 (3D) Figure 11. Computational domains and grids for three-dimensional problems.
(a). Geometry 1 (3D) (b). Geometry 2 (3D) (c). Geometry 3 (3D) Figure 12. Test 1 for a three-dimensional problem. Temperature, temperature clip x 2 = L y /2 m with the phase change interface, and displacement fields at the final time (from left to right). (c). Geometry 3 (3D) Figure 13. Test 2 for a three-dimensional problem. Temperature, temperature clip x 2 = L y /2 m with the phase change interface, and displacement fields at the final time (from left to right).

Conclusions
We have presented the mathematical model and finite element implementation for heat transfer and the mechanics of soils with phase change. The constructed simplified mathematical model is based on a definition of water and ice fraction volumes that are functions of temperature. The soil deformations occur due to porosity growth that happens due to ice and water density differences. The finite element approximation for the thermomechanical model is presented with implicit approximation by time and Picard iterations. We present an openly available implementation of the presented basic mathematical model using the FEniCS finite element library. We present the numerical investigation of the discrete model for the different time step sizes for algorithms with linearization from the previous time layer (one Picard iteration) and the Picard iterative method. We present computational time with the total number of nonlinear iterations for the two-dimensional and three-dimensional model problems for two test cases in three different geometries. Figures of the relative errors in L 2 norm for temperature and displacements between the given solution and the reference solution using a smaller time step size illustrate the convergence of the presented finite element model for different numbers of different time step sizes and linearization algorithms. Displacements and temperature fields with phasechange interface illustrate an influence of the porosity change due to phase-change of pore water into ice on the deformation of the soils. From the presented numerical results for the basic thermo-elasticity model with phase-change of porous water to ice, we observe good convergence with a few numbers of the nonlinear iterations. We saw that influence of the time step size is larger for the displacements field than for the temperature field. Furthermore, we observe that more nonlinear iterations are needed for a larger time step size that converges to the one iteration per time layer. The algorithms with one Picard iteration provide worse results than the algorithm with nonlinear iterations, where errors are larger for less number of time layers. We observed a good numerical convergence of the presented implementation with the small number of nonlinear iterations that depends on time step size for two-and three-dimensional model problems.