A phase field model for the freezing saturated porous medium

International Journal of Engineering Science 49 (2011) 768–780

Contents lists available at ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

A phase field model for the freezing saturated porous medium Jian-Fei Lu a,⇑, Yan-Pei Tan b, Jian-Hua Wang b a b

Department of Civil Engineering, Jiangsu University, Zhenjiang, Jiangsu 212013, PR China Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200030, PR China

a r t i c l e

i n f o

Article history: Received 22 February 2011 Accepted 24 March 2011 Available online 22 April 2011 Keywords: Porous medium Phase field Order parameter Mixture theory Memory effect

a b s t r a c t In terms of the mixture theory and phase field theory, a phase field model is developed for the saturated porous medium undergoing phase transition. In the proposed model, it is postulated that during the phase transition of the porous medium, both the solid skeleton and pore fluid will undergo phase transition. The phase states of the solid skeleton and pore fluid are characterized by respective order parameters. To simplify the proposed phase field model, the temperatures and order parameters of the solid skeleton and pore fluid are assumed to be equal. The balance laws of the porous medium are given by the conventional mixture theory. The order parameter and the porosity of the porous medium are considered as internal variables and their evolution equations are determined by the entropy inequality of the porous medium. The constitutive representations for the stresses, entropies, heat fluxes, drag force and the evolution equations for the order parameter and porosity are derived by exploitation of the entropy inequality. To illustrate the proposed model, a concrete phase field model for the freezing porous medium is established in the paper. In the model, the memory effect associated with phase transition of the porous medium is taken into account by assuming Stieltjes integral for the strain energy of the porous medium. The constitutive representations for the above variables are then derived according to the proposed free energy expression for the porous medium. Finally, the boundary condition associated with the proposed model and the determination of some parameters involved in our model are discussed in the paper briefly. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction When a fine-grained soil is freezing, a volume expansion named as frost heave which may damage the structures in contact with it usually occurs (Tsytovich, 1973, chaps. 1 and 2). Therefore, how to predict and control the heave of the freezing soil is a challenging problem for civil engineers. It is now widely recognized that the volume expansion of the freezing soil is mainly caused by the water transport from the unfrozen part of the soil to the freezing front rather than by the density difference between water and ice (Tsytovich, 1973). Furthermore, it has been acknowledged that the frozen fringe between the frozen region and unfrozen region (Akagawa, 1988; Gilpin, 1980; Miller, 1972; Radd & Oertle, 1973) is closely associated with the fluid migration in the freezing porous soil. Actually, the thickness of the film water surrounding the soil grains in the frozen fringe differs dramatically due to phase transition, which leads to a local fluid flow within the frozen fringe first, and the local fluid flow will in turn generate a long-range pressure gradient and thus a long-range pore fluid flow. Consequently, a correct phenomenological description of the drag force applied by the freezing film water to the pore fluid within the frozen fringe is very important for a successful description of fluid transport in the freezing porous medium.

⇑ Corresponding author. E-mail addresses: [email protected] (J.-F. Lu), [email protected] (Y.-P. Tan), [email protected] (J.-H. Wang). 0020-7225/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2011.03.014

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

769

For rocks, concretes and biological tissues, the stresses of the solid skeleton and pore fluid accumulated during the phase change process are directly relevant to their damage and fracture. For example, the stress state of pore ice of the saturated concrete is connected with the emergence of the internal microcracking and de-icer salt scaling (Marchand, Pleau, & Gagne, 1995); also, it is well-known that the thermal stresses developed in the freezing process of biological tissues may cause severe damage to biological tissues sometimes (Rabin & Steif, 2000). Moreover, for freezing porous medium, the stresses of the solid skeleton and pore fluid are also related to pore fluid transport. Consequently, to better understand the damage of the porous medium during phase transition and to have a complete description of pore fluid transport, it is crucial to have an exact evaluation of the stresses for the solid skeleton and pore fluid. However, as the moduli of the solid skeleton and pore fluid are heavily dependent on the degree of their phase transitions, thus, the stresses developed in the solid skeleton and pore fluid in phase transition process unavoidably depend on their strain histories. Therefore, a correct determination of the stresses of the freezing porous medium will inevitably involve memory effect. So far many freezing porous medium models have been established. Roughly, the freezing soil models can be categorized as the following three kinds: hydrodynamic models (Harlan, 1973; Sheppard, Kay, & Loch, 1978), rigid-ice models (Konrad & Morgenstern, 1980; O’Neill and Miller, 1985), and thermomechanical models (de Boer & Kowalski, 1995; Coussy, 2005; Duquennoi & Fremond, 1989; Fremond & Mikkola, 1991). Two severe shortcomings exist in the current models for the freezing porous medium. Firstly, the migration of the pore fluid resulting from the drag force from the freezing film water surrounding the solid grains is not phenomenologically accounted for properly. As noted previously, the remarkable thickness difference of the film water due to phase transition is the main cause for the fluid transport in the frozen fringe. Consequently, besides temperature gradient, the gradient of the degree of the phase transition of the film water also makes contribution to the drag force of the pore fluid in the frozen fringe. Nevertheless, the characterization of the gradient of the degree of the phase transition of the film water is beyond the capacity of the conventional freezing porous medium models. Secondly, the memory effect associated with the modulus variation of the solid skeleton and pore fluid during phase transition is ignored completely in the existing models. As the conventional freezing porous medium models are developed in the frame of the continuum mechanics and traditional thermodynamics, they lack effective method to specify the degree of the phase transition of the solid skeleton and pore fluid and relate the moduli of the solid skeleton and pore fluid to the degree of their phase transitions. As a result, it is difficult to account for the memory effect via conventional freezing porous medium models. To circumvent the above difficulties, in the context of the mixture theory, a phase field model is developed for the freezing saturated porous medium in this study. Central to the phase filed theory are the following important concepts (Boettinger, Warren, & Beckermann, 2002; Provatas & Elder, 2010). The phase transition of a material can be characterized by a order parameter representing the phase state of the material: 0 and 1 values of the order parameter represent the liquid and solid states, respectively, and the intermediate values denote the medium states; phase transition takes place within a layer with the order parameter taking on values between 0 and 1, rather than at a sharp interface. The most convincing physical justification for developing a phase field model for the freezing saturated soil is that the phase change of the saturated soil does occur in a transitional layer, i.e., in the above-mentioned frozen fringe (Akagawa, 1988; Gilpin, 1980; Miller, 1972; Radd & Oertle, 1973), which matches exactly with the above central conception of the phase field theory. Introducing the order parameter and its gradient into the expression for the free energy of the freezing porous medium within the context of the phase field theory will enable the gradient of the order parameter enter into the expression for the drag force. In this way, the above-mentioned local fluid flow due to the gradient of the degree of the phase transition of the film water can be accounted for phenomenologically. Moreover, with the order parameter, the memory effect of the porous medium during phase transition can be characterized by postulating the dependence of the moduli of the porous medium on the order parameter and introducing Stieltjes integral for the free energy and other variables in a straightforward way. Another advantage of the phase field model for the freezing porous medium is that the tracking of the freezing front for complex freezing problems is much more manageable than the conventional sharp interface model. Recently, the phase-field method has emerged as one of the most powerful methods for modeling many types of phase transition processes. Historically, the phase filed theory is based on a diffuse-interface description developed more than a century ago by Van der Waals (1893) and more than 50 years ago independently by Cahn and Hilliard (1958). In the phase field model, a free energy functional is constructed in terms of the phase order parameter as shown by Fix (1983), Langer (1986). Using the variational derivative of the free energy functional, the temporal evolution equations for non-conserved and conserved order parameters can be obtained, which are known as Cahn–Hilliard (1958) and Allen–Cahn (1979) equations, respectively. Conventionally, for the non-isothermal case, the temporal evolution equation for the order parameter have to be combined with a heat transfer equation obtained by adding an liberation of latent heat term to the standard heat transfer equation. As the above-mentioned heat transfer equation is not based on the thermodynamic arguments, thus, the combination of the evolution equation for the order parameter and heat transfer equation cannot guarantee a monotonically decreasing free energy functional. As pointed out by Penrose (1990) and Wang et al. (1993), the standard evolution equation for the order parameter is derived from a free energy functional that is only applicable to the isothermal case. Thus, to circumvent the difficulty of the phase field model, Penrose (1990) and Wang et al. (1993) introduced the functional of entropy and derived the governing equations for the phase field model via the entropy functional. More recently, the phase field theory has been extended to account for the convection of fluid by Anderson, McFadden, and Wheeler (2000) and Jeong, Goldenfeld, and Dantzig (2001). To develop phase-field models within the context of continuum physics, Fried and Gurtin (1993, 1996) introduce a balance law for microforces which are the forces expending power with the variation of the order

