1.Introduction

Biot’s theory of poroelasticity(Biot,1941)has gained new prominence in geotechnical Engineering to understand the coupled response of fluid f l owand deformation in deformable porous media(i.e.soils and rocks).The most common examples of this coupled hydro-mechanical(H-M)interaction are consolidation and subsidence induced by fluid extraction from underground formations(e.g.oil and gas from hydrocarbon reservoirs and water from aquifers).In these examples,the transient fluid f l ow affects deformation in the ground and vice versa(Wang,2000;Gutierrez and Lewis,2002;Neuzil,2003).Thus,consideration of this H-M interaction is essential for the safe design of structures built above or in saturated ground such as circular footing and deep tunnel.

Among various techniques to solve the coupled H-M equations,explicit technique is the most attractive one to be used(Minkoff et al.,2003;Dean et al.,2006).The technique not only is simple to be implemented,but also requires less time in building the separate mechanical and fluid f l ow codes for solving the coupled Biot’s equations that govern the response of the fluid and the porous medium.Unfortunately,explicit techniques come at a price:they are only conditionally stable.The nature of an explicit technique requires that small time step sizes must be used to maintain numerical stability,placing a strong limitation on the technique.Consequently,the efficiency of computer runtime for an explicit-type coupling approach cannot be fully exploited for largescale and long-term H-M simulations that may need enormously long computation time.For example,the H-M response of tunneling in low-permeability ground may need a consolidation time up to 250 years to reach the steady-state pore pressure distribution(Prassetyo and Gutierrez,2016a,b).

One of the most widely used explicit finite difference(FD)codes for H-M simulation in geotechnical engineering is the fast Lagrangian analysis of continua(FLAC)developed by Itasca(2011a).It is also conditionally stable.The time step size used for its f l ow calculation must be lower than a critical value in order that the pore pressure behavior remains stable and monotonic(Itasca,2011b).Unlike the fluid f l ow effect that takes place over a considerable amount of time,the geomechanical calculation in FLAC is not of particular concern because the response occurs almostinstantaneously particularly forstaticgeomechanical problems.

2.Need for an efficient explicit coupling technique

To improve the efficiency of an explicit coupling technique,one proposed method is to use an unconditionally stable fluid f l ow scheme that can be sequentially coupled with an existing geomechanical simulator such as FLAC.The alternating direction explicit(ADE)is one of such schemes.In an ADE scheme,two explicit FD equations are executed simultaneously in two physical directions:one in a forward sweep and the other in a reverse sweep(see Fig. 1).

To obtain the mathematical sense of how the ADE scheme works,the following diffusion equation is considered:

Having performed the operation presented in Eqs.(13)-(34),the novel high-order axisymmetric ADE scheme for the axisymmetric boundary points can be written for the forward sweep as

As presented in Barakat and Clark(1966),to solve Eq.(1),the ADE scheme executes two explicit FD equations in a forward sweep(Eq.(2))and in a reverse sweep(Eq.(3)):

Fig. 1.Illustration of(a)the forward and(b)the reverse sweeps of the standard plane strain ADE scheme.

where a,b and c are the constants defined as

4.2.2.Consistency

The final solution is the average of the solutions from the two sweeps defined as

It can be seen from Eqs.(2)and(3)and Fig. 1 that to solve a target point(xi,yj)at the current time level(n+1),both sweeping procedures use a combination of the solutions from the previous time level(n)and the current time level(n+1).The difference,however,is indicated by the location of the solutions from the current time level(n+1)with respect to the location of the target point.In the forward sweep(Eq.(2)),the calculation marches upward in an ascending order of x and y in which the solutions from the current time level(n+1)are located behind the target point(Fig. 1a).On the other hand,in the reverse sweep(Eq.(3)),the calculation marches downward in a descending order of x and y in which the solutions from the current time level(n+1)are located ahead of the target point(Fig. 1b).In this way,the explicitness of the scheme is preserved while improving the stability of the scheme.This is the merit of an ADE scheme.

From the state of geomechanical equilibrium,the solution to a coupled H-M problem is obtained sequentially by fi rst solving the fl ow problem and then the geomechanical problem.The new fourth-order axisymmetric ADE scheme is called to solve fl ow problem,resulting in the pore pressure solution pn+1.Note that the increment of pore pressure due to volumetric strain is evaluated in the geomechanical calculation as a zone value which is then distributed to the grid points.The total stress correction due to this pore pressure change is then calculated and updated into the geomechanical calculation.At the end of this geomechanical step,the displacement solution un+1is obtained.Several geomechanical steps may be taken to bring the system to equilibrium before the computation progresses to the next time step.This sequential procedure is repeated until the desired simulation time is reached.

Therefore,solving Eqs.(17)-(19)for∂p/∂r，∂2p/∂r2，and ∂2p/∂z2，respectively,gives

如今她把一切都告诉顾盼，最后她说，你知道了吧，你是警察，我却是个贼，你有大好前途，别和我搅在一起了。说完就把他推出去，喀一声锁了房门。任他在外面敲了良久，铁了心不开。

3.Governing equations for a coupled H-M problem

According to Biot’s theory of poroelasticity(Biot,1941),the coupled equations governing the responses of the fluid(pore pressure)and the solid(effective stress)are related to the changes in fluid content(ζ)and volumetric strain(εvol)of the porous medium.This relation can be written as

where M is the Biot modulus;α is the Biot coefficient;σijis the total stress tensor;δijis the Kronecker delta tensor;and K and G are the drained bulk and shear moduli of the porous medium,respectively.The change in fluid content can be related to the change in pore pressure bysubstituting Darcy’s law in Eq.(8)into the mass balance equation in Eq.(9).After neglecting the gravitational term,the resulting relation is given in Eq.(10).

where k is the mobility coefficient（k=kH/γw），kHis the hydraulic conductivity,γwis the unit weight of fluid,ρwis the fluid density,and g is the gravitational acceleration vector.

In this paper,the contribution from the mechanical deformation to the incrementof pore pressureis evaluated in the geomechanical calculation.Thus,during the f l ow calculation, δεvol=0.By substituting Eq.(10)into Eq.(6),the transient response of the fluid is then governed by the Laplacian of the pore pressure as

where cvis the fluid diffusivity defined as cv=kM.

