邊界單元法理論及程序設(shè)計1-位勢邊界元

BoundaryElementMethods—TheoryandProgramming

Xiao-WeiGAO,KaiYang(高效偉,楊愷)大連理工大學(xué)航空航天學(xué)院邊界單元法理論及程序設(shè)計SchoolofAeronauticsandAstronauticsDalianUniversityofTechnology1.Introduction

1.1Whytolearnboundaryelementmethods(BEM)？為什么要學(xué)邊界單元法(BEM)？Tosolveproblemsinwhichfiniteelementmethod(FEM)isinadequateorinefficient.GeometryisnotregularComputationalregionisInfiniteorsemi-infiniteFracturetipanalysis引言

BEMiseasytogeneratethedatarequiredtorunaproblemandcarryoutthemodificationsneededtoachieveanoptimumdesign.Whytolearnboundaryelementmethods(BEM)？為什么要學(xué)邊界單元法(BEM)？OnlysurfacesoftheproblemneedtobediscretizedintoboundaryelementsMoreconvenienttodooptimumdesigncomputationSuitabeforparameterbackanalysisBEMismoreaccuratethanothernumericalmethods.Whytolearnboundaryelementmethods(BEM)？為什么要學(xué)邊界單元法(BEM)？Semi-analysischaracteristicUseoffundamentalsolutionsGradientsofphysicalquantitieshavethesameaccuracylevelasthephysicalquantitiesselfBEMisveryefficienttosolveinfiniteorsemi-infiniteproblems.Whytolearnboundaryelementmethods(BEM)？為什么要學(xué)邊界單元法(BEM)？InfiniteboundaryconditionscanbeautomaticallymodeledUseofInfiniteelementsUndergroundopeningPilefoundationproblemsAerodynamicproblemsBEMisrobusttosolvestress(orflux)concentrationproblems.Whytolearnboundaryelementmethods(BEM)？為什么要學(xué)邊界單元法(BEM)？FundamentalsolutionsaresingularfunctionsFractureproblemsRigidfoundationproblemsBoundarylayerproblemsinfluidmechanicsReferences（參考文獻）BoundaryElementsAnIntroductoryCourse，

BrebbiaCA,DominguezJ，ComputationalMechanicsPublications，1992。BoundaryElementProgramminginMechanics，

GaoXW&DaviesTG，CambridgeUniversityPress，2002。邊界單元法的理論和工程應(yīng)用，(英)布瑞比亞等著，北京－國防工業(yè)出版社，1988。邊界元理論及應(yīng)用，楊德全,趙忠生，北京-理工大學(xué)出版社，2002。1.2MathematicalPreliminaries(數(shù)學(xué)預(yù)備知識)Summationnotation:

Repeatedsubscriptsimplysummationofalltermsintherange,e.g.,

for2Dproblems,andfor3Dproblems.Therelationshipbetweenthesurfacetractionsandstressescanbeexpressedaswhichcanbewritteninthecomponentformas(i

=

1):(i

=

2):(i

=

3):Kroneckerdeltasymbolwheni=jwhenijInmatrixnotation,itsequivalentistheidentitymatrix[I].

Diracdeltafunction（impulsefunction）

wherepisthesingularpointand

denotesavanishinglysmallradiusofintegrationaroundthissingularity.

Gauss’theoremwhere

istheproblemdomain,

istheboundaryandisthei-thcomponentoftheunitoutwardnormalvector.

ThereductionofcertaindomainintegralstosurfaceintegralsiscentraltoBEM.Themethodof‘integrationbyparts’isanotherbasicmathematicaltechniquewhichisusedofteninconjunctionwithGauss’theorem.

Differentiationoftheproductoftwofunctions(f&g),withrespecttoxi,yields:‘integrationbyparts’statement:UsingGauss’theoremweobtain:

BEMfor1DProblems一維問題的邊界單元法One-dimensionalsecond-orderdifferentialequation:Multiplyequationbygandintegratebyparts

Integrationbypartstothelastterm

Letg(x,p)bethefundamentalsolutionoftheequation:

wherecisaconstantand.UsingDiracdeltafunctionpropertySubstitutingitintopreviousequationandusingwhereItfollowsthatxabConsideringuboundaryconditions:SolvingforandyieldsitfollowsthatFinally,weobtainorwhere作業(yè)：推導(dǎo)左邊u已知，右邊q已知情況下的關(guān)系式.2.BEMforPotentialProblems位勢問題的邊界單元法

PotentialflowproblemsHeatconductionproblemsSeepageproblemsElectromagneticfieldproblemsAcousticproblems(Helmholtzequations)GoverningEquation(Laplaceequation):inBoundaryconditions:onon(EssentialCondition)(NaturalCondition)wherewhere2.1BasicIntegralEquationsMultiplygoverningequationbyweightedfunction

