COMSOL In A New Tensorial Formulation Of Non-Isothermal .

3y ago
43 Views
2 Downloads
457.65 KB
7 Pages
Last View : 16d ago
Last Download : 2m ago
Upload by : Grant Gall
Transcription

Excerpt from the Proceedings of the COMSOL Conference 2009 BostonCOMSOL in a New Tensorial Formulation of Non-IsothermalPoroelasticityMario-César Suárez A.*, 1 and Fernando Samaniego V.21Faculty of Sciences Michoacán University, 2Faculty of Engineering National University of Mexico* Edificio B Cd. Universitaria, 58060 Morelia, Mich., Mexico, email: mcsa50@gmail.comAbstract: The presence of a moving fluid in aporous rock modifies its mechanical response.Poroelasticity explains how the fluid inside thepores bears a portion of the total load supportedby the rock. The remaining part of the load issupported by the elastic skeleton, which containsa laminar fluid coupled to the framework byequilibrium and continuity conditions. This workintroduces an original tensorial formulation ofBiot’s theory and the experimental thermoporoelastic parameters that support the theory.By defining a 4-dimensional total stress tensorand three basic poroelastic coefficients, wededuce a system of differential equationscoupling the bulk rock and the fluid. Theinclusion of the fourth dimension allowsextending the theory of solid linear elasticity tothermoporoelastic rocks, taking into account theeffect of both the fluid phase and thetemperature. Introducing three volumetricthermal dilation coefficients, one for the fluidand two for the skeleton, a complete set ofparameters for geothermal poroelastic rocks areobtained.Keywords: Poroelasticity, thermoporoelasticity,Biot’s theory, finite elements.external forces. In this paper, we present anoriginal four-dimensional tensorial formulationof linear thermo-poroelasticity theory. Thisformulation makes more comprehensible thelinear Biot’s theory, rendering the resultingequations more convenient to be solved using theFinite Element Method. To illustrate practicalaspects of our model some classic applicationsare outlined and solved using the Earth ScienceModule of COMSOL-Multiphysics.1.1 Experimental backgroundIn classic elastic solids only the two Lamémoduli, (l, G) or Young’s elastic coefficient andPoisson’s ratio (E, n), are sufficient to describethe relations between strains and stresses. Inporoelasticity, we need five poroelastic modulifor the same relationships (Bundschuh andSuárez, 2009), but only three of these parametersare independent. The Biot’s field variables for anisotropic porous rock are the stress s acting inthe rock, the bulk volumetric strain εB, the porepressure pf and the variation of fluid masscontent ζ. The linear relations among thesevariables are the experimental foundations ofBiot’s poroelastic theory (Biot & Willis, 1957;Wang, 2000):1. IntroductionSeveral factors affect the geomechanicalbehavior of porous crustal rocks containingfluids: porosity, pressure, and temperature,characteristics of the fluids, fissures, and faults.Rocks in underground systems (aquifers,geothermal and hydrocarbon reservoirs) areporous, compressible, and elastic. The presenceof a moving fluid in the porous rock modifies itsmechanical response. Its elasticity is evidencedby the compression that results from the declineof the fluid pressure, which can shorten the porevolume. This reduction of the pore volume canbe the principal source of fluid released fromstorage. A rock mechanics model is a group ofequations capable of predicting the porousmedium deformation under different internal andεB σKB pfH ε C B B 1 ζ H, ζ σH pfR(1)H σ R 1 p f 1Where KB, H, and R are poroelastic coefficientsthat are experimentally measured as follows(Wang, 2000):εB Δ VB1 Δε , CB B , K B VBCB Δσ p f1 ΔεB H Δp f 1 Δζ Δζ , σ Δσ p f R Δp f σ(2)