770

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

parameter. Along the same line, Fabrizio et al. (2006) recently put forward a general continuum phase field model from basic laws of continuum mechanics and thermodynamics. In this theory, the phase field order parameter is treated as an internal variable of the constitutive functions and the evolution equation for the order parameter is derived from the entropy inequality. It should be stressed that treating of the order parameter as an internal variable and developing its phase field theory in the context of continuum physics is pivotal for the development of the phase field theory. In this manner, the phase field model can be established on a more rigorous theoretical basis. Also, the resulting phase field models are more general than the conventional phase filed models, and various coupling effects can thus be considered in phase field models in a systematic and thermodynamically consistent way. In the proposed phase field model for the freezing porous medium, as the film water surrounding the solid grains behaves quite differently from the free-flow pore water, it is thus viewed as a constituent of the solid skeleton. Moreover, as the phase transition of the film water entails dramatic variation of the properties of the solid skeleton, and furthermore, it is relevant to the pore fluid transport during phase transition, consequently, it is necessary for the solid frame to be endowed with an order parameter to specify its phase state. For the pore fluid, all the pore ‘‘fluids’’ including the pore fluid in the unfrozen region, the pore ice in the frozen region, as well as the pore fluid and pore ice in the frozen fringe are treated as one kind of ‘‘solid’’ (one phase) with its phase state specified by an order parameter. Although in the proposed model, pore fluid and solid skeleton possess their own order parameters and temperatures, respectively, to simplify the proposed phase field model, the order parameters and temperatures for the solid skeleton and pore fluid are assumed to be equal (Bowen, 1980, 1982). We follow the work of Fabrizio et al. (2006) and treat the order parameter for the solid skeleton and pore fluid as an internal variable. To describe the microstructure of the porous medium, the porosity of the porous medium is also introduced as an internal variable (Bedford & Drumheller, 1983; Bowen, 1982). The balance laws of the porous medium are given by the conventional mixture theory (Atkin & Craine, 1976; Bowen, 1982). The constitutive representations for the stresses, entropies, heat fluxes, drag force and the evolution equations for the order parameter and porosity are derived by exploitation of the entropy inequality of the porous medium. To illustrate the proposed model, a concrete phase field model for the freezing porous medium is established. In the model, free energy expressions with the memory effect associated with phase transition of the solid skeleton and pore fluid are constructed first. Accordingly, the above-mentioned constitutive functions and evolution equations are derived in terms of the free energies. Finally, the boundary condition associated with the proposed model and the determination of some parameters involved in our model are discussed briefly. 2. Kinematics and balance laws for the freezing porous medium In this paper, the porous medium is considered as consisting of the solid skeleton and the pore fluid. As noted above, the film water surrounding the solid grains is viewed as one part of the solid skeleton, so the solid skeleton here differs from that of the common porous medium. Obviously, the film water’s inclusion into the solid skeleton makes its properties highly dependent on the degree of phase transition, and the solid skeleton here is thus treated as a special kind of material with its properties characterized by an order parameter. The pore fluid here refers to the free-flow pore fluid of the saturated porous medium. As phase transition may occur in the pore fluid, hence, the pore fluid is regarded as a special kind of ‘‘solid’’ endowed with an order parameter. Consequently, both the solid skeleton and pore fluid here are different from those of the conventional porous medium free of phase transition. In what follows, the above solid skeleton is identified by the superscript ‘‘s’’ and the free-flow pore fluid is by the superscript ‘‘f’’; variable a is used to denote an individual bulk component, including both the solid and fluid components. The balance laws of the two components of the porous medium are formulated in the framework of the generalized mixture theory (Atkin & Craine, 1976; Bowen, 1980, 1982). Generally, the two components s and f move independently in the course of a thermodynamic process. Thus, the material particles of each component should be characterized by independent initial material position XðaÞ , and the motion of the material points of the a component, therefore, may be expressed as

^ðaÞ ðXðaÞ ; tÞ; x¼x

a ¼ s; f ;

ð1Þ

^ðaÞ ðXðaÞ ; tÞ is assumed to be invertible and differentiable as many times as necessary. The relative velocity of where function x the pore fluid with respect to the solid skeleton is defined as

vðf ;sÞ ¼ v ðf Þ  vðsÞ :

ð2Þ

The deformation gradient of the a component is defined by

FðaÞ ¼

^ðaÞ ðXðaÞ ; tÞ @x @XðaÞ

^ðaÞ ðXðaÞ ; tÞ; ¼ GRADx

a ¼ s; f ;

ð3Þ

where GRAD denotes the gradient with respect to the reference coordinate X(a). The Green strain tensor of the a component is represented by the following expression

EðaÞ ¼

1 ðaÞT ðaÞ ðF F  IÞ; 2

a ¼ s; f :

ð4Þ

In our model, the mass exchange between two components is excluded. Thus, for thermodynamic processes of the porous medium under consideration, the following mass balance equation for the a component is obtained

771

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

DðaÞ qðaÞ þ qðaÞ r  v ðaÞ ¼ 0;

a ¼ s; f ;

ð5Þ

(a)

( a)

where D denotes the material derivative with respect to the motion of the a component, v ponent. Similarly, the momentum balance equation for the a component has the form

is the velocity of the a com-

qðaÞ DðaÞ v ðaÞ  r  rðaÞ  qðaÞ bðaÞ ¼ TðaÞ ; a ¼ s; f ;

ð6Þ

where b(a) denotes external supply of linear momentum per unit mass, r(a) is the bulk stress of the a component and TðaÞ represent the rate of linear momentum transfer to the a component from the other component. For the porous medium under consideration, the momentum conservation law asserts that

Tðf Þ ¼ TðsÞ :

ð7Þ T

Moreover, assuming that there is no moment of momentum transfer between the components, it follows that rðaÞ ¼ rðaÞ . It is noted that the barycentric velocity of the porous medium is avoided in Eqs. (5) and (6), and thus, the diffusion term will not appear in the total stress of the porous medium (Bowen, 1980, 1982). The energy balance equation of the a component is as follows

qðaÞ DðaÞ eðaÞ þ r  qðaÞ  qðaÞ rðaÞ  rðaÞ : LðaÞ ¼ Q ðaÞ ; a ¼ s; f ; ( a)

ð8Þ (a)

(a)

where e denotes the internal energy density of the a component, q is the heat flux of the a component, r is the external ðaÞ ðaÞ  ðaÞ energy supply to the a component, L(a) is the velocity gradient of the a component, and L(a) = v(a)r, rðaÞ : LðaÞ ¼ rij Lij , Q is the rate of energy transfer to the a component from the other component. The energy conservation for the saturated porous medium implies the following relation

Wþ

X

 ðaÞ ¼ 0; Q

W :¼ Tðf Þ  v ðf ;sÞ ;

ð9Þ

a

where W denotes the work of the drag force between the solid skeleton and pore fluid. As noted previously, for the porous medium undergoing phase transition, the temperatures of the solid skeleton and pore fluid are assumed to be equal during the thermodynamic process (Bowen, 1982). Consequently, there is only one independent temperature for the porous medium and thus, energy balance equations (8) for two components can be added to form one independent energy balance equation. Thus, the Helmholtz free energy w(a) of the a component is defined by the following Legendre transform

wðaÞ ¼ eðaÞ  hgðaÞ ; where h, g

(a)

a ¼ s; f ;

ð10Þ

are the temperature and entropy of the a component, respectively.

3. Entropy inequality and constitutive relations for the freezing porous medium As the balance laws presented in the previous section are insufficient to characterize the system, thus, to obtain a closed system of equations, it is necessary to express some physical variables (the dependent variables) in terms of constitutive functions of the independent variables. In this study, the constitutive equations for the dependent variables are derived by the exploitation of the entropy inequality of the porous medium (Coleman & Noll, 1963). To describe the microstructure of the porous medium, the porosity (/) of the porous medium is also considered as a variable of the porous medium system. The fraction of the microscopic volume occupied by the fluid component is thus represented by the porosity n(f) = /, while the fraction of the microscopic volume occupied by the solid skeleton is n(s) = 1  /. Each individual component therefore, can be characterized by two densities, one of which is the true density (c(a)) relative to the volume it actually occupies, and the other is the bulk density (q(a) = n(a)c(a)) relative to the total volume occupied by the mixture of the porous medium. Introducing the porosity of the porous medium as a variable will inevitably entail the closure problem for the porous medium (Bedford & Drumheller, 1983), as there is no corresponding balance law for the porosity or volume fraction of each component. To close the system, two options are available now. The first method resolves the closure problem by introducing additional balance law for the volume fractions (Ahmadi, 1982; Passmann et al., 1984), while the second approach closes the system by treating volume fractions as internal variables (Bowen, 1982; Morland, 1972) and assuming relaxation type constitutive equations for volume fractions. In this study, we follow the approach of Morland (1972) and Bowen (1982) and treat the porosity as an internal variable to derive its relaxation type constitutive equation. To account for the phase transition of the porous medium, as noted previously, the order parameters for the solid skeleton and pore fluid are introduced as the internal variables for the system. Also, to fully accommodate the interfacial energy of the frozen fringe in the freezing porous medium, it is assumed that the free energy of the porous medium is dependent on the gradient of the order parameter as well. As the temperatures of the solid skeleton and pore fluid are assumed to be equal, thus, it is reasonable to suppose that the order parameters for the solid skeleton and pore fluid are also equal. Therefore, only one order parameter is independent in our model.