Eq.(11)is the partial differential equation(PDE)of the f l ow problem for a plane strain condition(i.e.for a coordinate system adopted in Plane A in Fig. 2).For an axisymmetric condition,the PDE in Eq.(11)may be rewritten in the form of Eq.(12),which corresponds to the adopted coordinate system shown in Plane B in Fig. 2:

The axisymmetric ADE scheme in Eqs.(33)and(34)is valid for all points inside the FD grid except for the points along the axis of symmetry or at r=0(see Fig. 5).The PDE in Eq.(12)has a singularity at r=0 because the equation encounters a division of 0/0 for the （1/r）∂p/∂r term.To remove this singularity,L’Hospital’s rule must be used.The numerator and the denominator of the term（1/r）∂p/∂r must be successively differentiated with respect to r until a finite limit is achieved.Eq.(35)shows the operation of L’Hospital’s rule for the （1/r）∂p/∂r term,in which the term simply becomes ∂2p/∂r2:

Fig. 2.Illustration of the coordinate systems used for plane strain(Plane A)and axisymmetric(Plane B)analyses.

where f= ∂p/∂t，the variable r is in the radial direction,while z is in the axial direction.The pore pressure p in Eq.(12)is assumed to be independent in the angularθdirection.

4.Novel high-order axisymmetric ADE scheme

4.1.Development of high-order axisymmetric ADE scheme

Todeveloptheproposedhigh-orderaxisymmetricADE scheme for a non-uniform grid,the non-uniform FD approximations of the axisymmetric diffusion equation in Eq.(12)must be developed first.To do this,we consider a general non-uniform FD grid that consists of four quadrilateral elements in the axisymmetric coordinate system in Fig. 3.From Fig. 3,writing the Taylor series expansion for pi+1，jand pi-1，jto the fourth-order gives

Multiplying pi+1，jin Eq.(13)by Δr-and pi-1，jin Eq.(14)by Δr+,and disregarding the remainder terms gives

Fig. 3.Non-uniform FD grid for the axisymmetric coordinate system.

To obtain ∂p/∂r,Eq.(14)is subtracted from Eq.(13);while to obtain ∂2p/∂r2，Eq.(16)is added to Eq.(15).To obtain ∂2p/∂z2，the steps to get Eqs.(13)-(16)are repeated for z and the resulting equations are added.These operations give

The main objective of this paper is to remedy the drawbacks of the standard ADE scheme so that the new ADE scheme can be coupled with a geomechanical simulator for an efficient simulation of axisymmetric coupled H-M problems in non-uniform grids.This will be achieved by formulating a novel high-order ADE scheme capable of solving the axisymmetric fluid f l ow problem in non-uniform grids.Derivation of the new scheme will be done by performing a fourth-order FD approximation of the spatial of the derivatives axisymmetric fluid-diffusion equation in a nonuniform grid configuration.The resulting fluid f l ow solution is then sequentially coupled with the existing geomechanical simulator in FLAC.Because the time step restriction is now removed,the H-M simulation can now be performed efficiently with less computation time while still maintaining numerical stability and high numerical accuracy.This paper is a major extension of the newly-derived high-order ADE scheme for solving plane strain fl ow problem in non-uniform grids(Prassetyo and Gutierrez,2017a,b).

In order to have the desired fourth-order approximations for∂p/∂r，∂2p/∂r2，and ∂2p/∂z2，the higher-order partial derivatives∂3p/∂r3，∂4p/∂r4，∂3p/∂z3，and ∂4p/∂z4in Eqs.(20)-(22)need to be approximated to the second-order accurate ones.To accomplish this,we perform repeated differentiations of Eq.(12)with respect to r and z until the higher-order terms appear.The differentiations give

The next step is to use the centered difference approximations in Eq.(23)for the lower-order partial derivatives in Eqs.(24)-(27)and substitute those approximations into Eqs.(20)-(22)and the governing equation in Eq.(12).By performing this operation,the desired fourth-order FD approximations of the axisymmetric diffusion equation can be obtained as

where

To assemble the fourth-order approximation of the axisymmetric ADE scheme,the well-known Crank-Nicolson implicit time integration scheme(Crank and Nicolson,1947)is applied to Eq.(28).The principle of the Crank-Nicolson technique is that difference approximations are developed at the midpoint of the time increment n(Chapra and Canale,2009).To do this,the temporal derivative ∂p/∂t is approximated at tn+1/2as

关于组织起来的意义，毛泽东指出：“在农民群众方面，几千年来都是个体经济，一家一户就是一个生产单位，这种分散的个体生产，就是封建统治的经济基础，而使农民自己陷于永远的穷苦。克服这种状况的唯一办法，就是逐渐地集体化。”［4］931毛泽东认为，通过合作社，“我们就可以把群众的力量组织成为一支劳动大军。这是人民群众得到解放的必由之路，由穷苦变富裕的必由之路，也是抗战胜利的必由之路”［4］931。

while the spatial derivatives ∂2p/∂r2and ∂2p/∂z2can be approximated at the midpoint by averaging the spatial derivatives at the beginning tnand at the end tn+1of the time increment as

Finally,having performed the operation,the novel high-order axisymmetric ADE scheme for a non-uniform grid can be written for the forward sweep as

and for the reverse sweep as whereα0toα8and β1to β4are the grid size-dependent constants.The expressions for these constants are given in Eqs.(A1)-(A5)in the Appendix to prevent encumbering the mainpresentation in the text.These constants are made up of more terms than those in Eq.(4)to account for non-uniformity in grid size.Fig. 4 illustrates the arrangement of the calculation sweeps for the newly-derived highorder axisymmetric ADE scheme.Compared to the standard scheme in Fig. 1,the new ADE scheme involves more contributing points in the calculation of the target point because of the fourthorder approximation,e.g.additional diagonal points from the cross-derivatives.Nevertheless,the new scheme still maintains the nine-point stencil,preserving its compactness.

“北接松径，南通峦雉，东以达虎角庵。游者之屦常满，然而素桷茅榱，了不异人意。”[3]426此亭构造简朴，一点也不吸引游者的眼球。然而登亭眺望，胜景扑面而来，使人油然而生山水鱼鸟之情。

After substituting ∂2p/∂r2for （1/r）∂p/∂r in Eq.(12),the PDE becomes that shown in Eq.(36).This is the governing PDE that needs tobe discretized toobtain the high-order ADE scheme for the boundary points located along the axis of symmetry(called the axisymmetric boundary points from now on).