The following figure (1) illustrates all theparts forming a poroelastic medium.Biot (1941) and (Biot & Willis, 1957)introduced three additional parameters, b, M andC, that are fundamental for the tensorialformulation herein presented. 1/M is called theconstrained specific storage, which is equal tothe change of ζ when pf changes measured atconstant strain. Both parameters M and C areexpressed in terms of the three fundamental onesdefined in equation (2):1 Δζ M Δ p f 1 KB R H2 ε B (4)KRH2; C BMM 2H KB RHFigure 1. Skeleton of sandstone showing its pores andsolid grains (3 3 3 mm3). (Piri, 2003).Here VB is the bulk volume, consisting of therock skeleton formed by the union of the volumeof the pores VF and the volume of the solidmatrix VS (Fig. 1). The control volume is DVB.The drained coefficients KB and CB are the bulkmodulus and the bulk compressibility of therock, respectively; 1/H is a poroelastic expansioncoefficient, which describes how much DVBchanges when pf changes while keeping theapplied stress s constant; 1/H also measures thechanges of ζ when s changes and pf remainsconstant. Finally 1/R is an unconstrained specificstorage coefficient, which represents the changesof ζ when pf changes. Inverting the matrixequation (1) and replacing the value of s in ζ weobtain:Kσ K Bε B B p f H(3)KB 1 KB ζ εB 2 pfH R H The sign conventions are stress s 0 in tensionand s 0 in compression; the volumetric strainεB 0 in expansion and εB 0 in contraction; thefluid content ζ 0 if fluid is added to the controlvolume DVB and ζ 0 if fluid is extracted fromDVB; the pore pressure pf 0 if it is larger thanthe atmospheric pressure.Let CS 1/KS be the compressibility of the solidmatrix. The Biot-Willis coefficient b is definedas the change of confining pressure pk withrespect to the fluid pressure change when thetotal volumetric strain remains constant: pkb pf KC KB 1 B KS MH ε(5)The coefficient C represents the coupling ofdeformations between the solid grains and thefluid. The coefficient M is the inverse of theconstrained specific storage, measured atconstant strain (Wang, 2000); this parametercharacterizes the elastic properties of the fluidbecause it measures how the fluid pressurechanges when ζ changes.These three parameters b, M and C are at thecore of the poroelastic partial differentialequations we introduce herein (Bundschuh andSuárez, 2009).2. Isothermal Poroelasticity ModelLet us and uf be the displacements of thesolid and fluid particles; let u uf – us be thedisplacement of the fluid phase relative to thesolid matrix respectively. Let εs, εf, js, j, Vs andVf be the volumetric dilatations, porosities andvolumes of each phase; – εV is the volumetricdeformation of the fluid phase relative to thesolid phase. The mathematical expressions ofthese variables are:

G G ΔV fG GΔVS ε s uS ; ε f u fVSVfG GGε V ε S ε f ; u u f uS ε V (6)This tensorial equation becomes identical tothe Hookean solids equation, when the rock haszero porosity and b 0. From equations (8) and(9), we deduce that:G GG G u u y u zG uS u f u x x y z()Biot and Willis (1957) introduced the strainvariable ζ (u, t), defined in equation (3), todescribe the volumetric deformation of the fluidrelative to the deformation of the solid withhomogeneous porosity:GGGGζ ( u ,t ) ϕ ( us u f ) ϕ ε s ϕ ε f ϕ ε V (7)The function ζ represents the variation offluid content in the pore during a poroelasticdeformation. The total applied stresses in theporous rock are similar to the equations ofclassic elasticity. However, we need to couplethe effect of the fluid in the pores. The linearcomponents of the global stresses, deducedexperimentally by Biot, (Biot, 1941; Biot andWillis, 1957; Wang, 2000) are:σ ij λU ε B δ ij 2 G ε ij C ζ δ ij1 ui u j 2 x j xi 1, if i j , δ ij 0, if i j λ C b ; for i , j x, y, zε ij ζ The fluid pressure is deduced from equation (3):KB R H 2 ζε B 2KH KB R B H (9)We define a two-order tensor sT (sij) infour dimensions, which includes the bulk stresstensor sB acting in the porous rock and the fluidstress sF acting in the fluid inside the pores,positive in compression:(σT σ B σ F )σ ij λ U ε B C ζ δ ij 2 G ε ijσ f p f M ζ C ε B ; i , j x, y , zτ ij λ eB δ ij 2 Gε ij(12)Tensor tij is called the Terzaghi (1943)effective stress that acts only in the solid matrix;bpf is the pore-fluid pressure. Since there are noshear tensions in the fluid, the pore fluid pressureaffects only the normal tensions si (i x, y, z).The functions sij are the applied stresses actingin the porous rock saturated with fluid. The solidmatrix (tij) supports one portion of the totalapplied tensions in the rock and the fluid in thepores (bpf ) supports the other part. This is amaximum for soils, when b º 1 and is minimumfor rocks with very low porosity where b º 0.For this reason, b is called the effective stresscoefficient.Inverting the matrices of equations (8) and(9), we arrive to the following tensorial form ofthe poroelastic strains:ε ii pf (11)(8)Where:λUσ ij τ ij b p f δ ij(10)σM σ ii2GσMH pσ3νσ M f ; ε ij ij (13a)3H2GEpfR CσM KU pfσ xx σ yy σ zz3M KU C2(13b) KB ε B b M ζ2K B λ G ; KU K B b 2 M3(14)The coefficient KU is the undrained bulkmodulus, which is related to the previous definedcoefficients. Note that both tensorial equations(10) and (13) only need four basic poroelasticconstants. The presence of fluid in the pores addsan extra tension due to the hydrostatic pressure,which is identified with the pore pressure,because it is supposed that all the pores areinterconnected. This linear theory is appropriatefor isothermal, homogeneous, and isotropicporous rocks.