772

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

When deriving the constitutive equations for the porous medium under phase transition, the following variables

C ¼ fEðaÞ ; /; h; rh; v; rv; Dv; v ðf ;sÞ g;

a ¼ s; f

ð11Þ

are chosen as the independent state variables. Herein, v is the order parameter for the solid skeleton and pore fluid, r and D denote the gradient operator and Laplacian. Moreover, as the densities of the solid skeleton and pore fluid can be determined using the strain tensors and mass balance laws, they are not included in the set of the independent variables. According to the above independent variables, the constitutive relations for the free energy, entropy, partial stress, heat flux, extra entropy flux and the drag force between the pore fluid and solid skeleton are assumed as follows

^ ðaÞ fEðaÞ ; h; /; v; rvg; wðaÞ ¼ w

gðaÞ ¼ g^ ðaÞ fEðaÞ ; h; /; v; rvg; rðaÞ ¼ r^ ðaÞ fEðaÞ ; h; /; v; rvg; ð12Þ

^ ðaÞ fh; rh; /; v; rv; v ðf ;sÞ g; qðaÞ ¼ q ðaÞ

k

^ ðaÞ fEðsÞ ; Eðf Þ ; h; rh; /; v; rvg; ¼k

^ ðf Þ fEðsÞ ; Eðf Þ ; h; /; v; rv; Dv; v ðf ;sÞ g; Tðf Þ ¼ T

a ¼ s; f ;

where k(a) is the extra entropy flux, which is due to the non-locality of the free energy (Fabrizio et al., 2006; Fosdick & Rajagopal, 1980; Wang et al., 1993). It is worthwhile noting that the curly bracket in above equation is different from the parenthesis used in the common constitutive assumptions, and it indicates that the variables listed above may depend on the histories of their arguments. However, it is assumed in our model, that the current free energy only depends on the current gradient of the order parameter. Also, it is noted that, for simplicity, in the above constitutive assumption, the free energies of the solid skeleton and pore fluid are postulated only to depend on the independent variables of respective phase, which violates the principles of equipresence and determinism (Truesdell & Toupin, 1960) and thus is different from the assumptions of Atkin and Craine (1976) and Bowen (1980, 1982) etc. However, it is well recognized that for mixtures consisting of immiscible constituents, some variables of each phase depend only upon the independent variables associated with the phase (Drumheller & Bedford, 1980; Eringen, 1994; Hassanizadeh & Gray, 1990; Thigpen & Berryman, 1986). As mentioned above, the order parameter and the porosity are considered as internal variables of the thermodynamic system, so to close the system, their temporal evolution equations should be derived via entropy inequality as well. Thus, the following relaxation type constitutive equations are assumed for the order parameter and porosity

^ ðvaÞ fEðaÞ ; h; /; v; rvg; DðaÞ v ¼ xðvaÞ ¼ x ðaÞ

ðaÞ

^ / fEðaÞ ; h; /; v; rvg; DðaÞ / ¼ x/ ¼ x

ð13Þ

a ¼ s; f :

Note that the material derivatives of an arbitrary variable and its gradient with respect to the motions of the solid skeleton and pore fluid have the following relations

DðsÞ g ¼ Dðf Þ g  v ðf ;sÞ  rg; DðsÞ rg ¼ Dðf Þ rg  ðrrgÞ  v ðf ;sÞ ¼ rDðf Þ g  rg  Lðf Þ  ðrrgÞ  v ðf ;sÞ ;

ð14Þ

where g is an arbitrary variable. In view of (14), it can be concluded that D(s)v and D(f)v as well as D(s)/ and D(f)/ are not independent. In the proposed model, we chose to derive the evolution equations for D(f)v and D(s)/, respectively. In this study, the second law of thermodynamics assumes the form of the Clausius–Duhem inequality. Thus, all the solutions of the balance equations and constitutive equations must satisfy the following Clausius–Duhem inequality (Coleman & Noll, 1963)



X

qðaÞ ðDðaÞ w^ ðaÞ þ gðaÞ DðaÞ hÞ þ

a

X

rðaÞ : LðaÞ 

a

X X 1 X ðaÞ  ðaÞ P 0: q  rh þ h r  kðaÞ þ Q h a a a

ð15Þ

The above Clausius–Duhem inequality differs from the standard Clausius–Duhem inequality for mixtures (Atkin and Craine, 1976) by the presence of the extra entropy flux k(a) (Fabrizio et al., 2006). The total stress, total mass density of the saturated porous medium, the total heat flux and total extra entropy flux are defined as follows:

q¼

X a

qðaÞ ; r ¼

X a

rðaÞ ; q ¼

X a

qðaÞ ;

k¼

X a

ðaÞ

k :

ð16Þ

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

773

By using constitutive assumption (12) as well as Eqs. (9) and (14), the inequality (15) leads to

" #    ^ ðrÞ @w ðr Þ ^ ðrÞ kh þ qðrÞ gðrÞ ÞDðf Þ h þ rðrÞ  qðrÞ FðrÞ  w ^ ðrÞ  : L ðr Þ  ðqðrÞ w  FðrÞT þ rv  qðrÞ E @ rv     ^ ðrÞ @w ðsÞ ^ ðr Þ  ^ ðrÞ   qðrÞ w  rxðvrÞ / x/  qðrÞ ðw vÞxðfv Þ  qðrÞ @ rv   ^ ðf Þ  ^ ðsÞ   ½Tðf Þ þ qðf Þ ðw /Þr/  qðsÞ ðw vÞrv  v ðf ;sÞ

ð17Þ

1 ^ ðsÞ kh þ qðsÞ gðsÞ Þv ðf ;sÞ   rh þ hr  k P 0;  ½q  hðqðsÞ w h wherein and in what follows, the terms with repeated (not limited to twice) superscript ‘‘ r ’’ denote the summation over the  ^ ðrÞ kh ¼ qðsÞ @ w ^ ðsÞ kh þ qðf Þ @ w ^ ðf Þ kh; the symbol w ^ ðaÞ  solid and fluid components, for example, qðrÞ @ w - (a = s, f) denotes the coef ðaÞ ^ ðaÞ ðaÞ  ^ ficient of the material derivative of - in the expanded expression for D w . Obviously, w - is different from the conven^ ðaÞ explicitly. However, if w ^ ðaÞ only has a parametric dependence on tional partial derivative, as -(t) may not be contained in w  ðaÞ  ðaÞ ^ ^ -, then w - is equivalent to @ w =@ -. Since the above inequality is linear with respect to D(f)h and L(a), thus, the following constitutive equations for the entropy and stress are obtained

qðrÞ gðrÞ ¼ qðrÞ w^ ðrÞ kh;

" !#    ^ ðaÞ ðaÞT ðþÞ ðaÞ  ðaÞ ðaÞ @ w ^ r ¼ q F  w E  F  S rv  q ; @ rv " !# ^ ðaÞ @w ¼ 0; a ¼ s; f SðÞ rv  qðaÞ @ rv ðaÞ ðaÞ

ðaÞ

ð18Þ

in which functions SðþÞ ½ and SðÞ ½ represent the operations for taking the symmetric and skew parts of their second order tensor arguments. Note that when derive the above constitutive relations, the symmetry of the stress r(a) and    qðaÞ FðaÞ  w^ ðaÞ EðaÞ  FðaÞT is used. Furthermore, Eq. (18)3 implies the following relation

^ ðaÞ @w ¼ jðaÞ hrv; @ rv

a ¼ s; f ;

ð19Þ