For the axisymmetric boundary points,the steps to develop the corresponding high-order axisymmetric ADE scheme are the same as those for the interior points.The steps will not be repeated here but the final arrangements will be given.

where p is the pore pressure.

Note that the constants α1，α5，α7，and β1are now equal to zero,meaning that fewer neighboring points are needed in the calculation.The expressions for grid size-dependent constants for the axisymmetric boundary points are given in Eq.(A6)in the Appendix.

Note that in the above derivations,the mobility coefficient k=kH/γwhas been assumed to be constant,causing the fluid diffusivity cvin the grid size-dependent constants to be constant as well.The permeability of a porous medium can certainly be made anisotropic by having different mobility coefficients in the radial and axial directions(krand kz,respectively).To accommodate this anisotropy,Eqs.(12)and(36)can be rewritten in more general forms as

Fig. 4.Illustration of(a)the forward and(b)the reverse sweeps of the new high-order axisymmetric ADE scheme.

4.2.Proof of convergence

Convergence is the attribute that the newly-derived ADE scheme must have to instill conf i dence in its results.According to the Lax-Richtmyer equivalence theorem,stability and consistency are the conditions to be satisf i ed by an FD scheme in order for it to be called a convergent scheme(Richtmyer and Morton,1967).Stability means that errors at any stage of the computation decay as the computation progresses,while consistency means that the FD scheme reduces to the governing PDE as the increments in the independent variables vanish.

Fig. 5.Illustration of an axisymmetric consolidation showing location of the axis of symmetry.

4.2.1.Stability

To prove the stability of the new ADE scheme,the von Neumann stability analysis is used.The new ADE scheme can be said to be stable if the amplif i cation factor of the solution meets the criterion|G（θ）| 1 when the solution is advanced one time step,expressing the decay of the frequency amplitude of the solution(Barakat and Clark,1966;Hoffman,2001;Strikwerda,2004;Mazumder,2016).To use the von Neumann stability analysis,the pore pressure solutions to the forward and reverse sweeps are rewritten in separable Fourier expansions as

For the interior points,the substitution yields the expression for calculating G(θ)for the forward sweep as

and for the reverse sweep as

Likewise,for the axisymmetric boundary points,the substitution yields the expressions for calculating G(θ)for the forward sweep as

and for the reverse sweep as

To check if the criterion for stability is met for each sweep,i.e.|G（θ）| 1，graphical examination of G(θ)will be helpful.To do this,the material properties from the axisymmetric consolidation problem in Section 5 are used to calculate the grid size-dependent constants from Eqs.(A1)-(A6).The resulting values of the constants are substituted into Eqs.(42)-(45)and the resulting values of G(θ)with varying θare then plotted.Fig. 6 shows the plots of G(θ)for the forward and reverse sweeps for the interior points(Fig. 6a)and for the axisymmetric boundary points(Fig. 6b).It can be seen that as θvaries,the set of points marked out by G(θ)lies within the unit circle with the radius of unity,ensuring that the FD solutions are bounded and that the new ADE scheme is unconditionally stable.In either interior or axisymmetric boundary points,the shape of the G(θ)curve for each sweep is not the same because the equations for G(θ)for each sweep are also different,resulting in different values of grid-size dependent constants.

据悉，中国质量奖是中国质量领域最高的政府性荣誉，于2012年设立，设中国质量奖和中国质量奖提名奖，每两年评选一次，旨在表彰国内质量领先、技术创新、品牌优秀、效益突出的组织和对促进质量发展作出突出贡献的个人和企业。（江源荐，黄筱鹂编辑）

正式审议阶段的主要任务是审查和讨论法案，这是大多数国家立法机关在审议法案中必经的步骤。由于我国不存在立法辩论制度，因而专门委员会的立法审查在正式审议阶段发挥着举足轻重的作用，他们负责对所提出的法案进行检查核对，目的在于发挥议员的专业性优势，保证立法的科学性与合法性。1979年后，逐渐恢复了1954年宪法确立的委员会制度，并将立法审查制度化——实行统一审议。统一审议的主要功能在于维护法制统一和扼制不当利益。[14]1982年通过的《全国人民大表大会组织法》首次在法律层面对专门委员会的分工予以确定，并建立了以法律委员会为主，专门委员会为辅的统一审议制度。

2.放养准备。认真进行池塘整理，去除过多淤泥、平整池底，在池塘坡底四周挖沟，利于早期虾苗培育和成虾捕获。在养殖水体设置埋入土下0.2m、上露0.5m内壁光滑的石棉瓦做防逃措施，尤其在暴雨洪涝时期防止冲垮池埂、漫坝引发逃虾。建好进排水系统。保持池水深0.8～1.5m，中间水深，四周有浅滩，水草丰富。

The consistency condition requires that the truncation errors of the newly-derived ADE scheme vanish asΔr,Δz,and Δt approach zero.The behavior of the truncation errors for the interior grids can be investigated by looking at Eq.(28)because the discretization of the high-order axisymmetric ADE scheme comes from this equation.Eq.(28)may be rewritten as

式中，ts和tf分别为硅基底和氧化硅薄膜的厚度，Es和νs分别为硅基底的杨氏模量和泊松比.当应力值呈负值时薄膜为压应力；当应力值呈正值时薄膜为张应力.对掺Ge氧化硅薄膜波导利用RIE刻蚀形成波导结构，用扫描电子显微镜(SEM, Hitachi S-4800型)获得了波导截面图与表面俯视图如图1所示，可见薄膜波导刻蚀角度较差，这是由RIE刻蚀的固有缺陷造成，但薄膜无裂纹无气泡质量较好.

where τi，jis the truncation error resulting from expanding the governing equation in Eq.(12)to the fourth-order one.The truncation error is defined as

Assuming that uniform grid sizes are used,Δr_= Δr+= Δr and Δz_= Δz+= Δz,the constants M,N,A,B,C,and D defined in Eq.(29)becomeΔr3,0,0,Δr2,0 and Δz2,respectively.Clearly,when these grid sizes approach zero, τi，jwill be zero and Eq.(46)reduces to the governing PDE defined in Eq.(12),guaranteeing the consistency of the new scheme.A similar approach can be taken for the axisymmetric boundary points in which when the truncation error goes to zero,the governing PDE in Eq.(36)will be recovered.

4.3.efficient sequential coupling technique