σ ij σ ij0 ( λ ε B b ( p p0 ) ) δ ij 3. ThermoPoroelasticity ModelThe equations of non-isothermal poroelasticprocesses are deduced using the Gibbs thermoporoelastic potential or available enthalpy perunit volume and the energy dissipation functionof the skeleton (Coussy, 1991).Analytic expressions are constructed in termsof the stresses, the porosity, the pore pressure,and the density of entropy per unit volume ofporous rock. As we did for the isothermalporoelasticity, we can write in a single fourdimensional tensor the thermoporoelasticequations relating stresses and strains. We havefor the pore pressure:p p 0 M ( ζ ζ 0 ) Cε B( Mϕ γϕ γ f) (T T 0 )(15)The volumetric thermal dilatation coefficientgB [1/K] measures the dilatation of the skeletonand gj [1/K] measures the dilatation of the pores:1 VB γB VB T pk 1 K 1 Vϕ 1 ϕ γϕ ϕ T p fVϕ T pf(16) 1 K (17)The fluid bulk modulus Kf and the thermalexpansivity of the fluid gf [1/K] are defined asfollows:11 ρ f (18) Cf ρ f p TKfγf 1 V f 1 ρ f ρ f T pV f T pff(19)The term pk is the confining pressure.Expanding the corresponding functions of theGibbs potential and equating to zero the energydissipation we obtain the 4D thermoporoelasticequations, which include the thermal tensions inthe total stress tensor (Bundschuh and Suárez,2009): 2 G ε ij K B γ B (T T0 )(20)In this case, an initial reference temperature T0and an initial pore pressure p0 are necessarybecause both thermodynamic variables T and pare going to change in non-isothermal processesoccurring in porous rock. The fluid stress isdeduced in a similar way:σ f p f M (ζ ζ 0 ) C ε B( Mϕ γϕ γ f) (T T0 )(21)4. Dynamic Poroelastic EquationsThe formulation we introduced herein is veryconvenient to be solved using the Finite ElementMethod. The fundamental poroelastic differentialequation is the tensorial form of Newton’ssecond law in continuum porous rock dynamics:GJJJG GGG 2u JJJG Gdiv σ T F ρ 2 ; div σ T T σ T tGGGGwhere : σ T B ε T ; ε T u(22)The terms sT and εT are the equivalent vectorialform of tensorial equations (20) and CB is thematrix of poroelastic constants. While F is thebody force acting on the rock and the tensordifferential operator is given by: x T 0 000 y z0 y0 x0 z0 z0 x y x y z (23) ux G u y ε T (ε x ε y ε z ε xy ε xz ε yz er ) u z Where u (ux, uy, uz) is the displacement vectorof equation (6). Using the operator inequation (22), the dynamic poroelastic equationbecomes:( TGGG 2u B u F ρ 2 t)(24)

