Zienkiewicz_TheFiniteElement_Volume 2\50559_02

2Solutionofnon-linearalgebraicequations2.1IntroductionInthesolutionoflinearproblemsbyafiniteelementmethodwealwaysneedtosolveasetofsimultaneousalgebraicequationsoftheformKa=f(2.1)Providedthecoefficientmatrixisnon-singularthesolutiontotheseequationsisunique.Inthesolutionofnon-linearproblemswewillalwaysobtainasetofalgebraicequations;however,theygenerallywillbenon-linear,whichweindicateas*(a)=f-P(a)=0whereaisthesetofdiscretizationparameters,favectorwhichisindependentoftheparametersandPavectordependentontheparameters.Theseequationsmayhavemultiplesolutions[i.e.morethanonesetofamaysatisfyEq.(2.2)].Thus,ifasolutionisachieveditmaynotnecessarilybethesolutionsought.Physicalinsightintothenatureoftheproblemand,usually,small-stepincrementalapproachesfromknownsolutionsareessentialtoobtainrealisticanswers.Suchincrementsareindeedalwaysrequirediftheconstitutivelawrelatingstressandstrainchangesispathdepen-dentoriftheload-displacementpathhasbifurcationsormultiplebranchesatcertainloadlevels.Thegeneralproblemshouldalwaysbeformulatedasthesolutionof*,,+I=*(%+l)=f,,+l-P(a,,+1)=o(2.3)whichstartsfromanearbysolutionata=a,,,*,,=o.f=f,,(2.4)andoftenarisesfromchangesintheforcingfunctionf,,tof,,iI=f,,+Af,,(2.5)ThedeterminationofthechangeAa,suchthata,,+I=a,,+4,Iterativetechniques23Fig.2.1PossibilityofmultiplesolutionswillbetheobjectiveandgenerallytheincrementsofAf,willbekeptreasonablysmallsothatpathdependencecanbefollowed.Further,suchincrementalprocedureswillbeusefulinavoidingexcessivenumbersofiterationsandinfollowingthephysicallycorrectpath.InFig.2.1weshowatypicalnon-uniquenesswhichmayoccurifthefunction+decreasesandsubsequentlyincreasesastheparameterauniformlyincreases.ItisclearthattofollowthepathAf,?willhavebothpositiveandnegativesignsduringacompletecomputationprocess.Itispossibletoobtainsolutionsinasingleincrementoffonlyinthecaseofmildnon-linearity(andnopathdependence),thatis,withf,,=0,Af,,=f,,+1=f(2.7)Theliteratureongeneralsolutionapproachesandonparticularapplicationsisextensiveand,inasinglechapter,itisnotpossibletoencompassfullyallthevariantswhichhavebeenintroduced.However,weshallattempttogiveacomprehensivepicturebyoutliningfirstthegenerulsolutionprocedures.Inlaterchaptersweshallfocusonproceduresassociatedwithrate-independentmaterialnon-linearity(plasticity),rate-dependentmaterialnon-linearity(creepandvisco-plasticity),somenon-linearfieldproblems,largedisplacmentsandotherspecialexamples.2.2Iterativetechniques2.2.1GeneralremarksThesolutionoftheproblemposedbyEqs(2.3)-(2.6)cannotbeapproacheddirectlyandsomeformofiterationwillalwaysberequired.Weshallconcentratehereonproceduresinwhichrepeatedsolutionoflinearequations(i.e.iteration)oftheformK'da:,=r;?+I(2.8)24Solutionofnon-linearalgebraicequationsinwhichasuperscriptiindicatestheiterationnumber.Intheseasolutionincrementdabiscomputed.*GaussianeliminationtechniquesofthetypediscussedinVolume1canbeusedtosolvethelinearequationsassociatedwitheachiteration.However,theapplicationofaniterativesolutionmethodmayprovetobemoreeconomical,andinlaterchaptersweshallfrequentlyrefertosuchpossibilitiesalthoughtheyhavenotbeenfullyexplored.Manyoftheiterativetechniquescurrentlyusedtosolvenon-linearproblemsorigi-natedbyintuitiveapplicationofphysicalreasoning.However,eachofsuchtech-niqueshasadirectassociationwithmethodsinnumericalanalysis,andinwhatfollowsweshallusethenomenclaturegenerallyacceptedintextsonthissubject.’-’Althoughwestateeachalgorithmforasetofnon-linearalgebraicequations,weshallillustrateeachprocedurebyusingasinglescalarequation.This,thoughusefulfromapedagogicalviewpoint,isdangerousasconvergenceofproblemswithnumerousdegreesoffreedommaydepartfromthesimplepatterninasingleequation.2.2.2TheNewton-RaphsonmethodTheNewton-Raphsonmethodisthemostrapidlyconvergentprocessforsolutionsofproblemsinwhichonlyoneevaluationof*ismadeineachiteration.Ofcourse,thisassumesthattheinitialsolutioniswithinthezoneofattractionand,thus,divergencedoesnotoccur.Indeed,theNewton-Raphsonmethodistheonlyprocessdescribedhereinwhichtheasymptoticrateofconvergenceisquadratic.ThemethodissometimessimplycalledNewton’smethodbutitappearstohavebeensimultaneouslyderivedbyRaphson,andaninterestinghistoryofitsoriginsisgiveninreference6.Inthisiterativemethodwenotethat,tothefirstorder,Eq.(2.3)canbeapproxi-matedasHeretheiterationcounteriusuallystartsbyassumingat?,]I=a,,(2.10)inwhicha,isaconvergedsolutionatapreviousloadlevelortimestep.Thejacobianmatrix(orinstructuraltermsthestiffnessmatrix)correspondingtoatangentdirec-tionisgivenbydPa*KT=-=--(2.11)dadaEquation(2.9)givesimmediatelytheiterativecorrectionas*NotethedifferencebetwecnasolutionincrementdaandadifferentialdaIterativetechniques25Fig.2.2TheNewton-Raphsonmethod.ordal=(K~)-'!€$+,(2.12)Aseriesofsuccessiveapproximationsgivesi+l-ia,,,-an+]fdd.1=a,+Aai(2.13)whereiAai=cdat(2.14)k=lTheprocessisillustratedinFig.2.2andshowstheveryrapidconvergencethatcanbeachieved.TheneedfortheintroductionofthetotalincrementAa:,isperhapsnotobviousherebutinfactitisessentialifthesolutionprocessispathdependent,asweshallseeinChapter3forsomenon-linearconstitutiveequationsofsolids.TheNewton-Raphsonprocess,despiteitsrapidconvergence,hassomenegativefeatures:1.anewK,matrixhastobecomputedateachiteration;2.ifdirectsolutionforEq.(2.12)isusedthematrixneedstobefactoredateachiteration;3.onsomeoccasionsthetangentmatrixissymmetricatasolutionstatebutunsym-metricotherwise(e.g.insomeschemesforintegratinglargerotationparameters'ornon-associatedplasticity).Inthesecasesanunsymmetricsolverisneededingeneral.Someofthesedrawbacksareabsentinalternativeprocedures,althoughgenerallythenaquadraticasymptoticrateofconvergenceislost.26Solutionofnon-linearalgebraicequations2.2.3ModifiedNewton-RaphsonmethodThismethodusesessentiallythesamealgorithmastheNewton-RaphsonprocessbutreplacesthevariablejacobianmatrixK;byaconstantapproximation.-K;MKT(2.15)givinginplaceofEq.(2.12),dat=KTI*L+~(2.16)Manypossiblechoicesexisthere.ForinstanceKTcanbechosenasthematrixcorrespondingtothefirstiterationKi[asshowninFig.2.3(a)]ormayevenbeonecorrespondingtosomeprevioustimesteporloadincrementKO[asshowninFig.2.3(b)].Inthecontextofsolvingproblemsinsolidmechanicsthemethodisalsoknownasthestresstransferorinitialstressmethod.Alternatively,theapproximationcanbechoseneveryfewiterationsasKT=KiwherejI,nowdroppingsubscripts,da'=(K;)-'w(2.20)where(K;)-'isdeterminedsothatda1-l=(~;)-'(*f-l-=(K;)-'~I-'(2.21)Iterativetechniques27Fig.2.3ThemodifiedNewton-Raphsonmethod:(a)withinitialtangentinincrement;(b)withinitialproblemtangent.ForthescalarsystemillustratedinFig.2.4thedeterminationofK:istrivialand,asshown,theconvergenceismuchmorerapidthaninthemodifiedNewton-Raphsonprocess(generallyasuper-linearasymptoticconvergencerateisachievedforanormoftheresidual).ForsystemswithmorethanonedegreeoffreedomthedeterminationofKioritsinverseismoredifficultandisnotunique.ManydifferentformsofthematrixKfcansatisfyrelation(2.1)and,asexpected,manyalternativesareusedinpractice.AlloftheseusesomeformofupdatingofapreviouslydeterminedmatrixorofitsinverseinamannerthatsatisfiesidenticallyEq.(2.21).Somesuchupdatespreservethematrixsymmetrywhereasothersdonot.Anyofthemethodswhichbeginwith28Solutionofnon-linearalgebraicequationsFig.2.4ThesecantmethodstartingfromaKOpredictionasymmetrictangentcanavoidthedifficultyofnon-symmetricmatrixformsthatariseintheNewton-RaphsonprocessandyetachieveafasterconvergencethanispossibleinthemodifiedNewton-Raphsonprocedures.SuchsecantupdatemethodsappeartostemfromideasintroducedfirstbyDavidon8anddevelopedlaterbyothers.DennisandMore'surveythefieldexten-sively,whileMatthiesandStrang"appeartobethefirsttousetheproceduresinthefiniteelementcontext.Furtherworkandassessmentoftheperformanceofvariousupdateproceduresisavailableinreferences11-14,TheBFGSupdate'(namedafterBroyden,Fletcher,GoldfarbandShanno)andtheDFPupdate'(Davidon,FletcherandPowell)preservematrixsymmetryandpositivedefinitenessandbotharewidelyused.WesummarizebelowastepoftheBFGSupdatefortheinverse,whichcanbewrittenas(K')-'=(I+~,~~)(K,-')-'(I+~,~~)(2.22)whereIisanidentitymatrixandv,=1-(da'-IlTy'-I]*l-1[-*Id(a')TqL-'(2.23)1&-Iw,=da(i-l)TY1-1whereyisdefinedbyEq.(2.21).SomealgebrawillreadilyverifythatsubstitutionofEqs(2.22)and(2.23)intoEq.(2.21)resultsinanidentity.Further,theformofEq.(2.22)guaranteespreservationofthesymmetryoftheoriginalmatrix.Thenatureoftheupdatedoesnotpreserveanysparsityintheoriginalmatrix.Forthisreasonitisconvenientateveryiterationtoreturntotheoriginal(sparse)matrixKA,usedinthefirstiterationandtoreapplythemultiplicationofEq.(2.22)throughIterativetechniques29Fig.2.5Direct(orPicard)iteration.allpreviousiterations.Thisgivesthealgorithmintheformbl=fi(1-tvjw~)!@'.I=2b2=(Kf)-'h,(2.24)i-2da'=n(I+w~-~v~-~)~~T/=0Thisnecessitatesthestorageofthevectorsviandwjforallpreviousiterationsandtheirsuccessivemultiplications.Furtherdetailsontheoperationsaredescribedwellinreference10.Whenthenumberofiterationsislarge(i>15)theefficiencyoftheupdatedecreasesasaresultofincipientinstability.Variousproceduresareopenatthisstage,themosteffectivebeingtherecomputationandfactorizationofatangentmatrixatthecurrentsolutionestimateandrestartingtheprocessagain.AnotherpossibilityistodisregardallthepreviousupdatesandreturntotheoriginalmatrixKB.SuchaprocedurewasfirstsuggestedbyCri~field'l.'~.'~in'thefiniteelementcontextandisillustratedinFig.2.5.Itisseentobeconvergentataslightlyslowerratebutavoidstotallythestabilitydifficultiespreviouslyencounteredandreducesthestorageandnumberofoperationsneeded.Obviouslyanyofthesecantupdatemethodscanbeusedhere.TheprocedureofFig.2.5isidenticaltothatgenerallyknownasdirect(orPicard)iteration'andisparticularlyusefulinthesolutionofnon-linearproblemswhichcanbewrittenas@(a)=f-K(a)a=0(2.25)Insuchacasea:+1=a,istakenandtheiterationproceedsasI+1Ia,,.,=[K(a:+I)]-frt+l(2.26)30Solutionofnon-linearalgebraicequations2.2.5Linesearchprocedures-accelerationofconvergenceAlltheiterativemethodsoftheprecedingsectionhaveanidenticalstructuredescribedbyEqs(2.12)-(2.14)inwhichvariousapproximationstotheNewtonmatrixKfareused.Foralloftheseaniterativevectorisdeterminedandthenewvalueoftheunknownsfoundasa:,:',=ai,l+da:,(2.27)startingfromI%+I=a,inwhicha,,istheknown(converged)solutionattheprevioustimesteporloadlevel.Theobjectiveistoachievethereductionof!PL:l,tozero,althoughthisisnotalwayseasilyachievedbyanyoftheproceduresdescribedeveninthescalarexampleillustrated.Togetasolutionapproximatelysatisfyingsuchascalarnon-linearproblemwouldhavebeeninfacteasierbysimplyevaluatingthescalarq:L'lforvariousvaluesofa,+landbysuitableinterpolationarrivingattherequiredanswer.Formulti-degree-of-freedomsystemssuchanapproachisobviouslynotpossibleunlesssomescalarnormoftheresidualisconsidered.Onepossibleapproachistowritel+l,/-ia,,.,-an+,+v,.,dai(2.28)anddeterminethestepsizevi,,sothataprojectionoftheresidualonthesearchdirec-tiondaiismadezero.WecoulddefinethisprojectionasG,,~=(daL)TqiA1,J(2.29)where$,>'i'=*(ai+1+v,,,da;),vj,o=1Here,ofcourse,othernormsoftheresidualcouldbeused.Thisprocessisknownasalinesearch,andvi,,canconvenientlybeobtainedbyusingaregulafulsi(orsecant)procedureasillustratedinFig2.6.AnobviousFig.2.6Regulafahappliedtolinesearch:(a)extrapolation;(b)interpolation.Iterativetechniques31disadvantageofalinesearchistheneedforseveralevaluationsof@.However,theaccelerationoftheoverallconvergencecanberemarkablewhenappliedtomodifiedorquasi-Newtonmethods.Indeed,linesearchisalsousefulinthefullNewtonmethodbymakingtheradiusofattractionlarger.Acompromisefrequentlyused"istoundertakethesearchonlyifG,>E(dab)TQ:,++I~I(2.30)wherethetoleranceEissetbetween0.5and0.8.ThismeansthatiftheiterationprocessdirectlyresultedinareductionoftheresidualtoEorlessofitsoriginalvaluealinesearchisnotused.2.2.6'Softening'behaviouranddisplacementcontrolInapplyingtheprecedingtoloadcontrolproblemswehaveimplicitlyassumedthattheiterationisassociatedwithpositiveincrementsoftheforcingvector,f,inEq.(2.5).Insomestructuralproblemsthisisasetofloadsthatcanbeassumedtobepropor-tionaltoeachother,sothatonecanwriteAf,,=AA,f,(2.31)Inmanyproblemsthesituationwillarisethatnosolutionexistsaboveacertainmax-imumvalueoffandthattherealsolutionisa'softening'branch,asshowninFig.2.1.InsuchcasesAx,willneedtobenegativeunlesstheproblemcanberecastasoneinwhichtheforcingcanbeappliedbydisplacementcontrol.Inasimplecaseofasingleloaditiseasytorecastthegeneralformulationtoincrementsofasingleprescribeddisplacementandmuchefforthasgoneintosuch~olutions.".'~-~'InallthesuccessfulapproachesofincrementationofAA,theoriginalproblemofEq.(2.3)isrewrittenasthesolutionofwith(2.32)beingincludedasvariablesinanyincrement.Nowanadditionalequation(constraint)needstobeprovidedtosolvefortheextravariableAx,,.Thisadditionalequationcantakevariousforms.Riksl'assumesthatineachincrementAaTAa,+AA2fifo=AI'(2.33)whereAIisaprescribed'length'inthespaceofn+1dimensions.Crisfieldl'.24pro-videsamorenaturalcontrolondisplacements,requiringthatAa;fAa,,=(2.34)Theseso-calledarc-lengthandsphericalpathcontrolsarebutsomeofthepossibleconstraints.32Solutionofnon-linearalgebraicequationsDirectadditionoftheconstraintEqs(2.33)or(2.34)tothesystemofEqs(2.32)isnowpossibleandthepreviouslydescribediterativemethodscouldagainbeused.However,the'tangent'equationsystemwouldalwaysloseitssymmetrysoanalter-nativeprocedureisgenerallyused.WenotethatforagiveniterationiwecanwritequitegenerallythesolutionasQ;+,=X;+,fo-P(aL+])(2.35)Qi'+'i=QL+I+dX;,fo-Kids;ThesolutionincrementforamaynowbegivenasdaL=(K\)-'[QL++dXLfo](2.36)dak=da:,+dXLdaiwhere(2.37)Nowanadditionalequationiscastusingtheconstraint.Thus,forinstance,withEq.(2.34)wehaveT(Ask-+dai)(Ask-+dak)=A12(2.38)whereAsh-'isdefinedbyEq.(2.14).OnsubstitutionofEq.(2.36)intoEq.(2.38)aquadraticequationisavailableforthesolutionoftheremainingunknowndXk(whichmaywellturnouttobenegative).Additionaldetailsmaybefoundinreferences11and24.AproceduresuggestedbyBergan2023issomewhatdifferentfromthosejustdescribed.HereafixedloadincrementAx,isfirstassumedandanyofthepreviouslyintroducediterativeproceduresareusedforcalculatingtheincrementdak.NowanewincrementAx,*,iscalculatedsothatitminimizesanormoftheresidual[(AX:,f0-P;;")'(nX;fo-P;;")]=A12(2.39)Theresultisthuscomputedfrom-=odA12dAX;andyieldsthesolution(2.40)Thisquantitymayagainwellbenegative,requiringaloaddecrease,anditindeedresultsinarapidresidualreductioninallcases,butprecisecontrolofdisplacementmagnitudesbecomesmoredifficult.TheinterpretationoftheBerganmethodinaone-dimensionalexample,showninFig.2.7,isilluminating.Hereitgivestheexactanswers-withadisplacementcontrol,themagnitudeofwhichisdeterminedbytheinitialAX,assumedtobetheslopeKTusedinthefirstiteration.Iterativetechniques33Fig.2.7One-dimensionalinterpretationoftheBerganprocedure.2.2.7ConvergencecriteriaInalltheiterativeprocessesdescribedthenumericalsolutionisonlyapproximatelyachievedandsometolerancelimitshavetobesettoterminatetheiteration.Sincefiniteprecisionarithmeticisusedinallcomputercalculations,onecanneverachieveabettersolutionthantheround-offlimitofthecalculations.Frequently,thecriteriausedinvolveanormofthedisplacementparameterchangesIIda:]Ior,morelogically,thatoftheresidualsIlQ:?+II1.InthelattercasethelimitcanoftenbeexpressedassometoleranceofthenormofforcesIlfn+lll.Thus,wemayrequirethatIlQL+llldEllfn+lll(2.41)whereEischosenasasmallnumber,andj(Q(l=(QTQ)I’*(2.42)Otheralternativesexistforchoosingthecomparisonnorm,andanotheroptionistousetheresidualofthefirstiterationasabasis.Thus,IlQf+llldEllQ:+lll(2.43)Theerrorduetotheincompletesolutionofthediscretenon-linearequationsisofcourseadditivetotheerrorofthediscretizationthatwefrequentlymeasureintheenergynorm(seeChapter14ofVolume1).Itispossiblethereforetousethesamenormforboundingoftheiterationprocess.Wecould,asathirdoption,requirethattheerrorintheenergynormsatisfydE‘=da>~lBb+l)”2