The coupled problems in this paper will be modeled and solved in FLAC.At every time step,the newly-derived axisymmetric ADE scheme will first be executed to solve the fluid f l ow problem,and then the system will be sequentially brought to geomechanical equilibrium using FLAC’s built-in geomechanical simulator.From now on,the use of the axisymmetric ADE scheme in FLAC is called the sequentially-explicit coupling technique based on the fourthorder axisymmetric ADE scheme or SEA-4-AXI.This coupling technique is written in FLAC’s programming language called FISH(Itasca,2011c).Fig. 7 shows the f l owchart of SEA-4-AXI to solve the coupled problems in this paper.

Despite its unconditional stability,the standard ADE scheme has three major drawbacks.First,it is only a second-order accurate scheme in time and space.This low-order accuracy will result in low numerical accuracy when a large time step size is used.Second,it is only valid for uniform grid sizes with the sameΔx andΔy for the entire grid model.This limitation prevents the standard ADE scheme from solving large-scale f l ow problems with increasing grid size towards the far-field boundaries.Third,it is only valid for plane strain conditions.This means that it cannot be used for solving f l ow problems under axisymmetric conditions.

It is worth mentioning that the geomechanical calculation in FLAC invokes the equations of motion to derive new velocities and displacements from stresses and forces.Strain rates are then derived from velocities,and new effective stresses from strain rates according to the adopted constitutive law.If the water bulk modulus is given,new increments of pore pressure due to pore volume change will also be evaluated from strain rates.This increment is called the mechanical generation of pore pressure.This cycle constitutes one geomechanical step.At the beginning of an H-M simulation,one fluid time step may putthesystem considerablyoutofequilibrium,inducing large unbalanced forces in the model.At this stage,many geomechanical steps will be taken for each fluid step.As the consolidation process continues,the changes in fluid pressure will become small,so that the system will remain in equilibrium(Itasca,2011b).

Fig. 6.Variation of G(θ)for the forward and reverse sweeps of the new axisymmetric ADE scheme for(a)the interior and(b)the axisymmetric boundary points.

In this paper,the accuracy and efficiency of SEA-4-AXI will be compared to FLAC’s built-in coupled technique called the basic scheme.The basic scheme uses the explicit f l ow formulation which requires that the fluid time step size should be smaller than a critical valueΔtcri.The value is calculated internally based on a discretization of the medium into zones composed of two overlaid triangular elements,in which the nodal fluid f l ux q is solved by using the Gauss divergence theorem(Itasca,2011b).The merit of SEA-4-AXI is that it removes this time step restriction,leading the explicit coupling technique to yield an efficient H-M simulation without suffering from numerical instability yet still maintaining high numerical accuracy.

Fig. 7.The sequential coupling technique in FLAC using the fourth-order axisymmetric ADE scheme(SEA-4-AXI).

5.Axisymmetric consolidation of a circular footing

Fig. 8.Geometry and non-uniform FD grid of the axisymmetric consolidation beneath a circular footing.

Table 1 Material properties for McNamee’s problem.

Hydraulic conductivity,kH(m/d)Porosity,φ Shear modulus,G(Pa)Drained bulk modulus,K(Pa)Fluid bulk modulus,Kw(Pa)Surface load,po(kPa)Diffusivity,cv(m2/d)Smallest grid size,Δx?Δy(m?m)Largest grid size,Δx?Δy(m?m)Number of zones 1.2?10-4 0.3 2.3?105 1.5?105 2?109 100 0.005 0.5?0.5 4?3.5 480

SEA-4-AXI is first validated to simulate the H-M response of an axisymmetric consolidation of a circular footing known as McNamee’s consolidation(McNamee and Gibson,1960),which is a commonproblem in geotechnicalengineering.A fullysaturated soil layer of height H=30 m and width W=40 m is loaded by a constant uniform pressure of po=100 kPa over a circular area of radius a=2 m.The soil is free to drain at the surface but the other boundaries are impermeable.The FD grid is gradually increasing in size toward the right and bottom boundaries of the model,testing the validity of the non-uniform FD approximation in the new axisymmetric ADE scheme(Fig. 8).The material properties for the model are given in Table 1.The model is run to a consolidation time of t=250 years,which is about 112 times the normalized consolidation time t*for this problem(t*=cvt/a2).The simulation in SEA-4-AXI uses Δt=20 d,while that in FLAC’s basic f l ow scheme uses Δtcri=5 d.

Note that the paper adopts the sign convention used in rock mechanics that stress and displacement are positive for compression.Unless stated otherwise,all pore pressures are total(initial hydrostatic plus excess pore pressure)and all stresses are effective.Further,the H-M simulations performed in this paper use constant Biot’s coefficientα =1,and the coefficient does not change with the mechanical deformation.This is done to justify the two simplif i ed assumptionsin thispaper:theground iselasticand the compressibility of grains is neglected compared to that of the drained bulk modulus K.The effects of usingα 1 and of varyingα with deformation on the H-M response of consolidation problems can be seen in the study of tunneling in poroplastic medium done by Bobet(2010).

Fig. 9.Contours of(a)pore pressure,(b)effective vertical stress,and(c)radial displacement after 25 years of consolidation(comparison between SEA-4-AXI and FLAC’s basic scheme).

Fig. 10.Histories of(a)pore pressure and(b)degree of consolidation at various monitoring points.

Fig. 9 shows the contour plots for the H-M responses after the simulation has run for 25 years(t*=11).As one can see,the comparison of the induced H-M responses,i.e.pore pressure,displacement and stress,between SEA-4-AXI and FLAC’s basic scheme is very satisfactory.To observe the long-term H-M behaviors,the pore pressure and settlement histories are monitored at various depths(Fig. 10a)and at various lateral distances from the axis of symmetry(Fig. 10b),respectively.The pore pressure is plotted as normalized pore pressure(p/po),while the settlement is plotted as the degree of consolidation,U(%).An excellent match is observed in the solutions between SEA-4-AXI and FLAC’s basic scheme and those between SEA-4-AXI and the analytical calculation given by McNamee and Gibson(1960).As an accuracy comparison,Table 2 compares the degree of consolidation between the numerical solution and the analytical solution at various normalized consolidation times t*.As can be seen,the average percentage error of the SEA-4-AXI solution is 0.17%while its maximum percentage error is less than 0.5%.

