Propagation of cohesive crack in unsaturated porous media by use of XFEM method

Various damage models have been presented to model the cracks. These models don't simulate the discontinuity in displacement field; just see its effect on an upper level on stress field through damage parameter. Damage parameter is used to modify other fields like pressure. This paper tries to model crack itself through fracture mechanics. The opening of the crack is modeled and is connected to stress field by use of cohesive crack model that relates plasticity of crack tip to its opening. When crack opens, there exists discharge of water and gas through crack that are calculated by use of multiscale methods and incorporated in weak form of governing equations.Due to complexity of proposed model, it can't be solved analytically and extended finite element method (XFEM) is used to solve it. Xfem method which is based on partition of unity property of FEM adds some terms to conventional interpolation functions to model discontinuity of fields adjacent to the crack. This method simulates crack in ordinary mesh of FEM and doesn't need remeshing around the crack when it propagates.


Introduction
In the present study the propagation of crack in unsaturated porous media is investigated.To describe the behavior of unsaturated porous media Gatmiri's equations [1][2][3] are used.The differential forms of these equations are written based on the two independent variables of net stress and suction.
Cohesive crack model [4][5][6] is used for mechanical interaction of porous media and crack.In this model, the resisting force against opening of the crack is determined based on an opening-force relation.Magnitude of this force varies from tensile strength of material for no opening to zero for critical opening ( c w ).Critical opening is the opening that cohesive force is zero after that.
It is assumed that water pressure has weak discontinuity in vicinity of the crack.Thus pressure gradient would be discontinues in neighborhood of the crack.This discontinuity leads to water flow through the opening.Deduction of this flow through the opening needs rewriting of mass equilibrium of water around the crack.Micro and macro coupling is obtained by adding the derived term for flow through the crack to the integral form of mass balance equation of water.
XFEM [7][8][9] method is used to descretize the resulted equations in space.This method which is based on partition of unity property of finite element shape functions adds some enrichment to ordinary shape functions and doesn't need any remeshing after crack propagation and allows modeling the discontinuity in a regular FEM mesh.Finite difference method is used for temporal descretization.The resulted equations are nonlinear and Newton-Raphson method is used to solve them.

Water Phase Transfer Equation
Considering the generalized Darcy equation for movement of Water phase in porous media, the equation of water flow is defined as: Where U is the liquid velocity and q w is the liquid flow vector.ρ w and K w are defined as water density and water permeability.s is the suction (P g -P w ) and Z is the gravity potential.
w K is total water permeability that is written in the following form: In above equation δ is identity tensor of order two, e the void ratio, w  a model parameter and is the relative permeability that shows fluid properties.It may be computed based on a state surface [10,11] or more simply a retention curve [12].The following formula is used for relative permeability: υ is the viscosity of fluid which is considered to be equal to 72 9.94 10 . .


 N S m .

Water Mass Conservation Equation:
Considering a representative volume element (RVE) of porous media, time derivative of included mass of water becomes equal to divergence of mass flow of water.Mathematical form of this equilibrium may be expressed as: Where n, S w , ρ w and u are defined as the porosity, saturation degree, water density and water velocity respectively.
The left side of mass balance equation can be expanded as: Solid particles are incompressible so the derivative of porosity is equal to the volumetric strain that is defined as sum of normal strains: The density of water is assumed to be function of water pressure, so its time derivative becomes dependable on pressure time derivative as: . ..
For calculation of saturation degree, the Van Genuchten [12] model that relates saturation and suction is considered here: S wr is the residual saturation degree that saturation degree can't reduce to a value less than it.s is the suction and α and n are the constants of the model.
The derivative form of the water saturation degree is written as: By substituting ( 6), ( 7), ( 8) and ( 10) for ( 5) water mass conservation equation on the basis of principle variables is concluded:

Solid Skeleton Behavior:
Considering balance of momentum equation for unsaturated media as: is the net stress and b is the volumetric force.
Constitutive law of solid skeleton in unsaturated porous medium can be presented on the basis of net stress and suction that are two independent variables:

Hydromechanical Coupling
The force which is exerted on the body of the crack is the combination of two factors: the first is the plasticity of cracked solid medium and the second is suction of unsaturated media.This force might be demonstrated as: .