:(weightedResidualFormulation)orwrittenas:whereInasimilarmannerPuttingthemtogetheryieldsLet

bethefundamentalsolutionoftheequation:

SubstitutingbackandfromI=0,weobtainUsingnotationswecanexpresstheintegralequationas

Fundamentalsolutionscanbederivedasfor2Dproblemsfor3DproblemsDerivationoftheFundamentalsolutionfor2DConsiderradialdistributionofand,wehaveTheLaplaceoperatorinthecylindricalcoordinatesystemcanbeexpressedasIntegratingthisequationtwiceyieldsTodetermine,integratingequationoveracircularareaaroundpointpwithradius,itfollowsthat(assuming)Finally,thefundamentalsolutionfor2DisderivedasSubstitutingintothisequationyieldsDerivationoftheFundamentalsolutionfor3DTheLaplaceoperatorinthesphericalcoordinatesystemcanbeexpressedasIntegratingthisequationtwiceyieldsConsiderradialdistributionofand,wehaveTodetermine,integratingequationoveraspherewithradius,itfollowsthatFinally,thefundamentalsolutionfor3DisderivedasSubstitutingintothisequationyieldsBoundaryIntegralEquationThederivedintegralequationisonlyvalidforinternalpoints.Tosetuptheintegralequationforboundarypoints,alimitprocessisperformed.Asimplewaytodothisistoconsiderthatthepointiisontheboundarybutthedomainitselfisaugmentedbyahemisphereofradius(in3D)asshowninthefollowingfigure.For3D(),theintegralaroundgives:NoticingthatitfollowsthatTheyproducewhatiscalledafreeterm.For2D(),theintegralaroundgives:NoticingthatitfollowsthatTheyproducewhatiscalledafreeterm.TheboundaryintegralequationcanbewrittenforsmoothboundarypointsasForboundarypointslocatedatacorner,sinceTheboundaryintegralequationcanbewrittenforsmoothboundarypointsastheboundaryintegralequationiswhere.Forthesakeofeasyunderstanding,theboundaryintegralequationcanbewrittenaswherePisthesourcepointandQthefieldpoint.BoundaryElementMethod(BEM)Tonumericallysolvetheboundaryintegralequation,theboundaryisdividedintoNsegmentsorelements.ConstantElementsFortheconstantelements,thevaluesofuandqareassumedtobeconstantovereachelementandequaltothevalueatthemid-elementnode.CollocatingthesourcepointPatthei-thnode,thediscretizedboundaryintegralequationbecomes:UsingthefollowingnotationsTheBEMequationcanbeexpressedasinwhich,andarevaluesofuandqatnodei.Letusnowcallwheni≠jwheni=jhencetheBEMequationcanbewrittenasIfthepositionofivariesfrom1toN,oneobtainsTheseequationscanbewritteninmatrixformaswhere[H]and[G]aretwoNXNmatricesand{u}and{q}arevectorsoflengthN,i.e.,Ifqisspecifiedatallnodes,itiseasytowritethesystemofequationasfollows:whereSimilarly,ifuisspecifiedatallnodes,wecanobtain:whereIfsomenodesarespecifiedwithuandotherwithq,wehavetorearrangethesystembymovingcolumnswithspecifiedvaluestotheright-handsideandcolumnswithunknownstotheleft-handsidetoformthefollowingsystemofequations.where{x}isavectorofunknownsu’sandq’svaluesandisfoundbymultiplyingthecorrespondingcolumnsbytheknownvaluesofu’sandq’s.Forexample,ifandaregivenandotheruareunknowns,wecanoperatethealgebraicequationsasMultiplyingtheright-handsidetogether,thefollowingmatrixequationcanbeobtained:whereEvaluationofInternalPotentialwherefor2Dproblemsfor3DproblemsiQrEvaluationofIntegralsIn2DProblems1.Evaluationofinfluencecoefficient2.EvaluationofinfluencecoefficientChangecoordinatestoalocaloneWhereistheelementlengthand.Thelastintegralisequalto1.So3.Evaluationofinfluencecoefficientwheredisthedistancefromtotheelement.12a)whenItfollowsthat12Forsegment:cSoForsegment:Itfollowsthatb)whenc)whend=0,d124.EvaluationofinfluencecoefficientOnedimensionalGaussquadratureformulasGaussquadraturewhere（）where

aretheGaussordinates;

aretheweights.nistheGaussorders;

Order(n)Ordinates()Weights(wk)1022

0.5773502692130

0.77459666920.88888888890.55555555564

0.3399810436

0.86113631160.65214515490.347854845150

0.5384693101

0.90617984590.56888888890.47862867050.23692688516

0.2386191861

0.6612093865

0.93246951420.46791393460.36076157300.1713244924(Stround&Secrest,1966)where