Other excellent agreements are observed in the distribution of pore pressure with depth and of settlement with distance from axisymmetric line.At various consolidation times from t*=0.5 to t*=11,the pore pressures with depth from SEA-4-AXI match very well with those from FLAC’s basic scheme(Fig. 11).Likewise,the agreement on the settlement profiles at t*=0.1 and 112 between SEA-4-AXI and FLAC’s basic scheme is excellent(Fig. 12).In the settlement profile,the otherhalf of the circular footingmodel(from0 m to-40 m)is projected symmetrically as the side being modeled.

Table 2 Comparison of the degree of consolidation U(in%)between SEA-4-AXI and the analytical solution at various normalized consolidation times t*.

t* Degree of consolidation,U(%) Error(%) Average error(%)SEA-4-AXI Analytical 0.1 43.01 43.03 0.03 0.17 0.5 72.33 72.23 0.14 2.8 89.32 89.09 0.26 5.6 92.97 92.54 0.47 11.2 95.14 94.97 0.18 22.5 96.75 96.81 0.06 45 98.13 98.35 0.22 100 100 100 0

Fig. 11.Profiles of pore pressure with depth at various consolidation times.

Fig. 12.Profiles of settlement at early(t*=0.1)and after long-term(t*=112)consolidations.

Fig. 13.The advancing tunnel problem showing(a)geometry and(b)excavation steps.

Because SEA-4-AXI uses a larger time step size than the basic scheme does,it is very efficient in performing the coupled simulation.The computer runtime for this coupled problem is significantly reduced by 50%,from 120 s in FLAC’s basic scheme to only 60 s in SEA-4-AXI.

6.H-M analysis of an advancing tunnel in deep saturated ground

The potential application of SEA-4-AXI to geotechnical engineering is now demonstrated by simulating the H-M response of an advancing tunnel in deep saturated ground.The advancing tunnel is represented by a two-dimensional(2D)axisymmetric model of an unlined circular excavation of radius R=2.5 m and at a depth of H=225 m.The tunnel is excavated through an elastic and saturated medium subjected to total isotropic in situ stress of σo = 4.5 MPa and initialhydrostatic pore pressure of po=2.25 MPa.Tunnel at this depth has been reported by Graziani and Boldini(2012)when investigating the influence of H-M coupling on deep tunnel response in clays.Note that an axisymmetric model with a circular boundary and isotropic stress is an appropriate approximation for very deep tunnels where the halfspace conditions corresponding to the ground surface can be disregarded.

The FD grid is also gradually increasing in size toward the right and top boundaries of the model(Fig. 13a).The model boundaries are impermeable except the tunnel wall and face which are made permeable.The advance of the tunnel starts from x/D=10 to the final tunnel face at x/D=0.The advance is simulated by progressively deactivating 40 FD regions of thicknessΔx=1.25 m(D/4)over an excavation period of 10 d,implying the excavation rate of va=5 m/d(Fig. 13b).

During deactivation of each regionΔx,a practically instantaneous face advance is being simulated for an excavation period of Δt=0(undrained loading),followed by a drained consolidation for a period ofΔt= Δx/va(drained loading).This step-by-step excavation,represented by the alternating stages of undrained excavation and drained consolidation,is performed until the final tunnel face at x/D=0 is reached.The response of the ground surrounding the tunnel shortly after the final tunnel face is reached is called the short-term response.This state also corresponds to the start of the standstill period(t*=0).This long-term consolidation is performed until the total period of 930 d(t*=100)is reached.The normalized consolidation time t*is calculated as t*=R2/cv(Giraud and Rousset,1996;Li,1999).The groundproperties for this problem are given in Table 3.

Table 3 Ground properties for advancing tunnel problem.

Hydraulic conductivity,kH(m/s)Porosity,φ Shear modulus,G(Pa)Drained bulk modulus,K(Pa)Fluid bulk modulus,Kw(Pa)Diffusivity,cv(m2/s)Smallest grid size,Δr? Δx(m?m)Largest grid size,Δr?Δx(m?m)Number of zones 5?10-10 0.4 1.3?108 1.3?108 2?109 2.6?10-4 0.5?1.25 5?2.5 1560

Fig. 14.Contours of(a)pore pressure,(b)radial displacement,(c)effective radial and(d)tangential stresses after the short-term consolidation(comparison between SEA-4-AXI and FLAC’s basic scheme).

The simulation in SEA-4-AXI usesΔt=20 minwhile that in FLAC usesΔtcri=5 min.As expected,SEA-4-AXI reduces the computer runtime by 42%,from 520 s in FLAC’s basic scheme to only 300 s in SEA-4-AXI.Discussion in this sectionwill be mainlyon the accuracy comparison between the results from SEA-4-AXI and those from FLAC’s basic scheme.

因此，家长要从全局出发，在生活中引导孩子自发地训练自控力，而不是纠结于满不满足孩子的要求。如果要求合理，就没有必要推迟，总是故意不满足孩子的要求，盲目地延迟满足，反而会伤害亲子关系。

要解决的技术问题是RFID阅读器清点标签时，多轮清点协同工作的问题，通过协同工作避免重复清点标签，从而减少标签识别时的查询次数，提高标签清点效率。

Fig. 14a-d showsthecontoursofporepressure,radial displacement and stress fields in the short term.The agreement between the contours from SEA-4-AXI and those from FLAC’s basic scheme is very satisfactory.Near the tunnel face,the instantaneous pore pressure develops within a short distance of-0.5D toward the advance core(Fig. 14a).For a detailed examination of the induced H-M response in the short-and longterm consolidations,the longitudinal profiles of the pore pressure,radial displacement and effective stresses at the tunnel wall are plotted along the axis of the tunnel(Fig. 15).As can be seen previously,there appears to be a short-term pore pressure buildup up to 2.65 MPa in the region just ahead of the tunnel face(Fig. 15a).In this region,the pore pressure ahead of the core does not have time to dissipate at a speed comparable to the excavation progress due to the nature of ground diffusivity.As drainage proceeds with time,the steep pore pressure gradient decreases,transitioning the pore pressure in the ground towards its steady state as represented by the long-term profile.This non-monotonic behavior and increase of pore pressure above the initial hydrostatic value are analogous to the so-called Mandel-Cryer effect(Mandel,1953;Cryer,1963;Abousleiman et al.,1996;Gutierrez and Lewis,2002).In tunneling,the Mandel-Cryer effect was also observed at the tunnel crown and was caused by the stress concentration in the stiffer zone that is still undrained(Prassetyo and Gutierrez,2016a,b).Due to the dissipation of pore pressure,the radial displacement also increases with time(Fig. 15b)as well as the effective stresses(Fig. 15c).