where the constant j(a) is associated with the thicknesses and the surface energy of the transition layers of the freezing a component material (Boettinger et al., 2002). Using Eq. (19) in Eq. (18)2, the expression for the stress r(a) can be simplified as follows:



rðaÞ ¼ qðaÞ FðaÞ  ðw^ ðaÞ EðaÞ Þ  FðaÞT  qðaÞ jðaÞ hðrv  rvÞ; a ¼ s; f :

ð20Þ

Using Eq. (18), the residual entropy inequality, i.e., the entropy production of the porous medium is obtained as follows:

" ! #  ^ ðsÞ 1 ðsÞ ðrÞ  ðsÞ @ w ðsÞ ðsÞ ðf ;sÞ ^  rh qh q  q ðw /Þx/  v þq g h @h   h    i  ^ ðf Þ  ^ ðsÞ   Tðf Þ þ qðf Þ w v rv  v ðf ;sÞ / r/  qðsÞ w ðr Þ

ð21Þ

   ^ ðrÞ @w ^ ðr Þ   qðrÞ w  rxðvrÞ þ hr  k P 0: v xðfv Þ  qðrÞ @ rv The above inequality can be further reduced to

" ! #    ^ ðsÞ 1 @w ðsÞ ^ ðr Þ  q  h qðsÞ  qðrÞ w þ qðsÞ gðsÞ v ðf ;sÞ  rh / x/  h @h " ! #       ^ ðsÞ @w ^ ðf Þ  ^ ðsÞ   Tðf Þ þ qðf Þ w rv  v ðf ;sÞ / r/  qðsÞ w v rv þ r  qðsÞ @ rv ! " !#    ^ ðrÞ ^ ðr Þ @w @w  0: þ r  hk  qðrÞ xðvrÞ  k  rh  xðfv Þ qðrÞ w^ ðrÞ v  r  qðrÞ @ rv @ rv

ð22Þ

If assuming the extra entropy flux has the following expression (Fabrizio et al., 2006; Wang et al., 1993)

!

qðrÞ @ w^ ðrÞ ðrÞ qðrÞ @ w^ ðrÞ ðf Þ qðsÞ @ w^ ðsÞ k¼ x ¼ x   rv  v ðf ;sÞ h @ rv v h @ rv v h @ rv

ð23Þ

774

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

and assuming the entropy for each component has the following form (Hassanizadeh & Gray, 1990; Wei & Muraleetharan, 2002)

qðaÞ gðaÞ ¼ qðaÞ w^ ðaÞ kh; a ¼ s; f

ð24Þ

then, the entropy production of the porous medium has the following simplified representation

" # " !#    ^ ðsÞ   qðrÞ  ^ ðrÞ   qðrÞ @ w^ ðrÞ 1 ðsÞ ðrÞ  ðf Þ ðsÞ @ w ðf ;sÞ ^  rh w / x/  h q xv  q  q rv  v w v  r  h h h @ rv @ rv " #       ^ ðsÞ @w ^ ðf Þ  ^ ðsÞ   Tðf Þ þ qðf Þ w Þrv  v ðf ;sÞ P 0: / r/  qðsÞ w v rv þ r  ðqðsÞ @ rv ðrÞ

ð25Þ

When Eq. (24) is supposed to hold, for the entropy inequality to assume the standard form for the entire concerned region, the following boundary condition should be fulfilled

Z

k  ds ¼ 0;

ð26Þ

@X

where oX is the boundary of the concerned region. If the porous medium under phase transition is assumed to deviate from the equilibrium state by a small amount, then, the linear expansion of the ‘‘fluxes’’ about the ‘‘forces’’ can be made. According to Curie-Prigogine principle (de Groot & Mazur, 1984), when making linear expansion of a ‘‘flux’’ about the ‘‘forces’’, the dependency of the ‘‘flux’’ tensor on the ‘‘force’’ tensors with odd number order differences between the ‘‘flux’’ tensor can be excluded. Consequently, in Eq. (25), ðsÞ ðf Þ the ‘‘fluxes’’ corresponding to the ‘‘forces’’ rh and v(f,s) cannot depend on the ‘‘forces’’ associated with x/ and xv . Likewise, ðsÞ ðf Þ ðsÞ the ‘‘fluxes’’ x/ and xv cannot depend on the ‘‘forces’’ rh and v(f,s). Also, for simplicity, the dependency of the ‘‘flux’’ x/ ðf Þ on the ‘‘force’’ of xv and vice versa are excluded in this paper. As a result, the constitutive relations for the time change rate of the order parameter, the porosity, the heat flux and drag force between the solid skeleton and pore fluid can be represented as follows:

 ðrÞ     q ^ ðrÞ  w v  r  ðqðrÞ jðrÞ rvÞ ; h    ðrÞ ^ ðrÞ  w / ; ¼ M / q

xðfv Þ ¼ Mv xðsÞ /

ð27Þ

q ¼ kh rh  ktv hv ðf ;sÞ þ qðsÞ jðsÞ hðrv  v ðf ;sÞ Þrv;         ^ ðf Þ  ^ ðsÞ  Tðf Þ ¼ bv v ðf ;sÞ  kv t rh  qðf Þ w / r/ þ qðsÞ w v rv  r  qðsÞ jðsÞ hrv rv; where Eq. (19) has been used in the derivation of the above equations. In (25), to guarantee the non-negative entropy production, constants Mv, M / , kh and bv should be non-negative. Constants Mv and M / are the coefficients for the relaxation of the order parameter and the porosity. Constant kh is the conventional coefficient of heat conduct and constant bv is determined by the permeability of the porous medium. Also, the Onsager reciprocity relation (de Groot & Mazur, 1984) and the nonnegative entropy production of the porous medium imply the following relations for the coefficients kh, ktv, kvt and bv

ktv ¼ kv t ;

2

kh bv  hkv t > 0:

ð28Þ

It is worth emphasizing that the constitutive relations given by Eqs. (18), (20), and (27) fulfill the frame indifference principle automatically (Truesdell & Noll, 1965). In addition, Eq. (27) shows that besides the contribution from the conventional diffusion-thermal coupling effect (Chapman & Cowling, 1970), the heat flux vector also involves the gradient of the order parameter and furthermore, it is also influenced by the scalar product of the gradient of the order parameter and the relative velocity v(f,s). Moreover, besides involving the conventional relative velocity term, gradient terms of the porosity and temperature (Bowen, 1980, 1982), the drag force also contains two gradient terms of the order parameter, one of which is associated with the non-locality of the free energy. As noted previously, the gradient terms of the order parameter can account for the drag force from the freezing film fluid in the frozen fringe phenomenologically, which is beyond the capacity of the existing freezing porous medium models. 4. A specific phase field model for the freezing porous medium Based on the general theory developed in the above section, a specific model for the saturated porous medium undergoing phase transition will be developed in this section. In the model, specific expressions for the free energy of the porous medium are constructed first. The expressions for other dependent variables can then be derived accordingly. For the thermodynamic process of the porous medium under consideration here, we assume the following reference state

n o n o ðaÞ EðaÞ ; cðaÞ ; /; v ðaÞ ; h; v ¼ 0; c0 ; /0 ; 0; h0 ; v0 ;

a ¼ s; f ;

ð29Þ

775

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780 ðaÞ

where c0 , /0 , h0, v0 are the reference density, porosity, temperature and order parameter for the initial state. When a disturbance is applied, the porous medium will deviate from the reference state and arrive at a new state. For the new state, it is assumed that the free energies for the solid skeleton and pore fluid can be expressed by the additive forms. Consequently, the free energy for the a component has the following form ðaÞ

ðaÞ

ðaÞ

ðaÞ

ðaÞ

ðaÞ

ðaÞ

ðaÞ

^ fEðaÞ ; /; v; hg þ w ^ fv; hg þ w fv; hg; ¼w M T I

wðaÞ ¼ wM þ wT þ wI

a ¼ s; f ;

ð30Þ

ðaÞ

where wM ; wT and wI denote the mechanical, thermal and interface free energies, respectively. Although strain energy has been taken into account in some phase field models (Chen, 2002), the memory effect associated with the variation of moduli due to phase transition has not been considered in the current literature so far. However, as noted previously, for the porous medium under phase transition, its moduli vary remarkably with the progress of phase transition. As a result, the same strain increment at different stages of phase transition will have different contribution to stress. Thus, to correctly evaluate the stress developed in the freezing porous medium, it is crucial to take into account of the memory effect associated with the variation of the moduli. It is noted that in this study only the memory effect associated with the variation of elastic moduli due to the phase transition is considered, while the memory effect related to the viscoelasticity (Christensen, 1982) of the porous medium is ignored. As the mechanical part free energy of each component of the porous medium involves above-mentioned memory effect, it thus can be considered as a functional of its independent variables. According to the Stone–Weierstrass theorem (Rudin, 1973), this part of free energy can be approximated uniformly by a polynomial in a set of real continuous linear scalar valued functionals of their variables (Chacon & Rivlin, 1964; Christensen, 1982). Therefore, for the isotropic porous medium, the mechanical part free energy for the a component is approximated by the following Stieltjes integral expression

