Numerical model of soil freezing

wlcable - Mon, 2010-03-01 15:57

Following finite element framework (Zienkiewicz and Taylor, 1991), we approximate

$\begin{displaymath}T(x,\tau)\approx\sum_{i=1}^n\psi_i(x)\bm{t}_i(\tau),\end{displaymath}$

where $\bm{t}_i(\tau)$ is a ``value'' of temperature at the -th finite element grid, and $\psi_i$ is the so-called basis function. After some standard manipulations, we derive a system of differential equations

$\begin{displaymath} \bm{M}(\bm{t})\frac{d}{d\tau}\bm{t}(\tau)=-\bm{K}(\bm{t})\bm{t}(\tau), \ \ \ \ \bm{t}=\bm{t}(\tau) \end{displaymath}$

(10)

where $\bm{t}(\tau)=\{\bm{t}_i(\tau)\}_{i=1}^n$ is the vector consisting of temperature values, $\bm{M}(\bm{t})=\{m_{ij}(\bm{t})\}_{ij=1}^n$ and $\bm{K}(\bm{t})=\{k_{ij}(\bm{t})\}_{ij=1}^n$ are the $n{\times}n$ capacitance and stiffness matrices, respectively. A further refinement, which is often used in finite element modeling of phase change problems, is to exploit a so-called ``lumped'' formulation, i.e. the capacitance matrix $\bm{M}$ is diagonal:

$\begin{displaymath} m_{ij}(\bm{t})=\delta_{ij}c_i(\bm{t})\int_0^l\psi_idx, \ \... ...bm{t}){\approx}C(\bm{t}_i,x_i)+L\frac{d\theta}{dT}(\bm{t}_i), \end{displaymath}$

(11)

where $\delta_{ij}$ is one if , or zero otherwise. Dalhuijsen and Segal (1986) provides justification for the lumped formulation on noting that it is computationally advantageous and avoids oscillations in numerical solutions when used in conjunction with the backward Euler scheme:

$\begin{displaymath} \begin{array}{rlc} \big[\bm{M}_{}^{(k)}+d\tau_k\bm{K}_{}^... ...\ \bm{t}_{}^{(k)}&=\bm{t}_0, & \ \ \ k=0. \\ \end{array} \end{displaymath}$

(12)

The main difficulty in numerical modeling of soil freezing/thawing is in consistent calculation of the derivative $d\theta/dT$ in (11), where $\theta(T)$ is not a continuously differentiable function defined by (2). In many reviews, it is proposed to employ the enthalpy temporal averaging to calculate $c_i(\bm{t})$ . We suggest an approach that incorporates ideas of temporal averaging just to evaluate the rapidly changing $\theta(T)$ by defining as

$\begin{displaymath} c_i(\bm{t}_{}^{(k)})=C(\bm{t}_{i}^{(k)},x_i)+L\frac{\theta... ..._r(\bm{t}_{i}^{(k-1)})}{\bm{t}_{i}^{(k)}-\bm{t}_{i}^{(k-1)}}. \end{displaymath}$

(13)

We note that an advantage of this definition is that it does not compute temporal averaging of the heat capacity, and hence reduces numerical computations, and at the same time preserves numerical accuracy of the original idea. Even though that we avoid direct computation of the derivative $d\theta/dT$ by replacing it by the finite difference in (13), the calculation of $d\theta/dT$ is still necessary to obtain a solution of the so-called adjoint model used to construct the gradient of the cost function.

We emphasize again that $\theta(T)$ is not a continuously differentiable function, and hence some regularization of $\theta$ is required to avoid lost of accuracy in solving the adjoint problem. We propose the following regularization of the unfrozen water content $\theta(T)$ in (2):

$\begin{displaymath} \theta_r(T)=\eta\vert T_*\vert^b\vert r(T)\vert^{-b}, \ \ \... ...T-T_*)^{-2}}} & T<T_* \\ 0 & T{\ge}T_* \end{array}\right., \end{displaymath}$

where if $\alpha\to0$ then $r(T)\to\min(T,T_*)$ , and hence $\theta_r(T)\to\theta(T)$ . Here, is an infinitely differentiable function, and consequently $\theta$ posses the same property. One of the advantages of this is that we increase a convergence rate of solving the non-linear heat equation (12) by Richardson iterations, precisely find a solution to the adjoint problem, and hence find the location of the global minimum more faster than without regularization.