Fig. 15.Longitudinal profiles of(a)pore pressure,(b)radial displacement,and(c)effective stresses at the tunnel wall after the short-and long-term consolidation.

It can also be seen from Fig. 15 that the agreement between the H-M responses from SEA-4-AXI and FLAC’s basic scheme is excellent with the average percentage difference well below 1.8%in both short-and long-term consolidations.As shown in Tables 4-7,the average percentage differences for pore pressure,radial displacement,effective radial stress,and effective tangential stress in the short term are 0.67%,0.58%,1.31%,and 0.41%,respectively,whilethose in the long term are 1.79%,0.78%,0.92%,and 0.61%,respectively.

Table 4 Comparison of the pore pressure p along the axis of the tunnel between SEA-4-AXI and FLAC’s basic scheme.

x/D p from SEA-4-AXI(MPa)p from FLAC’s basic scheme(MPa)Difference(%)Longterm Shortterm Shortterm Longterm Shortterm Longterm-10 2.25 2.18 2.25 2.16 0 0.83-9 2.25 2.17 2.25 2.15 0 0.81-8 2.25 2.16 2.25 2.14 0 0.81-6 2.25 2.11 2.26 2.09 0.44 0.93-5 2.26 2.06 2.26 2.04 0.2 1.02-4 2.26 1.99 2.26 1.97 0.09 1.14-3 2.27 1.88 2.27 1.86 0 1.29-1 2.61 1.36 2.65 1.33 1.51 2.12 0 1.5 0.54 1.5 0.51 0 6.05 1 1.04 0.35 1.07 0.35 3.15 1.87 1.5 0.99 0.33 1.01 0.33 1.37 1.63 2 0.95 0.32 0.96 0.32 0.95 1.52 2.5 0.91 0.31 0.91 0.31 0.26 1.48 3.75 0.82 0.3 0.83 0.3 0.67 1.48 5 0.76 0.29 0.77 0.29 0.86 1.64 6.25 0.71 0.29 0.72 0.29 1.09 1.8 7.5 0.68 0.29 0.69 0.28 0.68 2.12 8.75 0.66 0.29 0.66 0.28 0.64 2.48 10 0.65 0.29 0.65 0.28 0.76 2.9 Average 0.67 1.79

To capture the transient H-M response during the excavation period,a monitoring point is placed at x/D=5(see the top inset inFig. 16).The excavation starts at x/D=10(at t=0)and f i nishes at x/D=0(at t=10 d).The agreement in the induced H-M responses with time between SEA-4-AXI and FLAC’s basic scheme is again very satisfactory(Fig. 16).SEA-4-AXI can capture the excess pore pressure that develops before the tunnel face arrives at the monitoring point at t=3 d(at x/D=7)before it falls sharply as the tunnel face passes the point.Similar behavior is observed for the convergence history in Fig. 16b,while the corresponding effective stresses are shown in Fig. 16c.Note that the jagged shape in the plot of the H-M response is due to the application of the alternating undrained and drained loadings during the tunnel excavation.This wavy response is expected because the face is advanced at practically instantaneous excavation periodΔt=0(the vertical segment)before it is allowed to consolidate for the period ofΔt=0.25 d(the horizontal segment).The jagged shape has also been reported in the previous coupled analysis of tunneling in shallow ground(Callari and Casini,2005).

Table 5 Comparison of the radial displacement uralong the axis of the tunnel between SEA-4-AXI and FLAC’s basic scheme.

x/D urfrom SEA-4-AXI(cm)urfrom FLAC’s basic scheme(cm)Difference(%)Longterm Shortterm Shortterm Longterm Shortterm Longterm-10 0.1 0.2 0.1 0.2 0 0-9 0.1 0.2 0.1 0.2 0.00 0-8 0.1 0.24 0.1 0.24 0 0.57-6 0.1 0.29 0.1 0.28 1.14 2.44-5 0.11 0.3 0.11 0.3 0.74 0-4 0.1 0.38 0.1 0.37 0 2.70-3 0.1 0.46 0.1 0.45 0 2.22-1 0.28 0.73 0.27 0.71 3.16 2.63 0 1.15 1.37 1.15 1.36 0.48 0.59 1 3.83 4.02 3.9 4.07 1.93 1.25 1.5 4.02 4.2 4.06 4.23 1.06 0.71 2 4.1 4.28 4.13 4.3 0.68 0.49 2.5 4.14 4.33 4.16 4.34 0.42 0.34 3.75 4.19 4.4 4.2 4.41 0.23 0.13 5 4.2 4.45 4.21 4.46 0.27 0.2 6.25 4.2 4.48 4.21 4.49 0.22 0.15 7.5 4.2 4.50 4.21 4.51 0.23 0.17 8.75 4.19 4.51 4.21 4.52 0.25 0.18 10 4.22 4.55 4.23 4.55 0.11 0.04 Average 0.58 0.78

Table 6 Comparison of the effective radial stress σ′ralong the axis of the tunnel between SEA-4-AXI and FLAC’s basic scheme.

x/D σ′rfrom SEA-4-AXI(MPa) σ′rfrom FLAC’s basic scheme(MPa) Difference(%)Short-term Long-term Short-term Long-term Short-term Long-term-10 2.25 2.37 2.25 2.38 0.08 0.57-9 2.25 2.37 2.25 2.38 0.08 0.64-8 2.25 2.38 2.25 2.4 0.06 0.7-6 2.25 2.42 2.25 2.44 0.36 0.87-5 2.24 2.45 2.26 2.48 0.63 0.9-4 2.26 2.5 2.26 2.52 0.34 0.77-3 2.28 2.57 2.28 2.58 0.03 0.25-1 2.07 3.01 2 2.98 3.39 1.22 0 2.2 2.97 2.27 3.04 3.23 2.32 1 1.3 1.9 1.28 1.9 1.94 0.08 1.5 1.4 1.98 1.36 1.94 3.05 1.91 2 1.46 2.01 1.41 1.97 3.17 1.97 2.5 1.49 2.03 1.46 1.99 2.58 1.89 3.75 1.55 2.03 1.54 2.02 0.65 0.3 5 1.61 2.06 1.59 2.04 1.38 0.99 6.25 1.65 2.07 1.62 2.05 1.68 1.06 7.5 1.67 2.08 1.66 2.07 0.8 0.36 8.75 1.71 2.1 1.69 2.09 1.05 0.57 10 1.75 2.14 1.74 2.14 0.34 0.07 Average 1.31 0.92