Integration and Descritization of Governing Equations
To Obtain the integral form of governing equations, u * (allowable virtual displacement) and   * (allowable virtual liquid pressure) are multiplied in equations ( 11) and ( 12) respectively.The gas pressure is assumed to be equal to the atmospheric pressure.After using the Green theorem and some manipulation, final form of integrated equations are resulted: .( .(P P ) .).P ( ).P The The flow through the crack is involved in integral form of governing equations either.
Spatial discretization of governing equations is done through XFEM.XFEM method presumes the interpolation function of unknown variables consists of standard and enriched parts, so the unknown variables can be given as: After applying the spatial and temporal discretization, the matrix form of governing equations becomes: Since the existence of the terms representing internal force and discharge of water through crack, the system of equations become nonlinear and it becomes necessary to use a repetitive procedure like Newton-Raphson to solve the equations.

Example
In this section, the ability of model to simulate crack propagation in unsaturated porous media is investigated.A square plate with dimension of 250×250 mm and initial notch length of 50 mm along its centerline is considered.The suction distribution is shown in Figure 4. Suction increases around the crack body and because of the gradient of suction which is shown in Figure 5 water is absorbed to the crack.For having better assessment of the model performance, a comprehensive analyze is done.The influence of initial suction degree on water discharge at a specified time is shown in Figure 6 .Water exchange increases with decrease in initial saturation degree because when initial saturation degree decreases, suction increases and crack has more inclination to absorb water.Two mechanisms are at work in discharge through crack: The first is the porosity change and the second is plane Poiseuille flow which relates the cube of opening to discharge of water through crack.Here just the first mechanism is taken into account, therefore discharge of water increases in proximity of the crack.The influence of initial saturation degree on crack opening is shown in Figure 7.As it is shown in Figure 7, the opening of crack increases when the initial saturation degree decreases.Due to higher suction in lower initial saturation degree, the bearing capacity of the plate increases and it leads to higher stresses in the plate.These higher stresses cause the crack to propagate faster.Consequently at a specific time, length of crack and its opening is more for lower initial saturation degree.

Conclusion
It was shown the XFEM method has the potential to model crack in unsaturated porous media in a regular FEM mesh.A formulation is outlined to model unsaturated porous media and the terms related to crack were added to weak form of this formulation.A problem was solved and it was shown how the suction concentrates around the crack and not only helps to close the opening but also absorbs the water to the crack.
 and   * represent the average and difference of solid velocity at the two opposite sides of the crack, respectively.2h shows opening of the crack and subscript d indicates that the parameter is specific to the crack.Mechanical and hydraulic coupling terms are connected through d n that is porosity of damaged porous media and is defined as: w is opening of the crack.To calculate the coupling terms, two Gauss points are placed on every segment of crack path which is bounded in an element (Figure2 ).

Figure2.
Figure2.Finite element mesh.Blue circles show enriched nodes and crosses show Gauss points on the crack path.
suction forces exerted on the opened crack body (     u means opening of the crack).Rewriting the first term of the right hand side of equation (17) (b) as:

w
Pt are the enriched ones.The shape functions are substituted for equilibrium and mass conservation equation.Considering implicit finite difference method, temporal discretization can be done.

.
The upper and lower edges of plate are stretched with fixed vertical velocity ̇= 2.35 mm s Figure 3. shows the geometry of crack and the properties of media are mentioned in Table1.The edges of plate are impermeable to water and the plate is considered saturated, initially.

Figure 3 .
Figure 3.The geometry of the problem

Figure 4 .Figure 5 .Figure 7 .Figure 6 .
Figure 4. Suction distribution in plate is the intrinsic permeability,   [13,14]cohesive force versus openingDischarge of water through the crack is calculated similar to[13,14].It is supposed that the length of the crack is much more greater than its height and the changes of w S , and d n is negligible on the height of the crack.With these assumptions, flow of water through the crack is calculated as: d t is the cohesive force which is calculated based on the openingforce relation shown on Figure2, w S ()  gw PP is the suction force that is assumed to be similar on opposite faces of the crack, d Γ n is normal to the crack and d Γ shows the body of the crack.[[u]] td Wc Figure1.w  , w P , g P )