4.1 Solution of Thermoporoelastic Equations:The Finite Element MethodEquation (24) includes Biot’s poroelastic theory.It can be formulated and numerically solvedusing the Finite Element Method (FEM). Let Ωbe the bulk volume of the porous rock, and let Ω be its boundary, u is the set of admissibledisplacements in Eq. (22); fb is the volumetricforce and fs is the force acting on the surface Ω.After doing some algebra we arrive to a FEMfundamental equation for every element Ve in thediscretization:G e d e e G 2 d e t2G F e ; e 1, M VTBedrockstep(25)de is a vector containing the displacements of thenodes in each Ve. Equation (25) approximates thedisplacement u of the poroelastic rock. Fe is thevector of total nodal forces. Ke and Me are thestiffness and equivalent mass matrices for thefinite element Ve. The mathematical definitionsof both matrices are: e layers of the sequence are each 20 m thick. Thefirst and third layers are aquifers; the middlelayer is relatively impermeable to flow.Figure 2. Simplified geometry of the aquifer and theimpermeable bedrock in the basin. Initial state. B d V ; e 1, Me ; e ρ T dV(26)VeWhere is the matrix of shape functions thatinterpolate the displacements (Liu and Quek,2003). Matrix is called the strain poroelasticmatrix.5. Use of COMSOL MultiphysicsThis section contains two brief illustrations ofthe deformation of an aquifer (Leake & Hsieh,1997) and the form that a temperature changecan affect its poroelastic deformation. In the firstexample, we assume cold water at 20 C (1000kg/m3). After, we consider a temperature of250 C (50 bar, 800.4 kg/m3). Results are shownin figures (4) to (9). To simplify the discussionwe use the same model previously solved byCOMSOL-Multiphysics and described in theEarth Science Module (COMSOL AB, 2006):“Three sedimentary layers overlay impermeablebedrock in a basin where faulting creates abedrock step (BS) near the mountain front (Fig.2). The sediment stack totals 420 m at thedeepest point of the basin (x 0 m) but thins to120 m above the step (x 4000 m). The top twoFigure 3. The mesh of the basin with 2967 elements.As given by the problem statement, the materialshere are homogeneous and isotropic within alayer. The flow field is initially at steady state,but pumping from the lower aquifer reduceshydraulic head by 6 m per year at the basincenter (under isothermal conditions). The headdrop moves fluid away from the step. The fluidsupply in the upper reservoir is limitless. Theperiod of interest is 10 years”.The corresponding FE mesh has 2967elements excluding the bedrock step (Figure 3).The rock is Hookean, poroelastic and isothermal.In the first example we use the same data of theES Module user’s guide, except for the BiotWillis coefficient we assume that b 0.3.

Figure 4. Poroelastic deformation of the basin for theBS problem with cold water (20 C). Streamlinesrepresent the fluid to porous rock coupling.Figure 5. Poroelastic deformation of the basin for theBS problem with hot water (250 C). Streamlinesrepresent the fluid to porous rock coupling.Figure 6. Vertical strain at the basin with a BS. Caseof cold water (20 C).Figure 7. Vertical strain at the basin with a BS. Caseof hot water (250 C).Figure 8. Horizontal strain at the basin with a BS.Case of cold water (20 C).Figure 9. Horizontal strain at the basin with a BS.Case of hot water (250 C).

6. Discussion of ResultsThe two examples presented herein weresolved using COMSOL–Multiphysics for a wellknown problem of linked fluid flow and soliddeformation near a bedrock step in a sedimentarybasin described in the Earth Science module. Theproblem concerns the impact of pumping for abasin filled with sediments draping animpervious fault block. In the first example, weconsidered the water in the aquifer to be cold, at20 C. In the second example, the water is hot, at250 C. The basin is composed of three layershaving a total depth of 500 m and is 5000 m longin both cases.Figures (4) and (5) show simulation resultsfor years 1, 2, 5, and 10, respectively. Thedifference here with the results of COMSOLMultiphysics is that we have performed thesimulation of a thermoporoelastic - coupleddeformation when the water in the aquifercorresponds to geothermal conditions (fluiddensity of 800.4 kg/m3, temperature of 250 C,and pressure of 50 bar). Figures (6) and (7)compare the vertical strains and figures (8) and(9) compare the horizontal strains, in both casesrespectively. Figures (8) and (9) also illustratethe evolution of lateral deformations thatcompensate for the changing surface elevationabove the bedrock step. Note that vertical scalesare different in both examples for clarity.7. Conclusions All crustal rocks forming geothermalreservoirs are poroelastic and the fluid presenceinside the pores affects their geomechanicalproperties. The elasticity of aquifers and geothermalreservoirs is evidenced by the compressionresulting from the decline of the fluid pressure,which can shorten the pore volume. This reductionof the pore volume can be the principal source offluid released from storage The immediate physical experience showsthat the supply or extraction of heat producesdeformations in the rocks. Any variation oftemperature induces a thermo-poroelasticbehavior that influences the elastic response ofporous rocks. We introduced herein a general tensorialthermoporoelastic model that takes into accountboth the fluid and the temperature effects inlinear porous rock deformations, and presentingtwo practical examples solved with COMSOL. The second example illustrates the influenceof temperature changes on the poroelastic strains.For cold water, the estimated value of εz is about-1.5ä10-4, while for hot water εz is -7.5ä10-4.Therefore, the poroelastic deformations are muchhigher in geothermal reservoirs than inisothermal aquifers. In the first case the bulkmodulus of water Kw 0.45 GPa, correspondingto T 250 C. For cold aquifers Kw 2.5 GPa. Water bulk modulus affect other poroelasticcoefficients, including the expansivity of rocks,which is relatively small, but its effects canproduce severe structural damages in porousrocks subjected to strong temperature gradients,as happens during the injection of cold fluids.8. References1. Biot, M.A., General Theory of ThreeDimensional Consolidation, Journal of AppliedPhysics, 12, , 155-164, (1941).2. Biot, M.A. and Willis, D.G., The ElasticCoefficients of the Theory of Consolidation, J. ofApplied Mechanics, 24, 594-601, (1957).3. Bundschuh, J. and Suárez M. C., Introductionto the Modeling of Groundwater and GeothermalSystems: Fundamentals of Mass, Energy andSolute Transport in Poroelastic Rocks. Taylor &Francis CRC Press, (2009, in press).4. Coussy, O., Mécanique des Milieux Poreux,Ed. Technip, Paris, (1991).5. Leake, S.A. and Hsieh, P.A., Simulation ofDeformation of Sediments from Decline ofGround-Water Levels in an Aquifer Underlainby a Bedrock Step. US Geological Survey OpenFile Report 97 – 47, (1997).6. Liu, G.R. and Quek, S., The Finite ElementMethod – A Practical course, Butterworth –Heinemann, Bristol, (2003).7. Terzaghi, K., Theoretical Soil Mechanics.John Wiley, New York, (1943).8. Wang, H.F., Theory of Linear Poroelasticity with Applications to Geomechanics andHydrogeology. 287 pp. Princenton UniversityPress, (2000).