Table 7 Comparison of the effective tangential stress σ′θalong the axis of the tunnel between SEA-4-AXI and FLAC’s basic scheme.

x/D σ′θ from SEA-4-AXI(MPa) σ′θfrom FLAC’s basic scheme(MPa) Difference(%)Short-term Long-term Short-term Long-term Short-term Long-term-10 2.25 2.37 2.25 2.38 0.07 0.61-9 2.25 2.37 2.25 2.39 0.09 0.67-8 2.25 2.38 2.25 2.4 0.07 0.72-6 2.25 2.42 2.25 2.44 0.3 0.89-5 2.24 2.46 2.26 2.48 0.55 0.95-4 2.26 2.51 2.26 2.53 0.32 0.94-3 2.27 2.58 2.28 2.6 0.17 0.73-1 2.19 3.07 2.2 3.08 0.52 0.02 0 3.89 4.42 3.94 4.39 1.24 0.53 1 4.7 4.98 4.68 4.95 0.52 0.7 1.5 4.83 5.09 4.79 5.04 0.75 1.02 2 4.89 5.13 4.85 5.09 0.73 0.87 2.5 4.92 5.16 4.89 5.12 0.66 0.82 3.75 4.96 5.18 4.95 5.16 0.22 0.28 5 4.99 5.22 4.97 5.19 0.55 0.61 6.25 5.01 5.24 4.98 5.22 0.53 0.55 7.5 5.01 5.26 5 5.24 0.2 0.23 8.75 5.03 5.29 5.01 5.27 0.36 0.39 10 5.11 5.37 5.11 5.36 0 0.04 Average 0.41 0.61

During the standstill period,the convergence and the extrusion at various locations in the ground ahead of the face are also monitored(Fig. 17a and b,respectively).For graphical presentation,extrusion is regarded as positive deformation in this paper.As expected,the convergence and the extrusion located at the face(0D)are the largest.The convergenceδr0increases by 20%during the standstill period from 1.2 cm at t=10 d to 1.4 cm at t=930 d,while the face extrusion uex0increases up to 39%,from 3.3 cm to 4.5 cm.Fig. 17c shows the extrusion across the face height at both standstill periods.The face extrusions from SEA-4-AXI are very satisfactory in matching those from FLAC’s basic scheme with the average percentage differences of 0.85%and 0.63%at t=10 d and 930 d,respectively(Table 8).At other locations ahead of the face(at-1D and-2D),the convergence,δr1and δr2,and the core extrusion,uex1 and uex2,follow the same trend as those at the face but to a much less degree.

Fig. 16.Transient H-M response during the excavation period for(a)pore pressure,(b)radial displacement,and(c)effective stresses.

7.Conclusions

Fig. 17.Histories of(a)convergence and(b)extrusion during the excavation and standstill periods with the corresponding(c)face profiles.

Table 8 Comparison of the face extrusion uex0across the face height between SEA-4-AXI and FLAC’s basic scheme.

Face height(m)uex0from SEA-4-AXI(cm)uex0from FLAC’s basic scheme(cm)Difference(%)10 d 930 d 10 d 930 d 10 d 930 d 2.5 1.29 2.63 1.29 2.63 0.56 0.15 2 2.09 3.4 2.13 3.42 1.96 0.71 1.5 2.62 3.91 2.64 3.9 0.88 0.29 1 2.97 4.24 2.97 4.2 0.09 1.03 0.5 3.18 4.44 3.16 4.4 0.72 0.97 0 3.26 4.53 3.24 4.5 0.89 0.62-0.5 3.18 4.44 3.16 4.4 0.72 0.97-1 2.97 4.24 2.97 4.2 0.09 1.03-1.5 2.62 3.91 2.64 3.9 0.88 0.29-2 2.09 3.4 2.13 3.42 1.96 0.71-2.5 1.29 2.63 1.29 2.63 0.56 0.15 Average 0.85 0.63

The main goal of this paper was to eliminate the drawbacks in the standard ADE so that the new ADE scheme can be coupled with a geomechanical simulator for an ef fi cient simulation of axisymmetric coupled problems in non-uniform grid.This goal was achieved by presenting a novel high-order ADE scheme capable of solving the axisymmetric fluid fl ow problem in a nonuniform grid configuration.The new scheme was derived by performing a non-uniform fourth-order FD approximation of the spatial derivatives of the axisymmetric fluid-diffusion equation.The implicit Crank-Nicolson technique was then applied to the resulting approximation,and the subsequent equation was split into the forward and reverse sweeps,giving rise to the new axisymmetric ADE scheme.The pore pressure solutions obtained from the new scheme were then sequentially coupled with an existing geomechanical simulator in FLAC.This coupling procedure is called the sequentially-explicit coupling technique based on the fourth-order ADE scheme or SEA-4-AXI.By means of von Neumann and truncation error analyses,the new axisymmetric ADE scheme has proven to be unconditionally stable,consistent and fourth-order accurate in space.To demonstrate its application to geotechnical engineering,SEA-4-AXI was used to simulate the H-M interaction of a circular footing and of an advancing tunnel in deep saturated ground.

根据司法行政工作性质和队伍年龄、学历、职务等实际情况，河南司法行政系统教育培训工作坚持分级管理、分类实施，确保所属处级领导干部每5年参加培训累计3个月或者550学时以上、其他人员每年累计不少于12天或者90学时的培训，不同类别干部年度调训率、参训率达到100%，网络培训合格率达到99%，着力提升队伍的综合素质和能力，为推动司法行政工作改革创新提供人才保证和智力支持。

（3）灯光诱杀：每20亩悬挂频振式杀光灯一盏，诱杀小菜蛾、甘蓝夜蛾、斜纹夜蛾等蔬菜和金龟子、蝼蛄等地下害虫。