qðaÞ wðMaÞ ¼

Z t Z

1 2J ðaÞ  þ þ

1 J

0

1

0

Z t Z

ða Þ

 PðaÞ ðvðs1 Þ; vðs2 ÞÞ : DðaÞ EðaÞ ðs2 Þds2 : DðaÞ EðaÞ ðs1 Þds1

t

0

t

0

Z t Z

J ða Þ 1 2J ðaÞ

0

t

0

Z t Z 0

0



h

i

PðvaÞ ðvðs1 Þ; vðs2 ÞÞDðaÞ trðEðaÞ ðs2 ÞÞds2 DðaÞ tr EpðaÞ ðvðs1 ÞÞ þ EhðaÞ ðhðs1 ÞÞ ds1 h

i



Pðe/aÞ ðvðs1 Þ; vðs2 ÞÞDðaÞ tr EðaÞ ðs2 Þ  EpðaÞ ðvðs2 ÞÞ  EðhaÞ ðhðs2 ÞÞ ds2 DðaÞ /ðs1 Þds1 t



aÞ Pð// ðvðs1 Þ; vðs2 ÞÞDðaÞ /ðs2 Þds2 DðaÞ /ðs1 Þds1 ; a ¼ s; f ;

ð31Þ

where J(a) = |F(a)|, P(a) is the fourth order tensor related to the elastic moduli of the component; PðvaÞ is the material constant ðaÞ accounting for the bulk modulus of the a component; Pe/ is the material parameter related to the moduli accounting for the ðaÞ coupling between the strain and the variation of the porosity; P// is the material parameter accounting for the modulus for ðaÞ ðaÞ the variation of the porosity; Ep ; Eh are the strains due to phase transition and temperature variation. ðaÞ For the isotropic porous medium considered here, the strains EpðaÞ ; Eh due to phase transition and temperature change have the following expressions

i i 1 h ðaÞ2 1 h ðaÞ2 ðaÞ kp ðvÞ  1 I; Eh ¼ kh ðhÞ  1 I; 2 2 h i1=3 ðaÞ ðaÞ ðaÞ ; kh ðhÞ ¼ aðaÞ ðh  h0 Þ þ 1; kp ðvÞ ¼ pðvÞep þ 1   pðvÞ ¼ v3 6v2  15v þ 10 ; a ¼ s; f ;

EðpaÞ ðvÞ ¼

ð32Þ

ðaÞ

where ep is the volume change of per unit volume of the a component when it changes phase from the fluid state to the solid state, I is the identity tensor, a(a) is the coefficient of linear thermal expansion of the a component, h0 is the initial temperature of the porous medium and p(v) is the conventional interpolating function (Boettinger et al., 2002). For simplicity, the coefficient of linear thermal expansion a(a) here is viewed as a constant. As noted above, v = 0 and v = 1 indicate the fluid phase and solid phase, respectively. Also, the interpolating function p(v) vanishes at v = 0, and it equals to one at v = 1. ðaÞ It is worthwhile noting that, for simplicity, the coupling between EðpaÞ and Eh is neglected in Eq. (31). Also, to make the ðaÞ free energy has local minima at v = 0, 1, the derivative of Ep ðvÞ with respect to v should vanish at v = 0, 1. This requirement is fulfilled by the presence of the function p(v) in EðpaÞ ðvÞ as shown in Eq. (32). Besides, as the higher order terms are omitted in above representation for the free energy, thus, the above free energy representation is only applicable to the case for small deformation with large displacement. The fourth order tensor P(a) in (31) fulfils the following restriction aÞ ða Þ PðIJKL ðvðs1 Þ; vðs2 ÞÞ ¼ PKLIJ ðvðs2 Þ; vðs1 ÞÞ;

a ¼ s; f ;

ð33Þ

where the capital subscripts denote the components with respect to the Lagrangian (material) coordinate system. Also, to fully determine the stress, it is assumed that P(a) satisfies the following equation

776

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

aÞ ðaÞ ðaÞ PðIJKL ðvðtÞ; vðsÞÞ ¼ PKLIJ ðvðsÞ; vðtÞÞ ¼ C IJKL ðvðsÞÞ;

a ¼ s; f ;

ð34Þ

ðaÞ C IJKL

where in above equation is similar to the elastic tensor for the isotropic Kirchoff material (Truesdell & Noll, 1965) in the reference configuration. In the abstract notation, the elastic tensor for the a component has the following expression ð1Þ

CðaÞ ðvðsÞÞ ¼ kðaÞ ðvðsÞÞI  I þ 2lðaÞ ðvðsÞÞ I ; ð1Þ

I ¼

1 ðdIK dJL þ dIL dJK ÞeI  eJ  eK  eL ; 2

ð35Þ

a ¼ s; f ; ðaÞ

ðaÞ

where eI, eJ, eK and eL are the base vectors for the Lagrangian (material) coordinate system. Also, PðvaÞ , Pe/ and P// satisfy the following relations

PðvaÞ ðvðs1 Þ; vðs2 ÞÞ ¼ PðvaÞ ðvðs2 Þ; vðs1 ÞÞ;

ðaÞ ðaÞ Pe/ ðvðs1 Þ; vðs2 ÞÞ ¼ Pe/ ðvðs2 Þ; vðs1 ÞÞ;

aÞ ðaÞ Pð// ðvðs1 Þ; vðs2 ÞÞ ¼ P// ðvðs2 Þ; vðs1 ÞÞ;

2 3

PðvaÞ ðvðtÞ; vðsÞÞ ¼ PðvaÞ ðvðsÞ; vðtÞÞ ¼ K ðaÞ ðvðsÞÞ ¼ kðaÞ ðvðsÞÞ þ lðaÞ ðvðsÞÞ;

ð36Þ

ða Þ ðaÞ Pðe/aÞ ðvðtÞ; vðsÞÞ ¼ Pe/ ðvðsÞ; vðtÞÞ ¼ K e/ ðvðsÞÞ;

aÞ ða Þ ðaÞ Pð// ðvðtÞ; vðsÞÞ ¼ P// ðvðsÞ; vðtÞÞ ¼ K // ðvðsÞÞ;

a ¼ s; f ; ðaÞ

where K(a) is the bulk modulus of the a component, K e/ is the modulus accounting for the coupling between the strain and ðaÞ the variation of the porosity and K // is that for the variation of the porosity. In this paper, we exclude the thermal memory effect as considered by Gurtin and Pipkin (1968) and Golden (2008), so the thermal and interfacial part free energies can assume the standard forms. The thermal part of the free energy in Eq. (30) is the free energy of the pure a component free of deformation. Thus, free energy models in various existing phase field models can be borrowed. Many kinds of free energy model for freezing pure materials have been proposed so far (Boettinger et al., 2002). In this study, the following representation is assumed for the thermal part free energy of a unit mass of the a component (Wang et al., 1993) ðaÞ wT

" Z ðaÞ ^ ¼ wT ðv; hÞ ¼ h 

eðTfaÞ ðnÞ

hM

2

2

gðvÞ ¼ v ð1  vÞ ; ðaÞ Tf ðnÞ

h

n2

dn þ pðvÞ

Z

h

hM

LðaÞ ðnÞ n2

# dn þ W

ðaÞ

gðvÞ ;

ð37Þ

a ¼ s; f ;

( a)

where e and L (n) are thermal part internal energy for the fluid phase of the a component and the latent heat of the a component, W(a) is the potential height between the two states with minimum free energy, g(v) is the common double-well function, and hM is the melting temperature of the porous medium. Note that the melting temperature of the porous medium hM may depend on the pressure of the material under phase transition. However, we postulate that the stress variation of the porous medium is small compared with the pressure ranges in the phase diagrams of the solid skeleton and pore fluid and hence, it is thus treated as a constant here. In many phase field models, the latent heat is assumed to be constant (Boettinger et al., 2002; Wang et al., 1993). However, when dealing with materials with large heat capacity difference between the fluid and solid phase, it is necessary to account for the variation of its latent heat with temperature. In our model, the expressions for the thermal internal energy of the fluid phase and the latent heat are approximated as follows:

eðTfaÞ ðhÞ ¼ cðEfaÞ ðh  hM Þ;

  ðaÞ ðaÞ ða Þ ðaÞ ðaÞ LðaÞ ðhÞ ¼ LM þ cEf  cEs ðh  hM Þ ¼ LM þ DcE ðh  hM Þ;

ð38Þ

a ¼ s; f ; ðaÞ

where the internal energy at the reference state, i.e., the melting point, has been set to be zero, LM is the latent heat at the ðaÞ ðaÞ melting temperature, cEf ; cEs are the heat capacities of the fluid and solid phases of the a component at constant strain, ðaÞ ðaÞ ðaÞ which are assumed to be constants whether for the supercooled fluid or overheated solid, and DcE ¼ cEf  cEs . Substituting ðaÞ Eq. (38) for eTf and L(a) into Eq. (37), the thermal part free energy for per unit mass of the a component is obtained as follows:

h i ðaÞ ^ ðaÞ ðv; hÞ ¼ h cðaÞ ðvÞv ðaÞ ðhÞ þ bðaÞ ðvÞuðaÞ ðhÞ þ W ðaÞ gðvÞ ; wT ¼ w T ðaÞ

L ða Þ ða Þ ðaÞ ðaÞ ðaÞ cðaÞ ðvÞ ¼ cEf  pðvÞDcE ; b ðvÞ ¼ cEf þ pðvÞ M  pðvÞDcE ; hM  h h  hM v ðaÞ ðhÞ ¼ ln M ; uðaÞ ðhÞ ¼ ; a ¼ s; f : h h

ð39Þ

The advantage for using the above free energy model is that except W(a), all parameters can be determined by material constants of the freezing material. More importantly, it fulfils the fundamental requirements for the free energy representation of a freezing material. As the first, second derivatives of p(v), and the first derivative of g(v) vanish at v = 0, 1, while the second derivative of g(v) is positive at v = 0, 1, thus, the free energy given by Eq. (39) will have local minima at v = 0, 1.

777

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

Also, it can be shown that when h = hM, the two minima are equal; when h > hM, the free energy has the absolute minimum at v = 0; when h > hM, the free energy has the absolute minimum at v = 1. For the interface free energy of the pure a component material, the following expression is assumed (Boettinger et al., 2002; Provatas & Elder, 2010) ðaÞ

wI ðv; hÞ ¼

1 ða Þ j hjrvj2 ; 2

a ¼ s; f :

ð40Þ

It is worth noting that the above equation is in agreement with Eq. (19). Using the free energy defined by Eqs. (31), (39), and (40), the stress defined in Eq. (20) for the a component has the following expression

rðaÞ ¼

1

J

FðaÞ  ðaÞ

þ

1

Z t h 0

Z

J ðaÞ

 i

ðaÞ CðaÞ ðvÞ : DðaÞ EðaÞ ðsÞ  K ðaÞ ðvÞDðaÞ trEðpaÞ ðvÞ þ trEh ðhÞ I ds  FðaÞT

 ðaÞ K e/ ðvÞDðaÞ /ðsÞds BðaÞ ðtÞ  qðaÞ jðaÞ hðrv  rvÞ;

t

0

a ¼ s; f ;

ð41Þ

where and in what follows, the dependency of the order parameter v on the time variable s is dropped. The entropy for the a component is as follows:

gðaÞ ¼

Z

Z t  ðaÞ    dEðaÞ ðhÞ 1 dE ðhÞ ðaÞ ða Þ ðaÞ ðaÞ ðaÞ h tr tr h þ  ½cðaÞ ðvÞv ðaÞ ðhÞ K ð v ÞD tr E ð s Þ d s K ð v ÞD /ð s Þd s e/ ðaÞ ð a Þ dh dh q0 0 q0 0 h i 1 0 0 ðaÞ ðaÞ þ b ðvÞuðaÞ ðhÞ þ W ðaÞ g ðaÞ ðvÞ  h cðaÞ ðvÞv ðaÞ ðhÞ þ b ðvÞuðaÞ ðhÞ  jðaÞ jrvj2 ; a ¼ s; f ; 2 1

t

ð42Þ

ðaÞ

where q0 is the density of the a component for the reference state. Using Eqs. (27)1, (31) and (39), the evolution equation for the order parameter of the porous medium is obtained

" # Z t dEpðrÞ ðvÞ dEpðrÞ ðvÞ 1 ðrÞ ðrÞ xv ¼ Mv ðrÞ K ðvÞD trE ðsÞds  tr þ M v ðrÞ K e/ ðvÞD /ðsÞds  tr dv dv hJ hJ 0 0 h 0 i 0 0 ðr Þ  Mv qðrÞ cðrÞ ðvÞv ðrÞ ðhÞ þ b ðvÞuðrÞ ðhÞ þ W ðrÞ g ðrÞ ðvÞ þ M v r  ðqðrÞ jðrÞ rvÞ: Z

1

ðf Þ

t

ðrÞ

ðrÞ



ðrÞ

ð43Þ

Similarly, the evolution equation for the porosity of the porous medium can be determined by Eqs. (27)2, and (31) as follows:

(

xðsÞ / ¼ M /

Z

1

J ðrÞ

t

0

) Z t h i 1 ðr Þ ðr Þ ðrÞ K e/ ðvÞDðrÞ tr EðrÞ ðsÞ  EpðrÞ ðvÞ  Eh ðhÞ ds þ ðrÞ K // ðvÞDðrÞ /ðsÞds : J 0

ð44Þ

Likewise, using Eqs. (27)4, (31) and (39), the drag force between the solid skeleton and pore fluid has the expression

( ðf Þ

T

¼ bv v ( 

ðf ;sÞ

1

J ðsÞ

 kv t rh 

Z

Z

1

J ðf Þ

t 0

ðf Þ K e/ ð

t ðsÞ

ðsÞ

ðsÞ

K ðvÞD trE ðsÞds

0 0

ðsÞ0

ðf Þ

h

ðf Þ

vÞD tr E ðsÞ 



)

Epðf Þ ð

dEðsÞ p ðvÞ  tr rv  dv

vÞ 

(

1

J ðsÞ

i

ðf Þ Eh ðhÞ

Z 0

t

) ds r/ 

ðsÞ K e/ ð

ðsÞ

"

1

J ðf Þ 

vÞD /ðsÞds

Z 0

t

# ðf Þ K // ð

ðf Þ

vÞD /ðsÞds r/

) dEðsÞ p ðvÞ  tr rv dv

0

þ qðsÞ h½cðsÞ ðvÞv ðsÞ ðhÞ þ b ðvÞuðsÞ ðhÞ þ W ðsÞ g ðsÞ ðvÞrv  r  ðqðsÞ jðsÞ hrvÞrv:

ð45Þ

Eq. (43) determines the evolution for the order parameter of the porous medium. It implies that the evolution of the order parameter is coupled with the deformation of the porous medium. Eqs. (41)–(45) show that the expressions for the variables with the memory effect associated with phase transition are different from the conventional viscoelasticity case (Christensen, 1982). For viscoelasticity case, stresses are represented by the convolution form of the Stieltjes integral, while the Stieltjes integrals appeared in the expressions for the variables here are in summation form. The convolution type memory effect of the traditional viscoelastic material is associated with the relaxation of the material rather than the ‘‘aging’’ of the material. The constitutive relation of a traditional viscoelastic material, therefore, exhibits the time translation invariance characteristic (Christensen, 1982). However, for the freezing porous medium free of relaxation under consideration, as phase transition will result in modulus variations with time, it thus displays an ‘‘aging’’-like property. Consequently, as noted above, its constitutive relation is not in convolution form and thus it is time-translation variant. 5. Some discussions about the boundary conditions and the determination of some parameters in the model In this section, we shall discuss briefly the boundary conditions for the developed phase field model and the determination of some parameters involved in the model. As the focus of this contribution is to develop a phase field model for the freezing porous medium, thus, parameters need involved procedures for determination will not be addressed here.

778

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