processes are deduced using the Gibbs thermo-poroelastic potential or available enthalpy per unit volume and the energy dissipation function of the skeleton (Coussy, 1991). Analytic expressions are constructed in terms of the stresses, the porosity, the pore pressure, and the density of entropy per unit volume of porous rock.

Related Documents:

The COMSOL Multiphysics User’s Guide and COMSOL Multiphysics Reference Guide describe all interfaces and functionality included with the basic COMSOL Multiphysics license. These guides also have instructions about how to use COMSOL Multiphysics and how to access the documentation electronically through the COMSOL Multiphysics help desk.

The LiveLink interface is based on the COMSOL client/server architecture. A COMSOL thin client is running inside MA TLAB and has access to the COMSOL API through the MATLAB Java interface. Model information is stored in a model object available on the COMSOL server. The thin client communicates with the COMSOL

import and export functionality available in COMSOL Multiphysics and a detailed description of the import, export, and native COMSOL formats. For more detailed information on how to use COMSOL Multiphysics for import and export of data, see the COMSOL Multiphysics Reference Manual in the installed documentation set. Overview of Data Formats

COMSOL 5.3 Growing at a double-digit rate '14 COMSOL 5.0 '15 '16 5.1 '01 Started sales in Japan '17 COMSOL 5.2 COMSOL, Inc. General agency in Japan: Keisoku Engineering System Co., Ltd. (KESCO) 3 Multiphysics and Single-Physics Simulation Platform Mechanical, Fluid, Electrical, and Chemical Simulations Multiphysics -Coupled Phenomena

Introduction 1 Introduction Read this book if you are new to COMSOL Multiphysics .It provides a quick overview of the COMSOL environment with examples that show you how to use the COMSOL Desktop user interface. If you have not yet installed the software, install it now according to the

Use of COMSOL Multiphysics to Simulate RF Heating of Passive Conductive Implants in MRI Scanners Alan Leewood, PhD MED Institute, Inc. West Lafayette, IN . COMSOL Users’ Conference . Boston, MA . October 3, 2012 . Excerpt from the Proceedings of the 2012 COMSOL Conference in Milan

In this guide, the term COMSOL 4.4 refers not only to COMSOL Multiphysics 4.4 but also to any of the add-on products. Before You Begin 3. Software License Agreement. For the complete software license agreement with details on the license

Comsol Simple start 2D quasi steady RANS flow model Adding diluted species as progress variable Use Comsol Solve flow stationary Add time integration to stabilize the turbulent flame brush Comsol Conference 2013 31-10-2013 PAGE 10