Because it used large f l ow time steps,SEA-4-AXI has resulted in more computationally efficient H-M simulations than the fully coupled approach using FLAC’s basic f l ow scheme while still maintaining numerical stability and high numerical accuracy.Results from the axisymmetric consolidation of a circular footing and of an advancing tunnel in deep saturated ground showed that SEA-4-AXI reduced computer runtime up to 42%-50%that of FLAC’s basic scheme.Yet,the accuracy of the H-M response was well maintained as shown by the average percentage difference of only 0.5%-1.8%.This increased computational efficiency and high-order accuracy,together with its capability in solving axisymmetric H-M problems in non-uniform domains,suggest that SEA-4-AXI holds great promise for solving larger coupled problems in geotechnical engineering.

Conflicts of interest

The authors wish to confirm that there are no known Conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

The authors gratefully acknowledge the support from the University Transportation Center for Underground Transportation Infrastructure at the Colorado School of Mines for partially funding this research under Grant No.69A3551747118 of the Fixing America’s Surface Transportation Act(FAST Act)of U.S.DoT FY2016.The authors also wish to recognize the support from the Center for Underground Construction and Tunneling at the Colorado School of Mines for providing the FLAC key for the simulations performed in this paper.

Appendix

The grid size-dependent constants for Eqs.(33)and(34)are given in Eqs.(A1)-(A5):

where

The constantα0is defined as follows:

改革开放初期，中央就多次强调，要大力发展我国食品工业建设。但食品工业的发展并非一蹴而就，起点低、技术薄弱、缺乏创新、内动力不足等问题成为改革过程中的主要制约因素，食品工业发展依旧缓慢。

The grid size-dependent constants for Eqs.(37)and(38)are given in Eq.(A6):

References

Abousleiman Y,Cheng AHD,Cui L,Detournay E,Roegiers JC.Mandel’s problem revisited.Géotechnique 1996;46(2):187-95.

Barakat HZ,Clark JA.On the solution of the diffusion equations by numerical methods.Journal of Heat Transfer 1966;88(4):421-7.

Biot MA.General theory of three-dimensional consolidation.Journal of Applied Physics 1941;12(2):155-64.

Bobet A.Characteristic curves for deep circular tunnels in poroplastic rock.Rock Mechanics and Rock Engineering 2010;43(2):185-200.

Callari C,Casini S.Tunnels in saturated elasto-plastic soils:three-dimensional validation of a plane simulation procedure.In:Frémond M,Marceri F,editors.Mechanical modelling and computational issues in civil engineering.Berlin:Springer;2005.p.143-64.

Chapra S,Canale R.Numerical methods for engineers.6th ed.New York:McGraw-Hill;2009.

Crank J,Nicolson P.A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type.Mathematical Proceedings of the Cambridge Philosophical Society 1947;43(1):50-67.

Cryer CW.A comparison of the three-dimensional consolidation theories of Biot and Terzaghi.Quarterly Journal of Mechanics and Applied Mathematics 1963;16(4):401-12.

Dean RH,Gai X,Stone CM,Minkoff SE.A comparison of techniques for coupling porous f l ow and geomechanics.Society of Petroleum Engineers Journal 2006;11(1):132-40.

Giraud A,Rousset G.Time-dependent behaviour of deep clays.Engineering Geology 1996;41(1-4):181-95.

Graziani A,Boldini D.influence of hydro-mechanical coupling on tunnel response in clays.Journal of Geotechnical and Geoenvironmental Engineering 2012;138(3):415-8.

Gutierrez M,Lewis R.Coupling of fluid f l ow and deformation in underground formations.Journal of Engineering Mechanics 2002;128(7):779-87.

Hoffman JD.Numerical methods for engineers and scientists.2nd ed.New York:McGraw-Hill;2001.

Itasca.Fast Lagrangian analysis of continua(version 7.00):User’s guide.Minneapolis,USA:Itasca Consulting Group;2011a.

Itasca.Fast Lagrangian analysis of continua:fluid/mechanical interaction.Minneapolis,USA:Itasca Consulting Group;2011b.

Itasca.Fast Lagrangian analysis of continua:FISH in FLAC.Minneapolis,USA:Itasca Consulting Group;2011c.

Li X.Stress and displacement fields around a deep circular tunnel with partial sealing.Computers and Geotechnics 1999;24(2):125-40.

Mandel J.Consolidation des sols(Étude mathématique).Géotechnique 1953;3(7):287-99(in French).

Mazumder S.Numerical methods for partial differential equations:finite difference and finite volume methods.Oxford,UK:Academic Press;2016.

McNamee J,Gibson RE.Plane strain and axially symmetric problems of the consolidation of a semi-infinite clay stratum.Quarterly Journal of Mechanics and Applied Mathematics 1960;13(3):210-27.

Minkoff SE,Stone CM,Bryant S,Peszynska M,Wheeler MF.Coupled fluid f l ow and geomechanical deformation modeling.Journal of Petroleum Science and Engineering 2003;38(1-2):37-56.

Neuzil CE.Hydromechanical coupling in geologic processes.Hydrogeology Journal 2003;11(1):41-83.

Prassetyo SH,Gutierrez M.efficient sequential coupling technique for the simulation of hydro-mechanical interaction in rock engineering.In:Proceedings of the 51st U.S.rock mechanics/geomechanics symposium.Alexandria,USA:American Rock Mechanics Association;2017a.p.1-11.

Prassetyo SH,Gutierrez M.Explicit higher-order ADE solutions for fluid f l ow in the coupled Biot equations.In:Poromechanics VI-proceedings of the 6th biot conference on poromechanics.Reston,USA:American Society of Civil Engineers;2017b.

Prassetyo SH,Gutierrez M.Effect of surface loading on the hydro-mechanical response of a tunnel in saturated ground.Underground Space 2016a;1(1):1-19.

Prassetyo SH,Gutierrez M.influence of embankment Loading on the hydromechanical response of a NATM tunnel in saturated ground.In:Proceedings of the 50th U.S.Rock mechanics/geomechanics symposium.Alexandria,USA:American Rock Mechanics Association;2016b.p.1-8.

Richtmyer RD,Morton KW.Difference methods for initial-value problems.2nd ed.New York:Interscience Tracts in Pure and Applied Mathematics;1967.

Strikwerda JC.Finite difference schemes and partial differential equations.2nd ed.Philadelphia,USA:Society for Industrial and Applied Mathematics;2004.

Wang HF.Theory of linear poroelasticity with applications to geomechanics and hydrogeology.Princeton,USA:Princeton University Press;2000.