When applying our phase field model to a specific phase transition problem of the saturated porous medium, concrete boundary conditions should be prescribed a priori. For the boundary conditions associated with temperature, as the temperatures of the solid skeleton and pore fluid are assumed to be equal, two energy balance equations for the solid skeleton and pore fluid are not independent and thus they can be combined to form an independent one. As a result, it is straightforward to specify the two kinds of boundary condition (temperature or heat flux) associated with temperature at the boundary of a concerned domain. There are also two kinds of boundary condition associated with displacement, i.e., the displacement boundary condition and the traction boundary condition. The displacement boundary condition for the saturated porous medium can be imposed without any difficulty by prescribing the solid skeleton and pore fluid corresponding boundary values along the boundary. However, special care has to be taken when specifying the traction boundary condition for the saturated porous medium, as the solid skeleton and pore fluid have independent stresses and tractions, while only the total traction is given along the boundary. As a result, the manner for specifying the traction boundary condition is not unique (Rajagopal & Tao, 1995). Some kinds of traction boundary condition have now been proposed to avoid the above-mentioned non-unique problem, for example, the saturation boundary condition (Shi, Rajagopal, & Wineman, 1981), the traction splitting boundary condition (Rajagopal & Tao, 1995) and the continuity of chemical potential condition (Baek & Srinivasa, 2004), respectively. In terms of the research by Prasad and Rajagopal (2006), for simplicity, we choose to use the traction splitting boundary condition for the saturated porous medium. According to the traction splitting boundary condition (Rajagopal & Tao, 1995), the traction boundary condition for the saturated porous medium is as follows:

tðsÞ ¼ ð1  /Þt;

tðf Þ ¼ /t;

ð46Þ (s)

(f)

where t is the total traction on the boundary of the saturated porous medium; t and t are the tractions acting on the solid skeleton and pore fluid, respectively. For the moduli k(a), l(a) of the a component, it is assumed that their values for the complete unfrozen state (v = 0) and ðaÞ ðaÞ ðaÞ ðaÞ frozen state (v = 1) can be measured via experiments, which are represented by k0 ; l0 and k1 ; l1 , respectively. The moduli for the freezing a component thus have the following expressions ðaÞ

ðaÞ

ðaÞ

ðaÞ

kðaÞ ðvÞ ¼ k0 þ pðvÞðk1  k0 Þ ¼ k0 þ pðvÞDkðaÞ ;

lðaÞ ðvÞ ¼ lð0aÞ þ pðvÞðlð1aÞ  lð0aÞ Þ þ pðvÞDlðaÞ ; a ¼ s; f :

ð47Þ

As mentioned above, as the pore fluid is also treated as a ‘‘solid’’ in our model, thus, at the complete unfrozen state, the shear modulus of the pore fluid should vanish. However, for the stability of algorithm, it is better for the shear modulus to take on a ðaÞ comparative small value at the state. Other moduli, such as the bulk modulus K(a), the moduli K e/ ðvÞ and K a// (v) can also be ðaÞ determined in an analogous manner. Note that the method for the measurement of the moduli K e/ ðvÞ and K a// ðvÞ is similar to that for normal porous media (Lu & Hanyga, 2005; Lu, Hanyga, & Jeng, 2007a). The permeability of the porous medium (kv) is also dependent on the degree of the phase transition of the porous medium. We suppose that at the complete unfrozen state, the permeability of the porous medium takes on the value of the conventional permeability kv0, while at the complete frozen state, the porous medium is completely impermeable and its permeability, thus, vanishes. As a result, the permeability for the freezing porous medium and bv in Eq. (45) have the following representations

kv ðvÞ ¼ ½1  pðvÞkv 0 ; bv ðvÞ ¼

1 : ½1  pðvÞkv 0

ð48Þ

Note that when x ? 1, bv(x) given by Eq. (48) will approach infinity. In numerical calculation, this difficulty can be avoided by assuming comparatively large bv(x) when x ? 1 Likewise, for the coefficient kvt in Eq. (45) which accounting for the coupling between the relative motion between the solid skeleton and pore fluid and heat flux, the similar expression can be obtained

kv t ðvÞ ¼

1 kv t0 ½1  pðvÞ

ð49Þ

in which kvt0 is the value of the coefficient kvt at the complete unfrozen state. For the coefficient of thermal conduction, if its values for the solid skeleton and pore fluid at the complete unfrozen and ðaÞ ðaÞ frozen states are kh0 ; kh1 ða ¼ s; f Þ, respectively, then, in view of Eq. (16), the coefficient of thermal conduction for the freezing porous medium is given by

kh ðvÞ ¼

X a

ðaÞ kh0

þ pðvÞ

X a

ðaÞ kh1



X

!

ðaÞ kh0

:

ð50Þ

a

The parameters in Eqs. (38) and (39) can be obtained using the material constants of the concerned porous medium. The parameters w(a), k(a), Mx in Eqs. (39), (40), (43) and their relations with the thickness, excess free energy of the transition layer and the material parameters of the freezing material can be determined by the procedures outlined by Caginalp (1989) and Boettinger et al. (2002). The relaxation parameter for the porosity M / in Eq. (44) can be measured for the complete unfrozen state and frozen state first and its value for an arbitrary phase state can then be obtained by interpolation as shown in Eq. (47). The approach for the measurement of M / for the common porous medium is discussed by Lu, Hanyga and Jeng (2007b).

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

779

6. Conclusions In this paper, a phase field model for the freezing saturated porous medium is developed in the context of continuum physics in a thermodynamically consistent manner. The freezing saturated soil with a transition layer, i.e., the freezing fringe, between the frozen and unfrozen regions, is characterized by the modern phase field theory for the first time. The underlying fundamental concepts for the proposed model are as follows. Firstly, the film fluid surrounding the solid grains is regarded as part of the solid skeleton, and the pore fluid thus only includes free flow pore fluid. This definition for the solid skeleton and pore fluid is critical for the reasonable description of the freezing porous medium. Besides, introduction of the order parameter for the solid skeleton is also fundamental; it facilitates the description of the drag force acting on the pore fluid by the freezing film water in the freezing fringe. Moreover, introducing summation Stieltjes type integral to characterize the variables with the memory effect associated with phase transition is also a key element of our model. The evolution equation for the order parameter shows that the phase transition of the porous medium is completely coupled with the deformation of the porous medium. Also, the stress of the porous medium is related to the phase state of the porous medium. The expression for the drag force contains two terms with the gradient of the order parameter, which implies that the spatial variation of the film fluid thickness also contributes to the drag force. It is found that the heat flux vector of the freezing porous medium also involves the gradient of the order parameter. This characteristic for the heat flux vector is absent in the existing phase field models. Also, we find that the stresses, entropies, drag force and the time change rates of the order parameter and porosity all depend on the histories of their variables via summation type Stieltjes integral, which implies that the memory effect associated with the variation of the moduli during a phase transition is an ‘‘aging’’-like property other than viscoelasticity. This kind of memory effect displayed in the phase change materials has not been reported and investigated so far. Naturally, the developed model can be used in the predication and calculation of the heave and thawing of the freezing saturated soil. However, it is noted that although the developed model is mainly aimed at the freezing saturated soil, it remains quite general. Thus, the proposed theory can also be used in the analysis of the freezing saturated concrete and various saturated biological tissues. Based on the theoretical model established in this study, the corresponding finite element model can also be developed for numerical simulation of the freezing saturated porous medium. Also, the proposed model can be extended to the porous medium with solution pore fluids, which has implications in the salinization and salt expansion of the freezing saturated soil. Acknowledgements This research is supported by an open funding from State Key Laboratory for GeoMechanics and Deep Underground Engineering of China (No. SKLGDUE08012X). Also, the financial support received from NSFC is highly acknowledged (No. 51078171). Moreover, discussion with Polish Prof. Andrzej Hanyga about the research of this paper is highly appreciated by the first author. The constructive comments from the reviewers are highly acknowledged by the authors. References Ahmadi, G. (1982). A generalized continuum theory for granular materials. International Journal of Non-linear Mechanics, 17, 21–33. Akagawa, S. (1988). Experiment al study of frozen fringe characteristics. Cold Regions Science and Technology, 15, 209–223. Allen, S. M., & Cahn, J. W. (1979). A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metallurgica et Materialia, 27, 1085–1095. Anderson, D. M., McFadden, G. B., & Wheeler, A. A. (2000). A phase-field model of solidification with convection. Physica D, 135, 175–194. Atkin, R. J., & Craine, R. E. (1976). Continuum theories of mixtures: basic theory and historical development. Quarterly Journal of Mechanics and Applied Mathematics, 29, 209–244. Baek, S., & Srinivasa, A. R. (2004). Diffusion of a fluid through an elastic solid undergoing large deformation. International Journal of Non-linear Mechanics, 39, 201–218. Bedford, A., & Drumheller, D. S. (1983). Theories of immiscible and structured mixtures. International Journal of Engineering Science, 21, 863–960. Boer, R. D., & Kowalski, S. J. (1995). Thermodynamics of fluid-saturated porous media with a phase change. Acta Mechanica, 109, 167–189. Boettinger, W. J., Warren, J. A., & Beckermann, C. (2002). Phase-field simulation of solidification. Annual Review of Materials Research, 32, 163–194. Bowen, R. M. (1980). Incompressible porous-media models by use of the theory of mixtures. International Journal of Engineering Science, 18, 1129–1148. Bowen, R. M. (1982). Compressible porous-media models by use of the theory of mixtures. International Journal of Engineering Science, 20, 697–735. Caginalp, G. (1989). Stefan and Hele-Shaw type models as asymptotic limits of the phase-field equations. Physical Review A, 39, 5887–5896. Cahn, J. W., & Hilliard, J. E. (1958). Free energy of a nonuniform system. I. Interfacial free energy. Journal of Chemical Physics, 28, 258–267. Chacon, R. V. S., & Rivlin, R. S. (1964). Representation theorems in the mechanics of materials with memory. Zeitschrift Fur Angewandte Mathematik Und Physik, 15, 444–447. Chapman, S., & Cowling, T. G. (1970). The mathematical theory of non-uniform gases (3rd ed.). New York: Cambridge University Press. Chen, L. Q. (2002). Phase-field models for microstructure evolution. Annual Review of Materials Science, 32, 113–140. Christensen, R. M. (1982). Theory of viscoelasticity: An introduction. New York: Academic Press, Inc. Coleman, B. D., & Noll, W. (1963). The thermodynamics of elastic materials with heat conduction and viscosity. Archive for Rational Mechanics and Analysis, 13, 167–178. Coussy, O. (2005). Poromechanics of freezing materials. Journal of the Mechanics and Physics of Solids, 53, 1689–1718. de Groot, S. R., & Mazur, P. (1984). Non-equilibrium thermodynamics. New York: Dover. Drumheller, D., & Bedford, A. (1980). A thermomechanical theory for reacting immiscible mixtures. Archive for Rational Mechanics and Analysis, 73, 257–284. Duquennoi, C., & Fremond, M. (1989). Modelling of thermal soil behavior. In VTT symposium 94 (Vol. 2, pp. 895–915). Eringen, A. C. (1994). A continuum theory of swelling porous elastic soils. International Journal of Engineering Science, 32, 1337–1349. Fabrizio, M., Giorgib, C., & Morro, A. (2006). A thermodynamic approach to non-isothermal phase-field evolution in continuum physics. Physica D, 214, 144–156.