:theGaussordinates:theweightsn:theGaussorders

GaussintegrationruleforthelogarithmicallysingularfunctionsnOrdinates()Weights()10.25120.11200880610.60227690810.71853931900.281460680930.06389079310.36899706370.76688030390.51340455220.39198004120.094615406640.04144848010.24527491430.55616545350.84898239450.38346406810.38687531770.19043512690.039225487150.0291344722041170252050.67731417450.89477136100.29789347170.34977622650.23448829000.09893045950.0189115521ComputerCode(POCONBE)forPotentialProblemsusingConstantElementsMainprogramlist:RoutineOUTPTPCThisroutineoutputstheresults.Itfirstliststhecoordinatesoftheboundarynodesandthecorrespondingvaluesofpotentialanditsderivatives(orfluxes).Italsoprintsthevaluesofpotentialandfluxesatinternalpointsifanyhavebeenrequested.Example2.1:HeatFlowExampleDefinitionoftheproblemBoundaryconditionsInputdata(HEAT.INP):HEATFLOWEXAMPLE(12CONSTANTELEMENTS)125(Numbersofboundaryelementsandinternalnodes)0.0.2.0.4.0.6.0.6.2.6.4.(X(I),Y(I),I=1,N)6.6.4.6.2.6.0.6.0.4.0.2.10.(BoundaryconditionsKODE(I)andFI(I))10.10.00.00.00.10.10.10.0300.0300.0300.2.2.2.4.3.3.4.2.4.4.(Coordinatesofinternalpoints)Outputdata(HEAT.OUT):

BOUNDARYNODESXYPOTENTIALPOTENTIALDERIVATIVE0.10000E+010.00000E+000.25225E+030.00000E+000.30000E+010.00000E+000.15002E+030.00000E+000.50000E+010.00000E+000.47750E+020.00000E+000.60000E+010.10000E+010.00000E+00-0.52962E+020.60000E+010.30000E+010.00000E+00-0.48771E+020.60000E+010.50000E+010.00000E+00-0.52962E+020.50000E+010.60000E+010.47750E+020.00000E+000.30000E+010.60000E+010.15002E+030.00000E+000.10000E+010.60000E+010.25225E+030.00000E+000.00000E+000.50000E+010.30000E+030.52969E+020.00000E+000.30000E+010.30000E+030.48737E+020.00000E+000.10000E+010.30000E+030.52969E+02INTERNALPOINTSXYPOTENTIALFLUXXFLUXY0.20000E+010.20000E+010.20028E+03-0.50303E+02-0.14976E+000.20000E+010.40000E+010.20028E+03-0.50303E+020.14975E+000.30000E+010.30000E+010.15001E+03-0.50215E+02-0.25367E-050.40000E+010.20000E+010.99740E+02-0.50306E+020.14564E+000.40000E+010.40000E+010.99740E+02-0.50306E+02-0.14564E+00LinearElementsForthelinearelements,thevaluesofuandqareassumedtobelinearlyvaryingovereachelementandequaltothevaluesattheendsoftheelement.AfterdiscretizingtheboundaryintoaseriesofNelements,theboundaryintegralequationcanbewrittenasinwhich,theindexirepresentsthesourcepointcollocatingatnodei.whereistheinternalangleofthecornerinradians.Toevaluateboundaryintegrals,uandqareinterpolatedbytwonodalvaluesattheendsoftheelement.Takingthedimensionlesscoordinateasthevariablefrom-1to+1,thelinearvariationofucanbeexpressedasLettakevaluesatthetwoendsoftheelement.ItfollowsthatSolvingthisequationsetforkandbgiveswhereandarethenodalvaluesofu.UsingthefollowingnotationsSubstitutingthembacktothelineequationyieldsTheinterpolationformulationcanbewrittenasandarecalledshapefunctions.Similarly,whereNowtheintegralsinboundaryintegralequationscanbeexpressedaswhereAlgebraicSystemofEquationsFromWecanobtainforthei-thecollocationpointthat:whereorAswaspreviouslyshown,wecansimplywritewherewhenj

iwhenj=iandthewholesetinmatrixformbecomes[H]{u}=[G]{q}where[G]isnowanNX2Nrectangularmatrix.DeterminationoftheDiagonalTermsin

[H]forCloseDomainsThediagonaltermsin[H]

arecomputedimplicitly.Assumingaconstantpotentialoverthewholeboundary,thefluxmustbezeroandhence[H]{I}={0}Where{I}isavectorthatforallnodeshasaunitpotential.Thus,wehave(for)Whichgivesthediagonalcoefficientsintermsoftherestofthetermsofthe[H]matrix.DeterminationoftheDiagonalTermsin

