边界单元法理论及程序设计1-位势边界元

BoundaryElementMethods—TheoryandProgramming

Xiao-WeiGAO,KaiYang(高效伟,杨恺)大连理工大学航空航天学院边界单元法理论及程序设计SchoolofAeronauticsandAstronauticsDalianUniversityofTechnology1.Introduction

1.1Whytolearnboundaryelementmethods(BEM)？为什么要学边界单元法(BEM)？Tosolveproblemsinwhichfiniteelementmethod(FEM)isinadequateorinefficient.GeometryisnotregularComputationalregionisInfiniteorsemi-infiniteFracturetipanalysis引言

BEMiseasytogeneratethedatarequiredtorunaproblemandcarryoutthemodificationsneededtoachieveanoptimumdesign.Whytolearnboundaryelementmethods(BEM)？为什么要学边界单元法(BEM)？OnlysurfacesoftheproblemneedtobediscretizedintoboundaryelementsMoreconvenienttodooptimumdesigncomputationSuitabeforparameterbackanalysisBEMismoreaccuratethanothernumericalmethods.Whytolearnboundaryelementmethods(BEM)？为什么要学边界单元法(BEM)？Semi-analysischaracteristicUseoffundamentalsolutionsGradientsofphysicalquantitieshavethesameaccuracylevelasthephysicalquantitiesselfBEMisveryefficienttosolveinfiniteorsemi-infiniteproblems.Whytolearnboundaryelementmethods(BEM)？为什么要学边界单元法(BEM)？InfiniteboundaryconditionscanbeautomaticallymodeledUseofInfiniteelementsUndergroundopeningPilefoundationproblemsAerodynamicproblemsBEMisrobusttosolvestress(orflux)concentrationproblems.Whytolearnboundaryelementmethods(BEM)？为什么要学边界单元法(BEM)？FundamentalsolutionsaresingularfunctionsFractureproblemsRigidfoundationproblemsBoundarylayerproblemsinfluidmechanicsReferences（参考文献）BoundaryElementsAnIntroductoryCourse，

BrebbiaCA,DominguezJ，ComputationalMechanicsPublications，1992。BoundaryElementProgramminginMechanics，

GaoXW&DaviesTG，CambridgeUniversityPress，2002。边界单元法的理论和工程应用，(英)布瑞比亚等著，北京－国防工业出版社，1988。边界元理论及应用，杨德全,赵忠生，北京-理工大学出版社，2002。1.2MathematicalPreliminaries(数学预备知识)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作业：推导左边u已知，右边q已知情况下的关系式.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’svaluesand{b}isfoundbymultiplyingthecorrespondingcolumnsbytheknownvaluesofu’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‘after’thecorner.Unknownvalue:potential(b)Knownvalues:potential,andflux‘before’thecorner.Unknownvalue:flux‘after’thecorner(c)Knownvalues:potential,andflux‘after’thecorner.Unknownvalue:flux‘before’thecorner(d)Knownvalues:potential.Unknownvalues:flux‘before’and‘after’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高斯数值积分公式的误差:

致密封顶层高温陶瓷层过渡层氧化层粘结层基体氧化层过渡层粘结层9qm—高斯积分点数e—积分精度f—被积函数p为被积函数的奇异性阶数，由表征，表示边界单元在第个积分方向上的长度,R表示源点到边界单元的最小距离,为积分精度。在大量数值调查结果的基础上，得出了下列确定高斯点数的实用公式(Gao&Davies,2002)：式中重新整理后可得：3.单元子分法线单元1-2的子单元划分示意图几乎奇异积分计算方法1.解析积分法2.积分变换法算例分析具有两层涂层材料的圆柱形刚性基体外涂层杨氏模量与内涂层杨氏模量的比值为1/2，内外涂层的泊松比相同，采用平面应力条件计算。解析解:A点和B点的径向应力随外涂层厚度的变化

最大边界元所需划分的子单元数目QuadraticandHigherOrderElementsItisusuallymoreconvenientforarbitrarygeometriestoimplementsometypeofcurvilinearelements.Thesimplestofthesearethethreenodedquadraticelements.AfterdiscretizingtheboundaryintoaseriesofNelements,theboundaryintegralequationcanbewrittenasinwhich,theindexirepresentsthesourcepointcollocatingatnodei.Thepotentiala