780

J.-F. Lu et al. / International Journal of Engineering Science 49 (2011) 768–780

Fix, G. J. (1983). Phase field methods for free boundary problems. In A. Fasano & M. Primicerio (Eds.), Free boundary problems: Theory and applications (pp. 580–589). London: Pitman. Fosdick, R. L., & Rajagopal, K. R. (1980). Thermodynamics and stability of fluids of third grade. Proceedings of the Royal Society of London Series A, 339, 351–377. Fremond, M., & Mikkola, M. (1991). Thermomechanical modelling of freezing soil. In: Balkema (Ed.), Proceedings of the sixth international symposium on ground freezing, Rotterdam (pp. 17–24). Fried, E., & Gurtin, M. E. (1993). Continuum theory of thermally induced phase transitions based on an order parameter. Physica D, 68, 326–343. Fried, E., & Gurtin, M. E. (1996). A phase-field theory for solidification based on a general anisotropic sharp-interface theory with interfacial energy and entropy. Physica D, 91, 143–181. Gilpin, R. R. (1980). A model for the prediction of ice lensing and frost heave in soils. Water Resources Research, 16, 918–930. Golden, J. M. (2008). Phase transition in materials with thermal memory. Physica D, 237, 2499–2516. Gurtin, M. E., & Pipkin, A. C. (1968). A general theory of heat conduction with finite wave speeds. Archive for Rational Mechanics and Analysis, 31, 113–126. Harlan, R. L. (1973). Analysis of coupled heat–fluid transport in partially frozen soil. Water Resources Research, 9, 1314–1323. Hassanizadeh, S. M., & Gray, W. G. (1990). Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Advances in Water Resources, 13, 169–186. Jeong, J. H., Goldenfeld, N., & Dantzig, J. A. (2001). Phase field model for three-dimensional dendritic growth with fluid flow. Physical Review E, 64, 041602. Konrad, J. M., & Morgenstern, N. R. (1980). A mechanistic theory of ice lens formation in fine grained soils. Canadian Geotechnical Journal, 17, 476–486. Langer, J. S. (1986). Models of pattern formation in first-order phase transitions. In G. Grinstein & G. Mazenko (Eds.), Directions in condensed matter physics. Philadelphia: World Scientific. Lu, J.-F., & Hanyga, A. (2005). Linear dynamic model for porous media saturated by two immiscible fluids. International Journal of Solids and Structure, 42, 2689–2709. Lu, J.-F., Hanyga, A., & Jeng, D.-S. (2007a). A mixture-theory-based dynamic model for a porous medium saturated by two immiscible fluids. Journal of Applied Geophysics, 62, 89–106. Lu, J.-F., Hanyga, A., & Jeng, D.-S. (2007b). A linear dynamic model for a saturated porous medium. Transport in Porous Media, 68, 321–340. Marchand, J., Pleau, R., & Gagne, R. (1995). Deterioration of concrete due to freezing and thawing. In S. Mindess & J. Skalny (Eds.), Material science of concrete (pp. 283–354). Westerville, Ohio, USA: American Ceramic Society. Miller, R. D. (1972). Freezing and heaving of saturated and unsaturated soils. Highway Research Record, 393, 1–11. Morland, L. W. (1972). A simple constitutive theory for a fluid-saturated porous solid. Journal of Geophysical Research, 77, 890–900. O’Neill, K., & Miller, R. D. (1985). Exploration of a rigid ice model of frost heave. Water Resources Research, 21, 281–296. Passman, S. L., Nunziato, E. W., & Walsh, E. K. (1984). A theory of multiphase mixtures. In C. A. Truesdell (Ed.), Rational thermodynamics (3rd ed., pp. 286–325). Hopkins University Press (Appendix in Truesdell, C. A., Rational thermodynamics). Penrose, O. (1990). Thermodynamically consistent models of phase field type for the kinetics of phase transitions. Physica D, 43, 44–62. Prasad, S. C., & Rajagopal, K. R. (2006). On the diffusion of fluids through solids undergoing large deformations. Mathematics and Mechanics of Solids, 11, 291–305. Provatas, N., & Elder, K. (2010). Phase-field methods in material science and engineering. Wiley-VCH Verlag. Rabin, Y., & Steif, P. S. (2000). Thermal stress modeling in cryosurgery. International Journal of Solids and Structures, 37, 2363–2375. Radd, F. J., & Oertle, D. H. (1973). Experimental pressure studies of frost heave mechanism and the growth fusion behavior. In 2nd international conference on permafrost (pp. 377–384). Washington, DC: National Academy Press. Rajagopal, K. R., & Tao, L. (1995). Mechanics of mixtures. Singapore: World Scientific. Rudin, W. (1973). Functional analysis. New York: McGraw-Hill, Inc. Sheppard, M. T., Kay, B. D., & Loch, J. P. G. (1978). Development and testing of a computer model for heat and mass flow in freezing soils. In Proceedings of 3rd international conference on permafrost, Edmonton, Alberta, Canada (pp. 76–81). Shi, J. J., Rajagopal, K. R., & Wineman, A. S. (1981). Application of the theory of interacting continua to the diffusion of a fluid through a non-linear elastic media. International Journal of Engineering Science, 19, 871–889. Thigpen, L., & Berryman, J. G. (1986). Mechanics of porous elastic materials containing multiphase fluid. International Journal of Engineering Science, 23, 1203–1214. Truesdell, C., & Noll, W. (1965). The non-linear field theories of mechanics. In S. Flügge (Ed.), Encyclopedia of physics. Vol. III/3. Berlin: Springer (Sect. 19). Truesdell, C., & Toupin, R. (1960). The classical field theories. In S. Flügge (Ed.). Handbook der physik (Vol. III/1). Berlin: Springer. Tsytovich, N. A. (1973). Mechanics of frozen soils. Moskow: Vysshaya shkola. Van der Waals, J. D. (1893). Verhandel. Konik. Akad. Weten. Amsterdam (Sec. 1) Vol. 1, No. 8. Wang, S.-L., Sekerka, R. F., Wheeler, A. A., Murray, B. T., Coriell, S. R., Braun, R. J., & McFadden, G. B. (1993). Thermodynamically-consistent phase-field models for solidification. Physica D, 69, 189–200. Wei, C. F., & Muraleetharan, K. K. (2002). A continuum theory of porous media saturated by multiple immiscible fluids: I. Linear poroelasticity. International Journal of Engineering Science, 40, 1807–1833.