[H]forInfiniteDomainsIfaunitpotentialisprescribedforaboundlessdomain,itfollowsthatthediagonaltermsare(for)Whichgivesthediagonalcoefficientsintermsoftherestofthetermsofthe[H]matrix.Sinceforacircularboundary,DeterminationoftheDiagonalTermsin

[G]Fortheelementwhichincludesthesingularity,theinfluentcoefficientscanbeexpressedasUsingtheintegrationvariabletransformationitcanbederivedthatThenthediagonaltermsbecomeIntegratingthemanalyticallyyieldsApplyboundaryconditionsAccordingtothesmoothnessoftheboundarynodes(includingcorners),fourdifferentcasesarepossibledependingontheboundaryconditions:Knownvalues:fluxes‘before’and‘a(chǎn)fter’thecorner.Unknownvalue:potential(b)Knownvalues:potential,andflux‘before’thecorner.Unknownvalue:flux‘a(chǎn)fter’thecorner(c)Knownvalues:potential,andflux‘a(chǎn)fter’thecorner.Unknownvalue:flux‘before’thecorner(d)Knownvalues:potential.Unknownvalues:flux‘before’and‘a(chǎn)fter’thecorner.Thereisonlyoneunknownpernodeforthefirstthreecases,andtwounknownsforcase(d).Aslongasthereisonlyoneunknownpernode,systemcanbereorderedinsuchawaythatalltheunknownsaretakentothelefthandsideandobtaintheusualsystemofNXNequations,i.e.[A]{X}={F}Where{X}isthevectorofunknownboundarypotentialsandfluxes,and{F}istheknownvectorcomputedbytheproductoftheknowboundaryconditionsandthecorrespondingcoefficientsofthe[G]or[H]matrices.Whenthenumberofunknownsatacornernodeistwo(case(d)),oneextraequationisneededforthenode.Theproblemcanbesolvedusingauxiliaryequationsordiscontinuouselements.DiscontinuousElementsToavoidtheproblemofhavingtwounknownfluxesatacomernode(forwhichonlyoneboundaryelementequationcanbewritten),thenodesofthetwolinearelementswhichmeetatthecornercanbeshiftedinsidethetwoelements.Thenodesremainastwodistinctnodesandoneequationcanbewrittenforeachnode.Discontinuouselementsarealsousefulforsituationsinwhichoneofthevariablestakesaninfinitevalueattheendoftheelement(forinstanceatareentrycornerorinfracturemechanicsapplications).Thevaluesofuandqatanypointonalinearelementhavebeendefinedintermsoftheirvaluesattheextremepointsbyequation:DiscontinuousElementsIfthetwonodesofanelementhavebeenshiftedfromtheendsdistancesaandbrespectively,theaboveequationcanbeparticularizedforthenodes.whereandarethelocalcoordinatesofthenodalpoints.InvertingtheaboveequationgiveswhereThenewinterpolationformulationnowbecomesThesaverelationcanbewrittenfortheflux:Whensolvingapotentialproblem,continuousanddiscontinuouselementscanbeusedtogetherinthesamemesh.Thetotalnumberofnodeswillbeequaltothetotalnumberofelementsplusoneadditionalnodepereachdiscontinuouselement.Thecoefficientciisequalto0.5forthenodesondiscontinuouselements.NearlySingularIntegrals高斯數(shù)值積分公式的誤差:

致密封頂層高溫陶瓷層過渡層氧化層粘結(jié)層基體氧化層過渡層粘結(jié)層9qm—高斯積分點數(shù)e—積分精度f—被積函數(shù)p為被積函數(shù)的奇異性階數(shù)，由表征，表示邊界單元在第個積分方向上的長度,R表示源點到邊界單元的最小距離,為積分精度。在大量數(shù)值調(diào)查結(jié)果的基礎(chǔ)上，得出了下列確定高斯點數(shù)的實用公式(Gao&Davies,2002)：式中重新整理后可得：3.單元子分法線單元1-2的子單元劃分示意圖幾乎奇異積分計算方法1.解析積分法2.積分變換法算例分析具有兩層涂層材料的圓柱形剛性基體外涂層楊氏模量與內(nèi)涂層楊氏模量的比值為1/2，內(nèi)外涂層的泊松比相同，采用平面應(yīng)力條件計算。解析解:A點和B點的徑向應(yīng)力隨外涂層厚度的變化

最大邊界元所需劃分的子單元數(shù)目QuadraticandHigherOrderElementsItisusuallymoreconvenientforarbitrarygeometriestoimplementsometypeofcurvilinearelements.Thesimplestofthesearethethreenodedquadraticelements.AfterdiscretizingtheboundaryintoaseriesofNelements,theboundaryintegralequationcanbewrittenasinwhich,theindexirepresentsthesourcepointcollocatingatnodei.Thepotentiala