Numerical Analysis of Simulation

ORIGINALPAPERNumericalAnalysisofSimulationAccuracyforSeepageChangeinOne-DimensionalLayerUnderLow-FrequencyWaveLimingZheng.ZifengLiReceived:14March2017/Accepted:16August2017�SpringerInternationalPublishingAG2017AbstractTheparameteranalysisofpetrophysicalpropertiescouldbebeneficialforexplainingtheactionmechanismsofwavestimulationtechnique,whichcoupledtheoriginalDarcyflowfieldwiththelow-frequencywavefield.SpuriousreflectionoccurringattheartificialboundariesingeneralBiotconsolidationproblemcausedbythewavefieldwouldaffectthesimulationaccuracy.Theconsistencerequirementfordirectionsofwavepropagation,flowvelocity,andsoilspreadinginone-dimensionalmodelhadincreasedthedifficultyforobtainingreasonableresults.StrategiesasadditionofPMLs,modelsizeoptimization,approximationorderadjustmentofshapefunctions,andpropermeshgenerationwereinvestigatedinsequencewiththesystemCOMSOLtoimprovetheaccuracyofresults.Severalschemeswerecomparedineachstrategy.Periodicvariationovertimeandsharpdeclineoverdistanceofpetrophysicalpropertiesandwereobserved,whichagreewithgeneralexper-imentalcognitionofthewavestimulationtechnique.Thegridshapewouldinfluencesimulationaccuracyandabetterdistributionofporepressurewasfoundwithquadrilateralmeshgenerationintheone-dimen-sionalphysicalmodel.Consistenceinx-coordinateofnodescontributedtoaparallelfluctuationofproper-ties.Acomparisonbetweenthenumericalresultsandapproximateanalyticalsolutionswasmadeattheend,whichhadverifiedtheadvantageofnumericalsolu-tiononmechanismexplanation.Furtheranalysisforchoosingproperoperatingparametersofwavestim-ulationtechniquemightbehelped.KeywordsLow-frequencywave�Originalseepage�Finiteelement�Accuracy�Periodicvariation1IntroductionBiotporoelasticproblemwasinvolvedinmanygeologicalfieldssuchascivilengineering,oceanengineering,petroleumengineering,geophysicalprospecting,andgeothermalsystem.Itdenotedacouplingbetweensolidmechanicsandfluidflowinporousmedia.Thoughmanyresearchadvanceswereachieved,solutionofBiotporoelasticmodelinnumericalmethodwasstillachallengingworkduetothecouplingandinconsistencyofmultiphysicalfields.Instabilityandnon-convergenceofresultsappearedinmanypapers,especiallyaroundthedegreesofsoildeformationandfluidpercolationL.Zheng(&)�Z.LiCollegeofVehicleandEnergy,YanshanUniversity,Qinhuangdao066004,Chinae-mail:upczlm@sina.cnL.ZhengCollegeofPetroleumEngineering,ChinaUniversityofPetroleum(EastChina),Qingdao266580,Shandong,China123GeotechGeolEngDOI10.1007/s10706-017-0343-4Administrator高亮Administrator高亮Administrator高亮enhancement.Solutionanalysisandstabilityimprove-mentbecamethefocusofsomeresearchers.IntheBiotporoelasticproblemforthemechanismsofwavestimulationtechniquesinoilfielddevelop-ment,thephysicalpropertiessuchassoildeformation,flowvelocity,andporepressurewerealsoinfluencedbytheelasticwave.Changesinflowvelocityandporepressuredecidedtheoptimizationofvibrationfre-quencyandamplitude.DifferentfromthatingeneralBiotporoelasticproblem,thewavepropagationwasofsamedirectionwiththoseofseepageandsoilspread-inginwavestimulationtechnique.Whenanumericalsimulationwascarried,one-dimensional(1D)modelofporousmediawasbuiltinsystemCOMSOLtoverifytheexperimentaldataoranalyticalsolution.Theporousmediafullysaturatedwithfluidwasassumedwithacertainvelocityinitiallyandstimu-latedbylow-frequencywaveatoneedge.Stabilityandconvergenceofcalculationwereseriouslyinfluencedinactualsolutionbecauseofthespuriousreflectionsoccurringattheartificialboundaries.Strategiesforaccuracyincreasingwereessentialtofurtherlyprovideaninsighttotheactualmechanisms.ThestrategiescanrefertotheprogressingeneralBiotmodelling.Thereasoncausingnon-convergenceofsolutiontotheBiotpartialdifferentialequationswasconcludedbyFerronatoetal.(2010).Forexamples,hugenumberofunknownsofalgebraicsystemsinfullycoupledapproach,severelyill-conditionedcoefficientmatrixinthenumericaldiscretization,anddifferentformsofinstabilitiessufferedbythenumericalsolution.ThoughincompletesummarizedbyFerronato,acomprehensiononthestabilitycausedbyalgorithminBiotporoelasticproblemwasgiven.Thephysicalmodelbuilding,initialvalue,andmeshgenerationwouldalsoaffectthestabilityinFEsimulation.Certainly,thereexistedseveralcoupledapproaches(fully,iteratively,explicitly,andlooselycoupledmethods)insequentialstrategiesforsolution,andeachhaditsadvantagesanddisadvantages.Thefullycoupledmethod,inwhichthegoverningequationsofflowandsolidweresolvedsimultaneouslyateachtimestep,wasunconditionallystableandeasilyusedwhentheelementsnumberwassmall.Intheiterativecoupledmethod,oneofthegoverningequationswassolvedfirstandthentheotherwassolvedusingtheintermediatesolutioninformation.Explicitlycoupledmethodwasaspecialcaseofiterativecoupledmethodwithconsideringiterationforonlyonetime.Thelooselycoupledmethodwasofthefunctionasvariedgridorstepsizetosomeextent,oneequationwassolvedonlyaftermultiplestepsoftheotherequationweretaken.Rozaetal.(2014)performedacompre-hensivestudyonthestability,accuracy,andrateofconvergenceofiteratively,explicitly,andlooselycoupledmethodsforthesolutionofaBiotconsolida-tionmodel.Phametal.(2016)appliedthestaggerediterationprocessasapartitionedoriterativeapproachtosolvetheill-conditionedlinearequationsystems.Luoetal.(2016)employedaniterativelycouplingmethodandparallelcomputingtofullycaptureinteractionsbetweensolidandflowandhandlelargescaleproblems.Inthispaper,thesolverinCOMSOLsoftwarewassettoautomatic,whichcouldchoosethedirectoriterativesolverautomatically(KluczykandJacak2016).Forthesmallnumberofcalculationnodesinsimplephysicalmodel,theMUltifrontalMassivelyParallelSparsedirectSolver(MUMPS)belongingtofullycoupledmethodwaspriorconsidered.ThecoefficientmatrixofvariablesinalgorithmanalysisofBiotporoelasticproblem,whichmaybesparseandindefinite,wouldinfluencetherobustnessofsolution.Whendecomposingtheoriginalcoeffi-cientmatrixinFEsimulation,researchershadimprovedthestabilitythroughrevisingthematrix.Ferronatoetal.(2010)transformedthesparse,sym-metricandindefinitematrixofdifferentialvariablecoefficientsindiscretesolutionschemeintoblock-diagonalsymmetricpositivedefinitematrix.Badiaetal.(2009)adoptedbothmonolithicsolversandheterogeneousdomaindecompositionstrategiesinsolutionofN-SequationswiththeBiotsystem.WhiteandBorja(2008)introducedastabilizingsub-matrixintothe(2,2)positionofcoefficientmatrixinFEweakform,eliminatingthezero-block.Bothadditionofstabilityfunctionandomittingtermsinthematrixmightexpandthezero-blocktoimprovethecompu-tationefficiency.Thecoefficientmatrixinthispaperwasnotrevisedbecauseofthedirectsolverusedforsmallnumberofcalculationnodes.IndifferentcasesofgeneralBiotporoelasticproblem,theinteractionbetweensolidmechanicandwave-inducedflowwasvarious.Toadapttotheinteraction,directionofsoilspreadingmightbeparallelorperpendiculartothatofwavestimulationatoneendinthephysicalmodelwithadrainageorun-drainageboundaryattheoppositeend.AfterbuildingGeotechGeolEng123Administrator高亮Administrator高亮Administrator高亮Administrator高亮Administrator高亮Administrator高亮Administrator高亮thephysicalmodel,animportantprocedurewastosetupthegridsizeandmeshing.Freetriangular(orquadrilateral)gridin2Dmodelandfreetetrahedral(orhexahedrongrid)in3Dmodelwasgenerated.Forcasesthatcommandedforvariousspatialortemporalgridsizesfordifferentvariable,spatialmeshgener-ationcouldbeachievedinautonomyprogramming.Perfectlymatchedlayer(PML)wasintroducedinamountsofarticlesforwaveabsorption,reflectionandtransmission.Meshrefinement,perfectlymatchedlayer,andepitaxyofresearchdomain,etc.,werehelpfulforaccuracyimprovement.Forexample,Gessner(2009)usedstructuredmeshescontainingbrickshapedelementsinthesimulationofsoliddeformationandfluidflowrelatedtoMesoproterozoiccoppermineralization.TwodifferentspatialgridswereusedbyRozaetal.(2014)todiscretizethepressureandthedisplacementequations,respectively.KrishnamoorthyandKamal(2016)investigatednumericallytheeffectofpartialdrainageduringpenetrationwithsquaregridsonthemeasuredtipresistanceandthesubsequentporepressuredissipa-tionresponse.Zhaoetal.(2007)studiedtheperfor-manceofPMLforbodywaveswithvariousincidentanglesandvariousabsorbingthicknesses,aswellastheabsorptionoffreesurfaceRayleighwaves.NumericalresultsofKucukcobanandKallivokas(2013)demonstratedthestability,efficacy,andcost-effectivenessofthehybridformulation,withPMLsintroducedattheheterogeneitydomaintruncationsurface.Kimetal.(2014)carriedouttwosimulationcaseswithandwithoutPMLstocomparativelyevaluatetheeffectivenessofthePMLtechnique.Aparticularsizeandmeshgenerationwereanalyzedinthispaperconsideringthespecificdirectionsandboundariesforsimulationofwavestimulationtech-nique.WiththesystemCOMSOL,spatialmeshgenerationwasimprovedtosomeextentbyvaryingtheapproximateorderofshapefunctionsinthispaper.Sincethecalculatingtimewaslittleandsimilartothetimescaleofwavepropagation,temporalmeshgenerationwasnotconsideredhere.Whenthecalcu-latingtimeapproximatedthatofoilfielddevelopment,temporalmeshgenerationshouldnotbeneglected.Besidesmeshrefinement,bothlowrelativetoler-anceandhighapproximateorderofshapefunctionscouldincreasetheaccuracyofresultsbutdecreasethecalculationefficiency.Therelativetolerancecanbesetupwithalowervalueasmuchaspossible.Asmentionedabove,differentapproximateorderofshapefunctionscouldbeadoptedtochangethespatialgridsizeslightly,whichwasalsousedinmanypapers.Forexample,Reed(1984)observedthatusingadifferentapproximationorderfordisplacementandpressuremayhelpkeepthespuriousoscillationsundercontrol.ChungandChiou(2002)researchedtheapplicationsofSelf-AdjustedConvexApproximationinstructuralshapeoptimizationproblems.Shapefunctionsforamixedp-versionfiniteelementformu-lationbasedonaweakformofHamilton’sextendedprincipleinthetimedomainweredevelopedbyHodgesandHou(1991),producingahighpercentageofzerosincertaincoefficientmatrix.Serraetal.(2015)investigatedtheinfluenceofboundarycondi-tionsonwavepropagationinBiot–Allard’smodelofporoelasticitywhichdependedonthespatialdis-cretization,theorderofshapefunctions,andthechoiceoftheformulation.Todenotethecouplingofseepagefieldandwavefieldinlow-frequencywavestimulationtechnique,aspecialphysicalmodelwasbuilt.Darcyflowattheinitialstagewasintroducedtodemonstratethedevel-opmentofreservoirandafluctuatingdisplacementboundaryconditionwasintroducedtodemonstratetheinfluenceofdownholelow-frequencyartificialwave.Anadditionalrequirementforaconsistenceinthedirectionsofwavepropagation,soilspreadingandfluidflowhadmadethesimulationcomplicated.Ononehand,theaccurateofcalculationduetomultifieldscouplinghadtobeincreased.Ontheotherhand,themechanismoflow-frequencywavestimulationtechniqueinreservoirdevelopingprocessshouldbereflectedcorrectly.Withseveralstrategiesintroducedabove,theaccuracyofgoverningequationsofBiotporoelasticproblemforaparticularapplicationinoilfielddevelopmentmightbeimproved.Thereby,numericalstudyonthesimulationaccuracywasconducted.Firstly,PMLswereaddedinphysicalmodelsoastoweakenthespuriousreflectionsfromboundariesandkeeptheconsistenceofdirectionsoforiginalseepagefieldandwavefieldatthebeginning.Then,theapproximateorderofshapefunctionswasadjustedtodecreasetheerror.Thirdly,themodelsizeinthecomputationaldomainandmeshgenerationmethodwereoptimizedinseveralschemestoincreasetheaccuracyofphysicalproperties.Finally,acom-parisonwiththeapproximateanalyticalsolutionwasmadetoverifytheeffectivenessofnumericalGeotechGeolEng123Administrator高亮Administrator高亮Administrator高亮Administrator高亮simulation.Periodicvariationsovertimeandsharplydeclineoverdistanceofphysicalpropertiesfrominjectionendtoproductionendweresimulated.Aboveanalysisweredeterminedtolayfunctiononfurthermechanismexplanationofwavestimulationtechniques.2ModelEquationsBasedontheinteractionofsolidmechanicandlow-frequencywave-inducedfluidflowinlowpermeabil-ityreservoir,thegoverningequationswerebuiltandpartiallyassumedasthatofBiot(1956).Thegovern-ingequationshadalsoincludedspecificboundaryconditionsandinitialconditionsconsideringconstantpressureproductioninoilfielddevelopment.ThefluidfollowedDarcy’slawattheinitialstage.Thefluidwasinjectedfromaninjectionwellwithaconstantpressureandflowedoutfromaproductionwellwithanotherconstantpressure.ThewavevelocitywasalittlereducedbytheinertiaeffectoforiginalDarcyflowduetothecouplingofmatrix,fluidfieldandwavefield.Thevariationofwavevelocitymightrefertotheexperimentsaboutsoilconsolidationunderdrainagecondition(StroiszandFjær2011).Inturn,theoverallflowvelocitywasstimulatedbytheBiotflow(Zhengetal.2015).Animportantresearchobjectoflow-frequencywavestimulationtechniquewasaroundthevariationofflowvelocityaswellasitssensitivityunderdifferentphysicalpropertiesandoperatingparameters.Theaccuracyandstabilityofresultswasconcernedwiththeproperexplanationtovariationofflowvelocity.Thegoverningequationsaswellasitsdiscretizationformweregivenbelowforthefollowingaccuracyanalysis.2.1GoverningEquationsForanisotropic,homogenousporousmediasaturatedwithsingle-phaseflowingfluid,thegoverningequa-tionwasshownbelow.Thephysicalpropertieswereinfluencedbyeachotherintheformationdevelopmentprocess.Equations(1)–(3)formedacoupledpartialdifferentialsystemaboutu–pmethod.Equation(3)wasthegeneralizedDarcyequation,withtheinertiaforcesofsolid–fluidunderwave.l/ðÞr2uþgradl/ðÞþk/ðÞþa2M/ðÞ��divu�þaM/ðÞ�divwgþb¼q€uþqf€wð1ÞootaM/ðÞ�divuþPf��þM/ðÞdiv_w¼fð2Þ�gradPf¼qf€uþ€w/��þgxðÞkrPf��_w�k0kðÞð3Þwhereuandwwas,respectively,thesolidandrelativefluiddisplacement;ewasthevolumestrainofmatrix;Pfwasthefluidpressure;bwasthebodyforceonsolid;fwastheflowsourceorsink;qf,qsandqwas,respectively,thedensityoffluid,solidandporousmedium,q=/qf?(1-/)qs;gwasthedynamicfluidviscosity;/wastheporosity;kwasthepermeability;l,k,a,andMwereLamecoefficientsandBiotcoefficients,k¼Kb�23l,l=G,a=1-Kb/Ks;Kfwasthebulkmodulusoffluid;Ks,Kb,andGwasthebulkmodulusofrockinjacketedcompressiontest,bulkmodulusofmatrixinjacketedcompressiontestandshearmodulusofmatrix,respectively;1/Qwasacoefficientrepresentingthecouplingrelationshipbetweenvolumechangesoffluidandsolid,1/Q=(a-/)/Ks;k0wasthequasi-thresholdpressuregradient;xwasthevibrationfrequency.Certainly,themotionequationobeyedthegeneralDarcy’slawinitiallyandthefluidinlowpermeabilitylayerhadtoovercomethequasi-thresholdpressuregradientasinEq.(4)whenwithoutconsideringvibration.�gradPfþk0kðÞ¼gkrPf��_wð4ÞNeglectingtheminorchangesofporosityandpressureovertime,thetotalcontinuityequationwasderivedasinEq.(5),atransformofEq.(2)consid-eringthecompressibilityofsolid–fluid.�gradPf¼aMr2uþMr2wð5ÞWhensolvingthecoupledpartialdifferentialsystemwithu-wmethod,whichisprovedeasierandoflargerconvergencethanu-pmethod(Wang2002),Eq.(6)wasgotcombiningtherightsidesofEqs.(3)and(5).GeotechGeolEng123gradaM/ðÞ�divuþM/ðÞdivw½�¼qf€uþ€w/��þgxðÞkrPf��_w�k0kðÞð6ÞEquations(1),(5),and(6)formedacoupledpartialdifferentialsystemwithu-wmethoddefinedon1DdomainXCboundedbythefrontierC.ThesystemwassolvedwithappropriateboundaryconditionsasinEq.(7),andinitialconditionsinEq.(8).ux;0ðÞ¼u0x;tðÞ;Pfx;tðÞ¼Pf0xðÞ_wx;0ðÞ¼_w0xðÞ;wx;tðÞ¼w0xðÞ ð8ÞInEqs.(7)and(8),CI[CP[Cv=C,theright-handsidesareknownfunctions.CIwastheinjectionboundary;CPwastheproductionboundary;Cvwastheotherboundariesverticaltosoilspreading.Thedisplacementsofsolidandfluidattimetwerederivedfirstlyinfiniteelement(FE)solutionofgoverningequationswithdefinitecondition.Then,thepressurechangeaswellastheadjustedphysicalpropertieswasgotafterthedisplacementsweresubstitutedintothetotalcontinuityequation.Third,newvaluesofdisplacementsandpressureatnexttimet?1withadifferentcoefficientmatrixweredonethroughiterativecalculation.2.2FEDiscretizationAweakformofgoverningequationsconsideringtheboundaryvalueproblemwaspresentforFEsolution.ThevariablesasdisplacementsandpressureweresetinshapefunctionsNku,Nkw,andNkPwithdenotingthenumberofelementsofnu,nw,andnp.ui¼Xnuk¼1Nuk~uki;wi¼Xnwk¼1Nwk~wki;Pfi¼XnPk¼1NPk~Pkið9ÞTheweakformwithintegrationofpartsofEqs.(1)and(6)waswrittenasEq.(10),whenthebodyforcewasneglectedduringtheprocessinnertheporousmedia.M€UþC_UþKU¼Fð10ÞwhereU¼~u~wðÞT,M¼M1M2M3M4��,C¼000C4��,K¼K1K2K3K4��,F¼F1F2��,M1¼qRANulNukdA,M2¼qfRANulNwkdA,M3¼MT2¼qfRANwlNukdA,M4¼qf/RANwlNwkdA,K1¼lRANul;jNuk;jdAþlþkþa2MðÞRANul;iNuk;jdA,K2¼aMRANul;iNwk;jdA,K3¼KT2¼aMRANwl;iNuk;jdA,K4¼MRANwl;iNwk;jdA,C4¼gkRANwlNwkdA,F1¼�HNulrnð�aPnÞdsþRANwlk0dA,F2¼HNwlPnds.Theanytime-dependentfunctionwasconsideredtovarylinearlyintimebetweentandt?Dt,thentheanytime-derivativeitematanintermediateinstants=h(t?Dt)?(1-h)twasapproximatedbyafirst-orderincrementalratio.hwasascalarvaluecomprisedbetween0and1.AlinearalgebraicmatrixasEq.(11)wascomputedafterthediscretesolutionschemewasconducted.ZUtþDt¼Fð11ÞwhereZ¼MþCDtþhKDtðÞ2hi;F¼hDtðÞ2FtþDtþ2MðþCDtÞUt�MUt�Dtþ1�hðÞDtðÞ2Ft�ðKUtÞ.2.3StateEquationandDiscretizationThephysicalpropertiesinthecoefficientmatrixofEq.(10)wererelatedwithotherparameters(Maetal.2011;Liuetal.2012;Zhengetal.2007;Zhou2003;Dengetal.2012).Then,thespatialortemporaldiscretizationformsofpropertieswerewrittenasEq.(12)insolutionprocess.ux;tðÞ¼uIx;tðÞ;_wx;tðÞ¼_wIx;tðÞ;Pfx;tðÞ¼PIx;tðÞoverCIrux;tðÞ¼ruPx;tðÞ;rwx;tðÞ¼rwPx;tðÞ;Pfx;tðÞ¼PPx;tðÞoverCPrux;tðÞ¼ruvx;tðÞ;rwx;tðÞ¼rwvx;tðÞ;rPfx;tðÞ¼rPvx;tðÞoverCv8><>:ð7ÞGeotechGeolEng123Ptþ1fi¼Ptfi�aMutþ1iþut�1i�2uti���þMwtþ1iþwt�1i�2wti���=Dlð12aÞ/tþ1i¼/tiþn�aetþ1i�eti��þ1QcPtþ1fi�Ptfi��;ktþ1i¼k0/tþ1i/0��31�/01�/tþ1i!2ð12bÞqtþ1fi¼qf0exþ1fi�Pf0Kf!;qtþ1si¼qs0expa�1ðÞetþ1iþ1QcPtþ1fi�Pf0�1�/tþ1i2435ð12cÞKtþ1Si¼KS0exp�nS�/tþ1i��;Ktþ1bi¼Kb0expð�nb�/tþ1iÞ;Gtþ1i¼G0expð�nG�/tþ1iÞð12dÞktþ10i¼A�ktþ1i���Bð12eÞwhere/0wastheinitialporosityinreservoir;k0wastheinitialpermeabilityinreservoir;KS0,Kb0,G0wereapproximatelythemodulusofmatrix,andnS,nb,nGwerethecoefficientsobtainedfromfittingcurveswithporosity;nwasacoefficientabouteffectiveporositychangeundervibrationcausedbychangesinheightofadhesivelayerandeffectiveporevolume,1.1–1.5;u0wasthevibrationamplitude.2.4CoefficientPDEsinCOMSOLAfterthephysicalmodelwasbuiltinsystemCOMSOL,thegoverningequationwasachievedbyformingthreecoefficientPDEsasEq.(13)withspecificconditionsasEqs.(14)–(16).eio2<ot2þdio<otþr��cir<�ai<þciðÞþbi�r<þai<¼fi;<¼u;w;Pf��Tð13ÞTheboundaryconditionandinitialconditionwereassumedasEqs.(14)–(16)consideringconstantpres-suredevelopmentinreservoir.Where,XCrepresentedthecomputationaldomaininFig.1a;Cirepresentedtheboundaryi;PinandPoutwastheconstantpressureatinjectionendandproductionend,respectively.Thewaveformoflowfrequencyvibrationinseismicproductiontechniquewassimplifiedasatrigonomet-ricfunction(Bei2009),whichcouldberewrittenasuser-definedasthevibrationsourcechanged.a.SolidPDE:ujXC¼0;_ujXC¼0;t¼0ujCI¼u0expixtðÞ;t�0 ð14Þb.FluidPDE:c.PressurePDE:Pf��XC&fra