Introduction to the Finite Element Method in Electromagnetics

《Introduction to the Finite Element Method in Electromagnetics》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。

P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15IntroductiontotheFiniteElementMethodinElectromagneticsi P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15Copyright©2006byMorgan&ClaypoolAllrightsreserved.Nopartofthispublicationmaybereproduced,storedinaretrievalsystem,ortransmittedinanyformorbyanymeans—electronic,mechanical,photocopy,recording,oranyotherexceptforbriefquotationsinprintedreviews,withoutthepriorpermissionofthepublisher.IntroductiontotheFiniteElementMethodinElectromagneticsAnastasisC.Polycarpouwww.morganclaypool.com1598290460paperPolycarpou1598290479ebookPolycarpouDOI10.2200/S00019ED1V01Y200604CEM004APublicationintheMorgan&ClaypoolPublishers’seriesSYNTHESISLECTURESONCOMPUTATIONALELECTROMAGNETICSLecture#4FirstEdition10987654321PrintedintheUnitedStatesofAmericaii P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15IntroductiontotheFiniteElementMethodinElectromagneticsAnastasisC.PolycarpouIntercollege,CyprusSYNTHESISLECTURESONCOMPUTATIONALELECTROMAGNETICS#4MMorgan&ClaypoolPublishers&Ciii P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15Tomyparents,andtomywifeanddaughteriv P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15vABSTRACTThisserieslectureisanintroductiontothefiniteelementmethodwithapplicationsinelectro-magnetics.Thefiniteelementmethodisanumericalmethodthatisusedtosolveboundary-valueproblemscharacterizedbyapartialdifferentialequationandasetofboundaryconditions.Thegeometricaldomainofaboundary-valueproblemisdiscretizedusingsub-domainelements,calledthefiniteelements,andthedifferentialequationisappliedtoasingleelementafteritisbroughttoa“weak”integro-differentialform.Asetofshapefunctionsisusedtorepresenttheprimaryunknownvariableintheelementdomain.Asetoflinearequationsisobtainedforeachelementinthediscretizeddomain.Aglobalmatrixsystemisformedaftertheassemblyofallelements.Thislectureisdividedintotwochapters.Chapter1describesone-dimensionalboundary-valueproblemswithapplicationstoelectrostaticproblemsdescribedbythePoisson’sequation.Theaccuracyofthefiniteelementmethodisevaluatedforlinearandhigherorderelementsbycomputingthenumericalerrorbasedontwodifferentdefinitions.Chapter2describestwo-dimensionalboundary-valueproblemsintheareasofelectrostaticsandelectrodynamics(time-harmonicproblems).Forthesecondcategory,anabsorbingboundaryconditionwasimposedattheexteriorboundarytosimulateundisturbedwavepropagationtowardinfinity.Computationsofthenumericalerrorwereperformedinordertoevaluatetheaccuracyandeffectivenessofthemethodinsolvingelectromagneticproblems.Bothchaptersareaccompa-niedbyanumberofMatlabcodeswhichcanbeusedbythereadertosolveone-andtwo-dimensionalboundary-valueproblems.Thesecodescanbedownloadedfromthepublisher’sURL:www.morganclaypool.com/page/polycarpouThislectureiswrittenprimarilyforthenonexpertengineerortheundergraduateorgrad-uatestudentwhowantstolearn,forthefirsttime,thefiniteelementmethodwithapplicationstoelectromagnetics.Itisalsotargetedforresearchengineerswhohaveknowledgeofothernumericaltechniquesandwanttofamiliarizethemselveswiththefiniteelementmethod.Thelecturebeginswiththebasicsofthemethod,includingformulatingaboundary-valueproblemusingaweighted-residualmethodandtheGalerkinapproach,andcontinueswithimposingallthreetypesofboundaryconditionsincludingabsorbingboundaryconditions.Anotherimpor-tanttopicofemphasisisthedevelopmentofshapefunctionsincludingthoseofhigherorder.Insimplewords,thisserieslectureprovidesthereaderwithallinformationnecessaryforsomeonetoapplysuccessfullythefiniteelementmethodtoone-andtwo-dimensionalboundary-valueproblemsinelectromagnetics.Itissuitablefornewcomersinthefieldoffiniteelementsinelectromagnetics. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15KEYWORDSBoundary-valueproblems(BVPs),Finiteelementmethod(FEM),Galerkinapproach,Higherorderelements,Linearelements,Numericalmethods,Shape/interpolationfunctions,Weakformulationvi P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15viiContents1.One-DimensionalBoundary-ValueProblems...................................11.1Introduction.............................................................11.2ElectrostaticBVPandtheAnalyticalSolution...............................11.3TheFiniteElementMethod..............................................31.4DomainDiscretization....................................................51.5InterpolationFunctions...................................................51.6TheMethodofWeightedResidual:TheGalerkinApproach.................71.7AssemblyofElements...................................................131.8ImpositionofBoundaryConditions.......................................191.8.1DirichletBoundaryConditions...................................191.8.2MixedBoundaryConditions......................................221.9FiniteElementSolutionoftheElectrostaticBoundary-ValueProblem.......221.10One-DimensionalHigherOrderInterpolationFunctions...................291.10.1QuadraticElements..............................................301.10.2CubicElements.................................................331.11ElementMatrixandRight-Hand-SideVectorUsingQuadraticElements....361.12ElementMatrixandRight-Hand-SideVectorUsingCubicElements........391.13PostprocessingoftheSolution:QuadraticElements........................391.14PostprocessingoftheSolution:CubicElements............................471.15Software...............................................................482.Two-DimensionalBoundary-ValueProblems..................................512.1Introduction............................................................512.2DomainDiscretization...................................................522.3InterpolationFunctions..................................................542.3.1LinearTriangularElement........................................542.3.2BilinearQuadrilateralElement....................................592.4TheMethodofWeightedResidual:TheGalerkinApproach................612.5EvaluationofElementMatricesandVectors...............................662.5.1LinearTriangularElements.......................................672.5.2BilinearQuadrilateralElements...................................752.6AssemblyoftheGlobalMatrixSystem....................................86 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15viiiCONTENTS2.7ImpositionofBoundaryConditions.......................................902.8SolutionoftheGlobalMatrixSystem.....................................902.9PostprocessingoftheResults.............................................912.10ApplicationProblems....................................................922.10.1ElectrostaticBoundary-ValueProblem.............................922.10.2Two-DimensionalScatteringProblem.............................972.11HigherOrderElements................................................1052.11.1ANine-NodeQuadraticQuadrilateralElement...................1062.11.2ASix-NodeQuadraticTriangularElement........................1082.11.3ATen-NodeCubicTriangularElement...........................1102.12Software..............................................................111 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15ixPrefaceThisbookwaswrittenasanintroductorytexttothefiniteelementmethodinelectromagnetics.Thefiniteelementmethodhasbeenwidelyusedincomputationalelectromagneticsforthelast40–50yearswithanimpressivenumberofqualitypublicationsonthesubjectinthelate1980sand1990s.Itisahighlyversatilenumericalmethodthathasreceivedconsiderableattentionbyscientistsandresearchersaroundtheworldafterthelatesttechnologicaladvancementsandcomputerrevolutionofthetwentiethcentury.Themainconceptofthefiniteelementmethodisbasedonsubdividingthegeometricaldomainofaboundary-valueproblemintosmallersubdomains,calledthefiniteelements,andexpressingthegoverningdifferentialequationalongwiththeassociatedboundaryconditionsasasetoflinearequationsthatcanbesolvedcomputationallyusinglinearalgebratechniques.Thesubjectofthefiniteelementmethodinelectromagneticsisverybroadandcoversawiderangeoftopicsthatareimpossibletocoverinashortintroductorybook.Someofthesetopicsincludevectorelements,eigenvalueproblems,axisymmetricproblems,three-dimensionalscatteringandradiationproblems,microwaveandmillimeterwavecircuits,absorbingboundaryconditionsandtheperfectlymatchedlayer,hybridmethods,andafewmore.Thepurposeofthisbookisprimarilytheintroductionofthisnumericalmethodtotheundergraduatestudentandthenonexpertworkingengineerwhomaybeusingcommercialfiniteelementcodesorsimplyisinterestedinlearningthismethodforthefirsttime.Therefore,emphasiswasplacedonwritingabookthatislimitedinsizebutnotinsubstance,characterizedbysimplicityandclarity,freeofadvancedmathematicsandcomplexvariationalformulations,self-contained,andeffectiveinteachingthereaderthebasicsofthemethod.Itcanbeconsideredasafirstbookinlearningthefiniteelementmethodwithapplicationsinelectromagnetics.Moreadvancedbooksmayfollowtocoverspecifictopicsthatarenotdiscussedinthisintroductorybook.Thecontentofthebookisdividedintotwochapters.Chapter1presentsthefiniteelementformulationforone-dimensionalproblemsanditsspecificapplicationtoelectrostaticproblems.Initially,theformulationwascarriedoutusinglinearelements,whereastowardtheendofthechapter,theauthorintroduceshigherorderelementssuchasquadraticandcubic.ErroranalysisisalsopresentedinthischapterwherethenumericalerroriscomputedusingtwodifferentdefinitionsnamelythepercenterrorandtheerrorbasedontheL2norm.Itisimportanttoemphasizeherethatthroughouttheentirebook,alltheexpressionspresentedarederivedfromthebasics.Thereisnoexpressionthatispresentedwithoutderivation.Consequently,thereader P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-FMMOBK021-Polycarpou.clsApril29,200620:15xPREFACEcanfollowbetterthestepsinvolvedintheformulationofthemethodwithoutcreatinggapsanddoubts.Chapter2dealswiththefiniteelementformulationoftwo-dimensionalboundary-valueproblemsusingquadrilateralandtriangularelements.Thedevelopmentofhigherorderelementsisalsopresentedattheendofthechapter.ThefiniteelementformulationinvolvesimpositionofDirichlet,Neumann,ormixedboundaryconditions.Notethatafirst-orderabsorbingboundaryconditionthatisoftenusedtoterminatetheunboundeddomainofascatteringorradiationproblemisaspecialcaseofamixedboundarycondition.Theunderlinedformulationisappliedtoagenericsecond-orderpartialdifferentialequationwithasetofboundaryconditions:oneoftheDirichlettypeandoneofthemixedtype.Followingthismethodology,anytypeoftwo-dimensionalboundary-valueprobleminelectromagneticscanbetreatedusingtheunderlinedformulation.Thefiniteelementmethodintwodimensionswasappliedtoanelectrostaticproblemfirst,andthen,toascatteringproblemwhereafirst-orderabsorbingboundaryconditionwasusedtoterminatetheouterboundary.Anumberofplotsinthechapterillustratetheeffectivenessandtheaccuracyofthefiniteelementmethodascomparedtotheexactanalyticalsolution.ThenumericalerrorofthisformulationisquantifiedbyfollowingthesametypeoferroranalysisintroducedinChapter1.ThisbookonthefiniteelementmethodinelectromagneticsisaccompaniedbyanumberofcodeswrittenbytheauthorinMatlab.Thesearethefiniteelementcodesthatwereusedtogeneratemostofthegraphspresentedinthisbook.Specifically,therearethreeMatlabcodesfortheone-dimensionalcaseandtwoMatlabcodesforthetwo-dimensionalcasewhichcanbedownloadedfromthepublisher’sURL:www.morganclaypool.com/page/polycarpou.Thereadermayexecutethesecodes,modifycertainparameterssuchasmeshsizeorobjectdimensions,andvisualizetheresults.A.C.Polycarpou P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:141CHAPTER1One-DimensionalBoundary-ValueProblems1.1INTRODUCTIONInthischapter,thefiniteelementmethod(FEM)willbeappliedtoanelectrostaticboundary-valueproblem(BVP)inonedimension.Thereasonforchoosingtostartwithone-dimensional(1-D)problemsistohelpthereaderwalkthroughallthestepsoftheFEMwithouthavingtodealwithextensivemathematicalderivationsandgeometricalcomplexities.Thisway,thereaderwillgainabetterunderstandingoftheentirenumericalprocedureandgathersufficientknowledgetotackle2-Dand3-DBVPs.ThevalidityandaccuracyoftheFEMwillbeevaluated(aposteriorierroranalysis)bycomparingthenumericalresultwiththeexactanalyticalsolution.Therefore,beforeproceedingwithadetailedpresentationoftheFEManditsapplicationtothespecificBVP,itisinstructivethatwefirstputanefforttoobtainanalyticallythesolutiontotheproblemathand.Thiswillprovidethemeansforcomparisonandvalidationofthenumericalsolution.1.2ELECTROSTATICBVPANDTHEANALYTICALSOLUTIONProblemdefinition:Considertwoinfiniteinextentparallelconductingplatesthatarepositionednormaltothex-axisandseparatedbyadistanced,asshowninFigure1.1.OneplateismaintainedatafixedpotentialV=V0andthesecondplateismaintainedatV=0(ground).Theregionbetweentheplatesisfilledwithanonmagneticmediumhavingadielectricconstantεrandauniformelectronvolumechargedensityρv=−ρ0.Obtaintheanalyticalexpressionsfortheelectric(orelectrostatic)potentialandtheelectricfieldintheregionbetweenthetwoparallelplates.Analyticalsolution:ThepotentialdistributionatanypointbetweenthetwoplatesisgovernedbyPoisson’sequationρv∇(εr∇V)=−(1.1)ε0 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:142INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFIGURE1.1:GeometryoftheelectrostaticBVPsubjecttoasetofboundaryconditionsV(0)=V0(1.2)V(d)=0Forasimplemedium(homogeneous,linear,isotropic),Poisson’sequationinonedimensioncanbemoreconvenientlywrittenasd2Vρ0=(1.3)dx2εrε0whereρvwasreplacedby−ρ0.Integrating(1.3)twice,oneobtainsρ02V(x)=x+c1x+c0(1.4)2εrε0wherec1andc0areconstantstobedeterminedfromthesetofDirichletorotherwiseknownasessentialboundaryconditions.Thus,imposingthetwoboundaryconditionsin(1.2),theanalyticalsolutiontakestheformρ02ρ0dV0V(x)=x−+x+V0(1.5)2εrε02εrε0d P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS3TheelectricfieldexpressionisobtainedbytakingthenegativegradientoftheelectricpotentialdV(x)E(x)=−∇V=−aˆx(1.6)dxwhichresultsinE(x)=aˆV0ρ0dρ0xx+−(1.7)d2εrε0εrε0indicatingthattheelectricfieldisafunctionofx-coordinateandisdirectedalongthex-axis.ItisalsoimportanttonoticeherethattheelectricpotentialforthisparticularBVPisaquadraticfunctionofxwhereastheelectricfieldisalinearfunctionofx.1.3THEFINITEELEMENTMETHODTheFEM[1–5]isanumericaltechniquethatisusedtosolveBVPsgovernedbyadifferentialequationandasetofboundaryconditions.Themainideabehindthemethodistherepre-sentationofthedomainwithsmallersubdomainscalledthefiniteelements.Thedistributionoftheprimaryunknownquantityinsideanelementisinterpolatedbasedonthevaluesatthenodes,providednodalelementsareused,orthevaluesattheedges,incasevectorelementsareused.Theinterpolationorshapefunctionsmustbeacompletesetofpolynomials.Theaccuracyofthesolutiondepends,amongotherfactors,ontheorderofthesepolynomials,whichmaybelinear,quadratic,orhigherorder.Thenumericalsolutioncorrespondstothevaluesoftheprimaryunknownquantityatthenodesortheedgesofthediscretizeddomain.Thesolutionisobtainedaftersolvingasystemoflinearequations.Toformsuchalinearsystemofequations,thegoverningdifferentialequationandassociatedboundaryconditionsmustfirstbeconvertedtoanintegro-differentialformulationeitherbyminimizingafunctionalorusingaweighted-residualmethodsuchastheGalerkinapproach.Thisintegro-differentialformulationisappliedtoasingleelementandwiththeuseofproperweightandinterpolationfunctionstherespectiveelementequationsareobtained.TheassemblyofallelementsresultsinaglobalmatrixsystemthatrepresentstheentiredomainoftheBVP.Assaidinthepreviousparagraph,therearetwomethodsthatarewidelyusedtoobtainthefiniteelementequations:thevariationalmethodandtheweighted-residualmethod.ThevariationalapproachrequiresconstructionofafunctionalwhichrepresentstheenergyassociatedwiththeBVPathand.Afunctionalisafunctionexpressedinanintegralformandhasargumentsthatarefunctionsthemselves.Manyengineersandscientistsrefertoafunctionalasbeingafunctionoffunctions.AstableorstationarysolutiontoaBVPcanbeobtainedbyminimizingormaximizingthegoverningfunctional.Suchasolutioncorrespondstoeitheraminimumpoint,amaximumpoint,orasaddlepoint.Inthevicinityofsuchapoint,thenumericalsolutionis P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:144INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSstablemeaningthatitisratherinsensitivetosmallvariationsofdependentparameters.Thistranslatestoasmallernumericalerrorcomparedtoasolutionthatcorrespondstoanyotherpoint.Theprocessofminimizingormaximizingafunctionalinvolvestakingpartialderivativesofthefunctionalwithrespecttoeachofthedependentvariablesandsettingthemtozero.Thisformsasetofequationsthatcanbediscretizedwiththeproperchoiceofsubdomaininterpolationfunctionstogeneratethefiniteelementequations.Thesecondmethod,whichistheonefollowedinthisbook,isaweighted-residualmethodwidelyknownastheGalerkinmethod.ThismethodbeginsbyformingaresidualdirectlyfromthepartialdifferentialequationthatisassociatedwiththeBVPunderstudy.Simplystated,thismethoddoesnotrequiretheuseofafunctional.Theresidualisformedbytransferringalltermsofthepartialdifferentialequationononeside.Thisresidualisthenmultipliedbyaweightfunctionandintegratedoverthedomainofasingleelement.Thisisthereasonwhythemethodistermedasweighted-residualmethod.Ifthedifferentialequationisofsecondorder,asisthecasewithalltheproblemsconsideredinthisbook,itisrequiredthattheshapefunctionsusedtointerpolatetheprimaryunknownquantitybetwicedifferentiable.Thisrequirementisweakenedbyusingintegrationbypartsanddistributingthesecondderivativeequallybetweentheweightfunctionsandtheinterpolationfunctions.Inthisway,theassociatedweightfunctionsandinterpolationfunctionsarerequiredtobeonlyoncedifferentiable.Duetothisweakenedrequirement,theoutlinedformulationisalsoreferredtoastheweakformulation.Inaddition,iftheweightfunctionsarechosenfromthesamesetoffunctionsastheinterpolationfunctions,theunderlinedweighted-residualmethodiscalledGalerkinmethod.Inthisbook,wedecidedtofollowtheGalerkinapproachratherthanthevariationalapproach.ThereasonstemsfromthefactthattheGalerkinapproachissimpleandstartsdirectlyfromthegoverningdifferentialequation.Consequently,abeginnerwillhavelessdifficultyunderstandandcomprehendthestepsinvolvedintheformulationofthismethod.Onthecontrast,variationalmethodsrequireknowledgeofvariationalprinciples[6–8]inorderforsomeonetobeabletoconstructafunctional.Forsomewell-knownBVPs,thecorrespondingfunctionalisoftenavailable,buttherearecasesthatitisnecessarytoconstructoneusingvariationaltechniques.Toavoidthetediousprocedureofconstructingafunctionalandtheassociatedmathematicalcomplexities,itwasconsideredmoreappropriatetoimplementtheGalerkinapproachinsteadofthevariationalapproach.ThemajorstepsinvolvedintheapplicationoftheGalerkinFEMforthesolutionofaBVParethefollowing:•Discretizethedomainusingfiniteelements.•Chooseproperinterpolationfunctions(otherwiseknownasshapefunctionsorbasisfunctions). P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS5123Nex1324Nn–1NnFIGURE1.2:Discretizationofthe1-Ddomain•Obtainthecorrespondinglinearequationsforasingleelementbyfirstderivingtheweakformulationofthedifferentialequationsubjecttoasetofboundaryconditions.•Formtheglobalmatrixsystemofequationsthroughtheassemblyofallelements.•ImposeDirichletboundaryconditions.•Solvethelinearsystemofequationsusinglinearalgebratechniques.•Postprocesstheresults.Thesestepswillbefollowedone-by-oneinordertosolvetheelectrostaticBVPathand.1.4DOMAINDISCRETIZATIONThedomainoftheproblemathandcorrespondstoastraightlinealongthex-axisextendingfromx=0tox=d.AsshowninFigure1.2,thedomainissubdividedintoNelinesegmentscalledthefiniteelements.Theseelementsconstitutethefiniteelementmesh.Theelementnumberisshowncircledinthefigure.Eachelementhastwonodes;therefore,thetotalnumberofnodesinthedomainisNn=Ne+1.Dependingontheorderofshapefunctions,itmaybenecessarytointroduceadditionalnodesinsideeachelement.Suchelementsareknownashigherorderelements.Forlinearelements,thereareonlytwonodesthatarelocatedattheendpointsofthesegment.Thefiniteelementsdonothavetobeofthesamelength.Anelementisallowedtohaveanarbitrarylengthtoprovidetheabilitytogenerateadensermeshnearregionswherethesolutionisexpectedtohaverapidspatialvariations.Inaddition,thediscretizationofthedomainallowstheweakformulationoftheproblemtobeappliedtoeachelementseparately,thusallowingustodefinedistinctelementvaluesformaterialpropertiesandsources.Thisoffersgeneralityandversatilitytothemethod.1.5INTERPOLATIONFUNCTIONSTheapproximatesolutionobtainedfromtheweakformulationoftheproblemoverthedis-cretizeddomainisrepresentedinsideanelementbyasetofinterpolationfunctionsotherwiseknownasshapefunctions.Sincetheweakformulationoftheproblem,asitwillbeoutlinedinthefollowingsection,containsfirst-orderderivativesoftheprimaryunknownquantity,thechoseninterpolationfunctionsmustbecontinuouswithintheelementandatleastoncedifferentiable. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:146INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSee212xξxexe–1+112(a)(b)FIGURE1.3:(a)Elementalongx-axis.(b)Elementalongξ-axis(masterelement)Thesimplestchoiceisapolynomialofdegreeone.Theuseofapolynomialinsteadofanyothertypeoffunctionallowsfortheintegrations,whichareproductsoftheweakformulation,tobeevaluatedmoreeasily.Besidescontinuityanddifferentiability,anotherimportantrequirementisthatthesepolynomialsmustbecomplete.Inotherwords,theymustconsistofallthelowerorderterms.Thisisessentialinorderforthesolutiontobeaccuratelyrepresentedbytheinterpolationfunctionsinsideanelement.Letusnowconsiderafiniteelement(linesegment),asillustratedinFigure1.3(a).Thiselementhascoordinatesxeandxe,whichcorrespondtolocalnodes1and2,respectively.The12coordinatetransformation2(x−xe)ξ=1−1(1.8)xe−xe21canbeusedtotransformtheelementalongthex-axistothemasterelement,showninFigure1.3(b),whichresidesontheξ-axis.Theξ-coordinateisalsoknownasthenaturalcoordinate.Asillustratedinthefigure,themasterelementhasafixedpositionalongthenaturalcoordinateaxis.Theleftnodeoftheelement(node1)mapstoξ=−1whereastherightnode(node2)mapstoξ=+1.Itisthereforeeasierforustointegrateafunctiononthenaturalcoordinatesystemratherthanontheregularcoordinatesystem.Inotherwords,bymappinganelementontothenaturalcoordinateaxis,thelimitsofintegrationinvolvedintheweakformulationdonotchangeeverytimeanewelementisconsidered.Thiswouldbethecaseifanintegralweretobeevaluatedforelementsresidingontheregularcoordinateaxis.Duetothisobservation,itisinstructivethatinterpolationfunctionsbederivedbasedonthemasterelementratherthanthelocalelement.Atanypointinsidethemasterelement(−1≤ξ≤1),theprimaryunknownquantity(inourcasetheelectrostaticpotentialV)canbeexpressedasV(ξ)=VeN(ξ)+VeN(ξ)(1.9)1122whereN1(ξ),N2(ξ)aretheinterpolationfunctionsthatcorrespondtonodes1and2,respectively,andVe,Vearethevaluesoftheprimaryunknownquantityatthetwonodesoftheelement.12Notethatthenumberofinterpolationfunctionsusedperelementisequaltothenumberof P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS7N1N21xx–1+1–1+1(a)(b)VeV1eV2x–1+1(c)FIGURE1.4:(a)InterpolationfunctionN1.(b)InterpolationfunctionN2.(c)InterpolationofVnodesordegreesoffreedom(dof)thatbelongtotheelement.Inourcase,therearetwonodesperelementand,therefore,itisnecessarytousetwointerpolationfunctions.Fortheexpressionin(1.9)tobevalid,theinterpolationfunctionN1(ξ)mustbeequaltounityatnode1(i.e.,atξ=−1)andzeroatnode2(i.e.,atξ=+1),whereastheinterpolationfunctionN2(ξ)mustbeequaltounityatnode2(i.e.,atξ=+1)andzeroatnode1(i.e.,atξ=−1).Ifthisisnotthecase,then,theprimaryunknownquantityVwillnotbecontinuousacrosselementboundaries.Basedonthisclarification,thetwolinearinterpolationfunctionscanbeexpressedas1−ξN1(ξ)=(1.10)21+ξN2(ξ)=(1.11)2PlotsoftheseinterpolationfunctionsoverthemasterelementandalinearinterpolationoftheelectrostaticpotentialV,asdescribedby(1.9),areshowninFigure1.4.1.6THEMETHODOFWEIGHTEDRESIDUAL:THEGALERKINAPPROACHThefiniteelementequationsarederivedbyfirstconstructingtheweightedresidualofthegovern-ingdifferentialequationasappliedtoasingleelemente.Thisprocessisknownasthemethodofweightedresidualorweighted-residualmethod.GiventheelementshowninFigure1.5,withnodecoordinatesx=xeandx=xe,ourobjectiveistoobtainanapproximatesolutionsuch12thatthegoverningdifferentialequationissatisfiedinaweighted-integralsense.AsmentionedinSection1.5,asuitableformofsuchanapproximatesolutionspanningthedomainofthe P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:148INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSe12exex12FIGURE1.5:One-dimensionalfiniteelementelementeisapolynomialwhichsatisfiescertainrequirements.Thus,thesolutioninsideecanbewrittenasnV=veN(x)(1.12)jjj=1wherevearethesolutionvaluesatthenodesoftheelement,N(x)arethegoverninginterpo-jjlationfunctions,andnisthenumberofnodesinthedomainoftheelement;fora1-Dlinearelement,n=2.Theweightedresidualisformedbymovingalltermsofthedifferentialequationin(1.1)ononeside,multiplyingbyaweightfunctionw(x),andintegratingoverthedomainoftheelement.Theresultingweightedresidualforasingleelementhasthefollowingform:xe2ddVre=wεe+ρdx(1.13)vxedxdx1whereεeisthepermittivityofthemediuminsidetheelementgivenbyεe=εeε(1.14)r0WewillnowseekanapproximatesolutionforVbyforcingtheweightedresidualtobezero,thusobtainingaweighted-integralequationxe2ddVwεe+ρdx=0(1.15)vxedxdx1Foreachchoiceofweightfunctionw(x),anequationisformedwithunknownsbeingthevaluesoftheelectrostaticpotentialatthenodesoftheelement.Tobeabletosolvefortheunknowns,itisimportanttohaveasmanyindependentweightfunctionsasthenumberofnodesorunknowns.Itisoftenpreferredthattheweightfunctionsw(x)beidenticaltotheinterpolationorshapefunctionsN(x).ThisapproachisknownastheGalerkinfiniteelementmethodanditistheapproachthatwillbefollowedheretosolvethe1-DBVPbasedonPoisson’sequation.Fromtheweighted-integralequationin(1.15),itisevidentthattheinterpolationfunc-tionsmustbeatleastquadratic,andtherefore,twicedifferentiable.Toovercomethisstringent P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS9requirementandbeabletouselinearinterpolationfunctionsaswell,integrationbypartscanbeusedtotradethedoubledifferentiationonV(x)toanevenlydistributedsingledifferentiationonbothV(x)andw(x).IntegrationbypartsstatesthatbbUdV=UV|b−VdU(1.16)aaaUsingthisidentity,thefirsttermoftheweighted-integralequationcanbeconvenientlywrittenasxexex2exe2dedV2edVedV1dwedVwεdx=wdε=wε−εdxxedxdxxedxdxxexedxdx1111(1.17)Thus,theweakformulationofthegoverningdifferentialequationcanbeexpressedasxexex2e2dwdV2dVεedx−wρdx−wεe=0(1.18)vxedxdxxedxxe111Theaboveweakformulationrequiresthattheinterpolationfunctionsbeatleastoncedifferen-tiable.ThelasttermisevaluatedatthetwoendpointsoftheelementresultinginxedV2dVdVwεe=wxeεe−wxeεe(1.19)dx2dx1dxxex=xex=xe121Fromelectrostatics[9],itisknownthattheelectricfluxdensityorelectricdisplacementisdefinedalongthex-axisasdVDe=−εe(1.20)xdxSubstituting(1.20)into(1.19),thelatterbecomesxedV2wεe=−wxeDexe+wxeDexe(1.21)dx2x21x1xe1Asaresultofthisobservation,theweakformofthedifferentialequationbecomesxexe2dwdV2εedx−wρdx+wxeDexe−wxeDexe=0(1.22)v2x21x1xedxdxxe11Thefiniteelementequationsforasinglefiniteelemente,accordingtotheGalerkinapproach,canbeobtainedbysubstitutingtheapproximatesolution,givenby(1.12),intheweakformofthedifferentialequationandtakingtheweightfunctionsasbeingidenticaltotheinterpolationfunctions.Thus,for1-Dlinearelements,thefiniteelementequations,afterwesubstitutethe P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1410INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSapproximatesolutionforn=2andtheappropriateweightfunctions,takethefollowingform:xe2xe2dN1eedNj2eeeeeeεvjdx=N1ρvdx+N1x1Dxx1−N1x2Dxx2xedxdxxe1j=11xe2xe(1.23)2dN2eedNj2eeeeeeεvjdx=N2ρvdx+N2x1Dxx1−N2x2Dxx2xedxdxxe1j=11Basedonthepropertiesoftheinterpolationfunctionsintroducedintheprevioussection,itisdeducedthatNxe=111Nxe=021(1.24)Nxe=012Nxe=122Asaresult,thefiniteelementequationsin(1.23)canbeexpressedasxe2xe2dN1eedNj2eeεvjdx=N1ρvdx+Dxx1xedxdxxe1j=11xe2xe(1.25)2dN2eedNj2eeεvjdx=N2ρvdx−Dxx2xedxdxxe1j=11Thetwoalgebraicequationsin(1.25)canbeconvenientlywritteninmatrixformKeKevefeDe1112111=+(1.26)KeKevefe−De2122222wherexe2dNidNjKe=εedx(1.27)ijxedxdx1andxe2fe=Nρdx(1.28)iivxe1fori=1,2andj=1,2.Thesecondelementright-hand-sidevectorin(1.26)isdenotedbyed;i.e.,Dee1d=(1.29)−De2 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS11Thematrixsystemin(1.26)correspondstoasingleelementandnottotheentiredomainoftheBVP.Rememberthattheentiredomainisacollectionofelements.Theelementmatricesandvectorsrepresentingeachfiniteelementinthedomainmustbeassembledaccordingtotheconnectivityoftheelementsinordertoformtheglobalmatrixsystem.Theassemblyprocesswillbeexplainedinthefollowingsection.TheentriesoftheelementmatrixKeareobtainedbyevaluatingtheintegralin(1.27)alongthex-axiswithlimitsx=xeandx=xe.However,inSection1.5,weintroducedanew12coordinatesystem,namelythenaturalcoordinatesystem(ξ-coordinate),wheretheintegrationcanbeevaluatedmoreconveniently.Specifically,itcanbeshownfrom(1.8)that22dξ=dx=dx(1.30)xe−xele21whereleisthelengthofthefiniteelement(linesegment).Thus,ledx=dξ(1.31)2Inaddition,usingthechainruleofdifferentiation,thefirstderivativeoftheinterpolationfunctionNi(x)withrespecttoxcanbewrittenasdNidNidξ2dNi==(1.32)dxdξdxledξThus,theintegrationshownin(1.27)canbeequivalentlywrittenas+12dN2dNleeiejK=εdξijledξledξ2−1(1.33)2+1dNdNiej=εdξle−1dξdξConcerningtheelectrostaticproblemathand,themediumpermittivityεoftheregionbetweenthetwoparallelplatesisconstant,andtherefore,itcanbeplacedoutsidetheintegral.Thus,onecanwritethate+1e2εdNidNjKij=dξfori,j=1,2(1.34)le−1dξdξwhere1−ξN1(ξ)=2(1.35)1+ξN2(ξ)=2 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1412INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFromtheabovedefinitionofthetwolinearinterpolationfunctionswithrespecttoξ,itisevidentthatdN11=−(1.36)dξ2anddN21=+(1.37)dξ2Asaresult,thecorrespondingentriesoftheelementmatrixKearegivenbyεe+1−1Ke=(1.38)le−1+1Inasimilarway,theentriesoftheelementright-hand-sidevectorfecanbeevaluatedusing(1.28).FortheelectrostaticBVPathand,thevolumechargedensityρvwasassumedconstantandequalto−ρ0.Thus,usingthecoordinatetransformationfromtheregularcoordinatesystem(x-axis)tothenaturalcoordinatesystem(ξ-axis),theintegralin(1.28)canbemoreconvenientlywrittenase+1elρ0fi=−Ni(ξ)dξ(1.39)2−1Substituting(1.35)into(1.39),oneobtains+1leρ+11−ξleρξ2leρe000f1=−dξ=−ξ−=−(1.40)2−1242−12+1leρ+11+ξleρξ2leρe000f2=−dξ=−ξ+=−(1.41)2−1242−12Therefore,vectorfecanbeexpressedasleρ1e0f=−(1.42)21Thegoverningmatrixsystemforasingleelementisgivenbyεe+1−1veleρ1De1=−0+1(1.43)le−1+1ve21−De22 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS131.7ASSEMBLYOFELEMENTSIntheprevioussection,wepresentedtheweakformulationofthegoverningdifferentialequationbasedontheGalerkinapproachandthedevelopmentoftherespectivefiniteelementequationsforasingleelement.Solution,however,oftheBVPmandatestheassemblyofallfiniteelements,basedontheelementconnectivityinformation,andtheformationofaglobalmatrixsystemoflinearequations.ThisglobalmatrixsystemcanbesolvedtoobtainanumericalsolutiontotheproblemoncetheDirichletboundaryconditionsareimposed.TheimpositionofDirichletboundaryconditionsisthetopicofthefollowingsection.Toexplaintheprocessofassembly,considerFigure1.2,whichshowsthediscretizationofthedomainintoNefiniteelements.Eachelementehastwonodeswithlocalnodenumbers1and2.ThetotalnumberofnodesinthedomainisNn=Ne+1,whicharenumberedusingaglobalnodenumberschemestartingfrom1allthewaytoNn.Accordingtotheweakformulationpresentedintheprevioussectionandassuminglinearfiniteelements,asetoftwoequationsisproducedforeachelementinthedomain.Ingeneral,foranarbitraryelemente,theresultingtwoequationsaregivenbyKeve+Keve=fe+De(1.44)11112211Keve+Keve=fe−De(1.45)21122222Specifically,fore=(1):(1)(1)(1)(1)(1)(1)Kv+Kv=f+D(1.46)11112211(1)(1)(1)(1)(1)(1)Kv+Kv=f−D(1.47)21122222e=(2):(2)(2)(2)(2)(2)(2)Kv+Kv=f+D(1.48)11112211(2)(2)(2)(2)(2)(2)Kv+Kv=f−D(1.49)21122222e=(3):(3)(3)(3)(3)(3)(3)Kv+Kv=f+D(1.50)11112211(3)(3)(3)(3)(3)(3)Kv+Kv=f−D(1.51)21122222...e=(Ne):K(Ne)v(Ne)+K(Ne)v(Ne)=f(Ne)+D(Ne)(1.52)11112211K(Ne)v(Ne)+K(Ne)v(Ne)=f(Ne)−D(Ne)(1.53)21122222 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1414INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFromFigure1.2,itcanbeseenthatthesecondnodeofelementeisthesameasthefirstnodeofelemente+1.Inotherwords,ve=ve+1(1.54)21ormorespecifically(1)(2)v2=v1=V2(2)(3)v2=v1=V3(3)(4)(1.55)v2=v1=V4...whereVncorrespondstothenumericalsolutionattheglobalnodenumbern.Substitutingthisnumericalsolutioninto(1.46)–(1.53)resultsinthefollowingsystemofequations:e=(1):(1)(1)(1)(1)K11V1+K12V2=f1+D1(1.56)(1)(1)(1)(1)K21V1+K22V2=f2−D2(1.57)e=(2):(2)(2)(2)(2)K11V2+K12V3=f1+D1(1.58)(2)(2)(2)(2)K21V2+K22V3=f2−D2(1.59)e=(3):(3)(3)(3)(3)K11V3+K12V4=f1+D1(1.60)(3)(3)(3)(3)K21V3+K22V4=f2−D2(1.61)...e=(Ne):K(Ne)V+K(Ne)V=f(Ne)+D(Ne)(1.62)11Ne12Ne+111K(Ne)V+K(Ne)V=f(Ne)−D(Ne)(1.63)21Ne22Ne+122Adding(1.57)to(1.58)yields(1)(1)(2)(2)(1)(2)(1)(2)K21V1+K22+K11V2+K12V3=f2+f1+−D2+D1(1.64)Similarly,adding(1.59)to(1.60)yields(2)(2)(3)(3)(2)(3)(2)(3)K21V2+K22+K11V3+K12V4=f2+f1+−D2+D1(1.65) P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS15andsoon.Equations(1.56)and(1.63),i.e.,thefirstandthelastone,willnotbeaddedtoanyotherequation.Thus,thetotalnumberofequationsreducesfrom2NetoNe+1.Theresultingsystemofequationsisthefollowing:(1)(1)(1)(1)K11V1+K12V2=f1+D1(1)(1)(2)(2)(1)(2)(1)(2)K21V1+K22+K11V2+K12V3=f2+f1+−D2+D1(2)(2)(3)(3)(2)(3)(2)(3)K21V2+K22+K11V3+K12V4=f2+f1+−D2+D1...K(Ne−1)V+K(Ne−1)+K(Ne)V+K(Ne)V=f(Ne−1)+f(Ne)21Ne−12211Ne12Ne+121+−D(Ne−1)+D(Ne)21K(Ne)V+K(Ne)V=f(Ne)−D(Ne)21Ne22Ne+122(1.66)Thesealgebraicequationscanbewrittenmoreconvenientlyinmatrixform:⎡(1)(1)⎤⎧⎫K11K12⎪⎪V1⎪⎪⎢⎥⎪⎪⎪⎪⎢(1)(1)(2)(2)⎥⎪⎪⎪⎪⎢K21K22+K11K12⎥⎪⎪V2⎪⎪⎢⎥⎪⎪⎪⎪⎢K(2)K(2)+K(3)K(3)⎥⎪⎨V⎪⎬⎢21221112⎥3⎢⎥⎢⎢...............................................................⎥⎥⎪⎪..⎪⎪⎢⎥⎪⎪.⎪⎪⎢(Ne−1)(Ne−1)(Ne)(Ne)⎥⎪⎪⎪⎪⎪⎪⎪⎪⎢⎣K21K22+K11K12⎥⎦⎪⎪VNe⎪⎪⎪⎩⎪⎭K(Ne)K(Ne)VN+12122e(1.67)⎧⎫⎧⎫⎪⎪f(1)⎪⎪⎪⎪D(1)⎪⎪⎪⎪1⎪⎪⎪⎪1⎪⎪⎪⎪⎪⎪(1)(2)⎪⎪⎪⎪⎪⎪⎪⎪(1)(2)⎪⎪⎪⎪f+f−D+D⎪⎪21⎪⎪⎪⎪21⎪⎪⎪⎪(2)(3)⎪⎪⎪⎪(2)(3)⎪⎪⎨f+f⎬⎨−D+D⎬2121=+⎪⎪..⎪⎪⎪⎪..⎪⎪⎪⎪.⎪⎪⎪⎪.⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪f(Ne−1)+f(Ne)⎪⎪⎪⎪−D(Ne−1)+D(Ne)⎪⎪⎪⎪21⎪⎪⎪⎪21⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩f(Ne)⎭⎩−D(Ne)⎭22Equation(1.67)correspondstotheglobalmatrixsystemofequationsresultedfromtheassemblyprocessofallelementsinthefiniteelementdomain.Thesecondglobalright-hand-sidevectorin(1.67),denotedbyd,representstheassembledelectricfluxdensityatthenodes P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1416INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSofthefiniteelementmesh.Forexample,thesecondentryofthisvector,whichrepresentstheassembledelectricfluxdensityatglobalnodenumber2,isequalto(1)(2)(1)dV(2)dV−D+D=ε−ε=0(1.68)21dxx=x(1)dxx=x(2)21Theexpression(1.68)wouldhavebeenequaltozeroprovidedthatVcorrespondstotheexactsolutionoftheproblemwhichhasacontinuousderivativeacrosstheentiredomain.Onthecontrary,Vcorrespondstotheapproximatenumericalsolutionwhoseaccuracydependsonthechoiceofsubdomaininterpolationfunctions.ThefirstderivativeofVisnotnecessarilycontinuousacrosselements,regardlessofthefactthatmaterialparameters,suchasdielectricconstant,arethesameforallelements.Thecontinuityofthesecondaryvariableacrosselements,whichisdirectlyproportionaltothefirstderivativeoftheprimaryvariable,wasneverimposedinthefiniteelementformulationpresentedinSection1.6.Inparticular,thefirstderivativeofVisconstantoverthelengthofalinearelementand,asaresult,therewillbeastepdiscontinuityinthedistributionofthesecondaryvariable.Thisdiscontinuitywillappearatelementboundaries.However,asthenumberofelementsinthedomainisincreased,thisstepdiscontinuitytendstodecrease.Consequently,forasufficientlylargenumberofelementsinthedomain,itmakesabsolutesensetoclaimthat−D(1)+D(2)∼=0(1.69)21Asaresultofthisobservation,thesecondglobalright-hand-sidevectordcanbewrittenas⎧⎫(1)⎪⎪D1⎪⎪⎪⎪⎪⎪⎪⎪0⎪⎪⎪⎪⎪⎪⎨0⎬d=.⎪⎪..⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪0⎪⎪⎩(Ne)⎭−D2Foragivenfiniteelementmesh,theassemblyprocessmaybeautomatedbyforminganarray(ortable)holdingtheelementconnectivityinformation.Theelementconnectivityinformationrelatesthelocalnodenumbersofanelementtothecorrespondingglobalnodenumbers,asshowninTable1.1.Aglobalnodenumber,startingfrom1endingatNn(whereNnisthetotalnumberofnodesinthefiniteelementmesh),isuniquelyassignedtoeachnode.Fromtheweakformulationpresentedintheprevioussection,theelementcoefficientmatrixKewasfoundtobeεe+1−1Ke=(1.70)le−1+1 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS17TABLE1.1:ElementconnectivityinformationELEMENTLOCALNODE1LOCALNODE2112223334.........NeNn−1NnThiselementcoefficientmatrixmustbemapped,accordingtotheelementconnectivityinfor-mation,totheglobalcoefficientmatrixwithdimensionsNn×Nn.Heretheword“mapping”meansaddingtheentriesoftheelementcoefficientmatrixtothecorrespondingentriesoftheglobalcoefficientmatrix.NotethattheglobalcoefficientmatrixKmustbeinitializedtothezeromatrixbeforewebegintheassemblyprocess.Therefore,accordingtoTable1.1,entryKe11ofelement2mapstoentryKoftheglobalmatrix;entryKeofelement2mapstoentry2212K23oftheglobalmatrix,andsoon.Toillustratetheassemblyprocess,consideranexamplewherewehaveonlythreeelementsinthedomain,i.e.,Ne=3andNn=4.Asaresult,thedimensionsoftheglobalmatrixis4×4andcanbeformedbyaddingthecontributionofeachelementcoefficientmatrixtoaninitialized-to-zeroglobalKmatrix.Thus,⎡⎤⎡⎤⎡⎤+1−10000000000(1)⎢⎥(2)⎢⎥(3)⎢⎥ε⎢−1+100⎥ε⎢0+1−10⎥ε⎢0000⎥K=⎢⎥+⎢⎥+⎢⎥(1.71)l(1)⎣0000⎦l(2)⎣0−1+10⎦l(3)⎣00+1−1⎦0000000000−1+1whichisequalto⎡⎤ε(1)ε(1)−00⎢⎢l(1)l(1)⎥⎥⎢ε(1)ε(1)ε(2)ε(2)⎥⎢⎢−+−0⎥⎥⎢l(1)l(1)l(2)l(2)⎥K=⎢(2)(2)(3)(3)⎥(1.72)⎢εεεε⎥⎢0−+−⎥⎢l(2)l(2)l(3)l(3)⎥⎣(3)(3)⎦εε00−l(3)l(3) P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1418INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSTheelementright-hand-sidevectorisassembledtotheglobalright-hand-sidevectorinaverysimilarway.Eachentryfromtheelementright-hand-sidevectorisaddedtothecor-respondingentryoftheglobalright-hand-sidevectoraccordingtotheelementconnectivityinformationtabulatedinTable1.1.Notethattheglobalright-hand-sidevectormustbeinitial-izedtozerobeforethebeginningoftheassemblyprocess.Ignoringfornowthecontributionofelementvectordeandtakingintoaccountonlythecontributionbyelementvectorfe,whichisgivenby(1.42),theassembledglobalright-hand-sidevectorbisgivenby⎧⎫⎧⎫⎧⎫⎪⎪1⎪⎪⎪⎪0⎪⎪⎪⎪0⎪⎪⎪⎨⎪⎬⎪⎨⎪⎬⎪⎨⎪⎬l(1)ρ1l(2)ρ1l(3)ρ0000b=−−−(1.73)2⎪⎪0⎪⎪2⎪⎪1⎪⎪2⎪⎪1⎪⎪⎪⎩⎪⎭⎪⎩⎪⎭⎪⎩⎪⎭001whichcanalsobewrittenas⎧⎫⎪⎪l(1)⎪⎪⎪⎨⎪⎬ρl(1)+l(2)0b=−(1.74)2⎪⎪l(2)+l(3)⎪⎪⎪⎩⎪⎭l(3)Notethatρ0denotesthechargedistributionintheregionbetweentheplates,anditwasassumedconstantforallelements.Usingacomputer,theassemblyprocessoftheglobalcoefficientmatrixandtheglobalright-hand-sidevectorbeginsbyformingfirsta2-Darray,namedhereelmconn,withdi-mensionsNe×2,whichholdstheelementconnectivityinformationgiveninTable1.1.AsimplealgorithmwritteninMatlabillustratestheassemblyprocessoftheglobalmatrixandright-hand-sidevectorbyloopingthroughtheelementsofthediscretizeddomain:#InitializetheglobalmatrixK=sparse(Nn,Nn);b=zeros(Nn,1);#Loopthroughtheelementsfore=1:Ne#Doubleloopthroughthelocalnodesofeachelementfori=1:2forj=1:2K(elmconn(e,i),elmconn(e,j))=K(elmconn(e,i),elmconn(e,j))+Ke(i,j);endb(elmconn(e,i))=b(elmconn(e,i))+fe(i);endend P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS19Thenamesofthevariablesusedareself-explanatory,e.g.,elmconnstandsforthearrayholdingtheelementconnectivityinformation.Inaddition,forthefillingoftheglobalcoefficientmatrix,theMatlabcommandsparsewasusedinordertosaveoncomputermemorybyavoidingstoringthezeroentries.Notethatthemajorityofentriesintheglobalcoefficientmatrixwillbezeroafterthecompletionoftheassemblyprocess.Wecallsuchamatrixasparsematrix.Thesparsityoftheglobalcoefficientmatrixisattributedtothesubdomainnatureoftheshapefunctions.Inotherwords,agivenshapefunction,whichcorrespondstoaspecificnode,isnonzeroonlyinsidetheelementswherethisnodebelongsto.Fortheremainingelementsinthediscretizeddomain,thisspecificshapefunctioniszeroand,therefore,thereisnointeractionbetweenthespecificnodeandtheassociatednodesofthoseelements.Asaresult,thecorrespondingentriesoftheglobalcoefficientmatrixwillbezero.1.8IMPOSITIONOFBOUNDARYCONDITIONSThissectionoutlinestheprocedureusedtoimposeboundaryconditionsonthesetoflinearequationsobtainedfromtheweakformulationofthegoverningdifferentialequation.Priortoimposingboundaryconditions,theglobalmatrixsystemissingularand,thus,cannotbesolvedtoobtainauniquesolution.AnonsingularmatrixsystemisobtainedafterimposingtheboundaryconditionsassociatedwithagivenBVP.ThetwotypesofboundaryconditionsthatarediscussedinthissectionaretheDirichletoressentialboundaryconditionsandthemixedboundaryconditions.TheDirichletboundaryconditionsinvolveonlytheprimaryunknownvariablewhereasthemixedboundaryconditionsinvolveboththeprimaryunknownvariableanditsderivative.AnothertypeofboundaryconditionsistheNeumannboundaryconditionswhichcanbeconsideredasaspecialcaseofthemixedboundaryconditions.1.8.1DirichletBoundaryConditionsTheweakformulationofthegoverningdifferentialequationovertheentirefiniteelementdomainresultsinasystemofNlinearequationswithNunknowns.ForthenodalFEM,theNunknownscorrespondtotheNnodesofthedomain.Thus,thereisonedegreeoffreedom(dof)orunknownpernode.Foragenericfiniteelementmesh,notnecessarily1-D,theresultingmatrixsystemoflinearequationsis⎡⎤⎧⎫⎧⎫K11K12K13···K1N⎪⎪V1⎪⎪⎪⎪b1⎪⎪⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪⎢K21K22K23···K2N⎥⎪⎪⎨V2⎪⎪⎬⎪⎪⎨b2⎪⎪⎬⎢⎥⎢⎢K31K32K33···K3N⎥⎥V3=b3(1.75)⎢..........⎥⎪⎪⎪⎪..⎪⎪⎪⎪⎪⎪⎪⎪..⎪⎪⎪⎪⎣.....⎦⎪⎪.⎪⎪⎪⎪.⎪⎪⎩⎭⎩⎭KN1KN2KN3···KNNVNbN P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1420INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwheretheunknownquantityistheelectricpotentialVatthenodesofthefiniteelementmesh.Thematrixsystemin(1.75)canalsobewrittenasasetofNlinearequationswithNunknownsasshownbelow:K11V1+K12V2+K13V3+···+K1NVN=b1K21V1+K22V2+K23V3+···+K2NVN=b2K31V1+K32V2+K33V3+···+K3NVN=b3(1.76)...KN1V1+KN2V2+KN3V3+···+KNNVN=bNADirichletboundaryconditionisimposedataspecificnodeofthefiniteelementmesh.Forexample,theBVPathandhastwoDirichletboundaryconditionstobeimposed:V1=V0(1.77)andVN=0(1.78)Toimposeforexample(1.77),whichisassociatedwithnode1,thecorrespondinglinearequationin(1.76)mustbeeliminated;inthisparticularcase,thefirstequationmustbeeliminated.Then,wesubstitutethevalueV1=V0inalltheremainingN−1equations.Thisresultsinthefollowinglinearsystemofequations:K21V0+K22V2+K23V3+···+K2NVN=b2K31V0+K32V2+K33V3+···+K3NVN=b3.(1.79)..KN1V0+KN2V2+KN3V3+···+KNNVN=bNThefirstproducttermineachoftheaboveequationsisaconstanttermwhichcanbetransferredtotheright-handsideandobtainK22V2+K23V3+···+K2NVN=b2−K21V0K32V2+K33V3+···+K3NVN=b3−K31V0.(1.80)..KN2V2+KN3V3+···+KNNVN=bN−KN1V0Thetotalnumberoflinearequationshasbeenreducedbyone.Inotherwords,ifwehavetoimposeMDirichletboundaryconditions,thesizeofthefinallinearsystemofequationswillbereducedtoN−M.Thesetoflinearequationsin(1.80)canbeequivalentlyexpressedin P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS21matrixformas⎡⎤⎧⎫⎧⎫K22K23···K2N⎪⎪⎪⎪V2⎪⎪⎪⎪⎪⎪⎪⎪b2−K21V0⎪⎪⎪⎪⎢⎢KK···K⎥⎥⎨V⎬⎨b−KV⎬32333N33310⎢⎢...⎥⎥.=.(1.81)⎣....···..⎦⎪⎪⎪⎪..⎪⎪⎪⎪⎪⎪⎪⎪..⎪⎪⎪⎪⎩⎭⎩⎭KN2KN3···KNNVNbN−KN1V0Comparingtheglobalmatrixsystemin(1.81)withtheglobalmatrixsystemin(1.75),whichcorrespondstotheproblembeforeimposingthefirstDirichletboundarycondition,itisobviousthatthefirstrowofthematrixandright-hand-sidevector,aswellasthefirstcolumnofthematrix,wereeliminated.Theremainingentriesoftheright-hand-sidevectorwereupdatedaccordingtobi=bi−Ki1V0fori=2,3,...,N(1.82)Ingeneral,toimposeaDirichletboundaryconditionatagivennoden,letussayVn=V0(1.83)thenthrowofthesystemmatrix,thenthrowoftheright-hand-sidevector,andthenthrowoftheunknownvector,aswellasthenthcolumnofthesystemmatrixmustbeeliminated,whereastheremainingentriesoftheright-hand-sidevectormustbeupdatedaccordingtobi=bi−KinV0fori=1,2,...,N;i=n(1.84)Aftereliminatingagivenrow,alltherowsunderneathmustbeshiftedupwardbyoneposition.Similarly,aftereliminatingagivencolumn,allthecolumnsontherightmustbeshiftedtowardtheleftbyoneposition.Oncethematrixsystemissolved,onlyN−Munknownswillbedetermined;thevaluesoftheprimaryunknownquantityattheremainingMnodesareknownfromthesetofDirichletboundaryconditions.ThismethodofimposingDirichletboundaryconditionsisknownastheMethodofEliminationbecausethealgebraicequationthatcorrespondstotheglobalnodeatwhichtheDirichletboundaryconditionisimposedmustbeeliminated,thusreducingthesizeofthegoverningmatrixsystembyone.Fromaprogrammingpointofview,itismoreconvenienttonumbertheglobalnodesofthefiniteelementmeshinsuchawaythatthenodeswhichcorrespondtoDirichletboundaryconditionsappearlast.Thus,usingthemethodofeliminationinimposingDirichletboundaryconditions,itisnotnecessaryeverytimewedeletearowtohavetoshiftthebottomrowsupward.Thesameargumentappliestocolumnsaswell.Thiswillsavecomputationaltimeandmakeprogrammingsimplerandmorestraightforward. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1422INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS1.8.2MixedBoundaryConditionsAmixedboundarycondition,unlikeaDirichletorNeumannboundarycondition,involvesthevariableofthequantityunderstudyanditsderivative.AgenericformofamixedboundaryconditioncanbeexpressedasdVε+αV=β(1.85)dxwhereαandβareconstants.ForaNeumannboundarycondition,i.e.,dV/dx=0atagivenx-coordinate,constantsαandβaresettozero.AsitwillbeshowninChapter2,afirst-orderabsorbingboundarycondition(ABC),whichisusedinscatteringandradiationproblemstotruncatetheunboundedfreespace,isaformofamixedboundarycondition.Toimposethismixedboundaryconditionatanexteriornode,letussaytherightmostnodeofthedomainwithglobalnumberN,weneedtoremindourselvesthattheglobalright-hand-sidevectorbisasuperpositionoftwootherglobalvectorsnamelyfanddwhere⎧⎫⎧⎫⎪⎪(1)dV⎪⎪(1)⎪⎪−ε⎪⎪⎪⎪D1⎪⎪⎪⎪dxx=x(1)⎪⎪⎪⎪⎪⎪⎪⎪1⎪⎪⎪⎪0⎪⎪⎪⎪0⎪⎪⎨⎬⎨⎬d=..=..(1.86)..⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪0⎪⎪⎪⎪⎪⎪⎪⎪0⎪⎪⎪⎪⎩−D(Ne)⎭⎪⎪dV⎪⎪2⎪⎪⎩ε(Ne)⎪⎪⎭dxx=x(Ne)2Toimposethemixedboundaryconditiongivenby(1.85)attherightmostnodeofthefiniteelementmesh,wemustreplacethelastentryofvectordbydVε(Ne)=β−αV(1.87)Ndxx=x(Ne)2NotethatVNistheunknownpotentialattheNthnodeofthedomain.ThetermαVNisthereforetransferredtotheleft-handsidewhereasβstaysatthelastentrypositionofvectord.TransferringαVNtotheleft-handsideofthematrixsystemisequivalenttoaddingconstantαtotheKNNentryoftheglobalcoefficientmatrix.Adetailedexplanationonhowafirst-orderABCisimplementedina2-DscatteringproblemwillbeprovidedinChapter2.1.9FINITEELEMENTSOLUTIONOFTHEELECTROSTATICBOUNDARY-VALUEPROBLEMAftergoingthroughthemajorstepsoftheFEM,wearenowinthepositiontousethispowerfulnumericalmethodtosolvetheelectrostaticBVPathand.Rememberthattheobjectiveisto P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS2312341234345x(cm)x1=0x2=2x3=4x4=6x5=8FIGURE1.6:Discretizationofthedomainusingfourlinearfiniteelementscomputetheelectricpotentialdistributionbetweentwoparallelplatesseparatedbyadistancedandpositionednormaltothex-axis.TheleftmostplateismaintainedataconstantpotentialV0whereastherightmostplateisgrounded.Theregionbetweentheplatesischaracterizedbyadielectricconstantεrandauniformelectronchargedensity−ρ0.TheanalyticalsolutiontotheproblemwaspresentedinSection1.2.Theexactanalyticalexpressionsfortheelectricpotentialandelectricfieldthatexistintheregionbetweenthetwoplatesaregivenby(1.5)and(1.7),respectively.Forafiniteelementsimulationandcomparisonofthenumericalsolutionwiththeexactanalyticalsolution,itisrequiredthatcertainparametersbedefined.Forthesakeofcomputation,itisassumedthatεr=1V0=1V(1.88)d=8cmρ=10−8C/m30Now,considerthatthedomainbetweentheplatesisequallydividedintofourlinearfiniteelements,asshowninFigure1.6.Allelementsinthedomainarecharacterizedbythesamelengthleandthesamedielectricconstantεe.Thus,accordingtoSection1.6,theelementrcoefficientmatrixKeisgivenby8.85×10−12+1−1+1−1Ke==4.425×10−10(1.89)2×10−2−1+1−1+1whereεe=εeε=8.85×10−12F/mandle=2×10−2m.Allunitsareexpressedinthemet-r0ricsystem.Theelementright-hand-sidevectorfebecomes2×10−2×10−811fe=−=−10−10(1.90)211Thecontributionoftheright-hand-sidevectordetotheglobalright-hand-sidevector,aswasshownintheprevioussection,iszeroforallnodesexceptforthetwoendnodesofthedomain.However,atthesetwoendnodes,Dirichletboundaryconditionsmustbeimposedand, P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1424INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICStherefore,thecontributionbyvectordeiseffectivelydiscarded.Thus,basedontheassemblyprocesspresentedinSection1.7,theglobalmatrixsystemforthefiniteelementmesh,depictedinFigure1.6,becomes⎡⎤⎧⎫⎧⎫1−1000⎪⎪V1⎪⎪⎪⎪1⎪⎪⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪⎢−12−100⎥⎪⎨V2⎪⎬⎪⎨2⎪⎬−10⎢⎥−104.425×10⎢0−12−10⎥V3=−102(1.91)⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪⎣00−12−1⎦⎪⎪V4⎪⎪⎪⎪2⎪⎪⎪⎩⎪⎭⎪⎩⎪⎭000−11V51Dividingbothsidesby4.425×10−10,thematrixsystemcanbeequivalentlywrittenas⎡⎤⎧⎫⎧⎫1−1000⎪⎪V1⎪⎪⎪⎪−0.2259887⎪⎪⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪⎢−12−100⎥⎪⎨V2⎪⎬⎪⎨−0.4519774⎪⎬⎢⎥⎢0−12−10⎥V3=−0.4519774(1.92)⎢⎥⎪⎪⎪⎪⎪⎪⎪⎪⎣00−12−1⎦⎪⎪V4⎪⎪⎪⎪−0.4519774⎪⎪⎪⎩⎪⎭⎪⎩⎪⎭000−11V5−0.2259887ImposingtheDirichletboundaryconditionV=1atnode1eliminatestheentirefirstrow,includingthefirstrowoftheright-hand-sidevector,andthefirstcolumnofthecoefficientmatrix.Oncethisisdone,theright-hand-sidevectormustbeupdatedaccordingto(1.82),thusresultinginthefollowingreducedmatrixsystem:⎡⎤⎧⎫⎧⎫2−100⎪⎪V2⎪⎪⎪⎪0.5480226⎪⎪⎢⎥⎪⎨⎪⎬⎪⎨⎪⎬⎢−12−10⎥V3−0.4519774⎢⎥=(1.93)⎣0−12−1⎦⎪⎪V4⎪⎪⎪⎪−0.4519774⎪⎪⎪⎩⎪⎭⎪⎩⎪⎭00−11V5−0.2259887ThesecondboundaryconditionV=0atnode5isimposedbyeliminatingtheentirelastrowofthematrixsystemandthelastcolumnofthecoefficientmatrix.Updatingtheright-hand-sidevectorisneedlesssinceV5=0.Thus,thefinalglobalmatrixsystembecomes⎡⎤⎧⎫⎧⎫2−10⎪⎨V2⎪⎬⎪⎨0.5480226⎪⎬⎢⎥⎣−12−1⎦V3=−0.4519774(1.94)⎪⎩⎪⎭⎪⎩⎪⎭0−12V4−0.4519774Thisglobalmatrixsystemcanbesolvedusingseveraltechniquesmostofwhichcanbefoundincommonlinearalgebrabooks.OnesuchtechniqueisCramer’srule,whichisalso P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS25outlinedinvariousintroductorybooksonlinearalgebra[10].Thus,usingCramer’srule,onecansolvefortheelectricpotentialatthethreeinteriornodesofthefiniteelementdomain.Specifically,0.5480226−10−0.45197742−1−0.4519774−120.2881356V2===0.0720339(1.95)2−104−12−10−1220.54802260−1−0.4519774−10−0.45197742−1.6158192V3===−0.4039548(1.96)2−104−12−10−122−10.5480226−12−0.45197740−1−0.4519774−1.7118644V4===−0.4279661(1.97)2−104−12−10−12NotealsothatV1=1(1.98)andV5=0(1.99)Havingsolvedtheglobalmatrixsystem,theunknownelectricpotentialatthenodesofthefiniteelementmeshisfound.Toplottheelectricpotentialatintermediatepointsrequirestheuseoftheinterpolationorshapefunctionsemployedforeachfiniteelement.FortheBVPathand,linearinterpolationfunctionswereusedand,thus,thenumericalsolutionatintermediatepointsinsideanelementisgivenbyV(ξ)=VeN(ξ)+VeN(ξ)(1.100)1122 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1426INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwhere1−ξN1(ξ)=2(1.101)1+ξN2(ξ)=2and2x−xe1ξ=−1(1.102)xe−xe21Substituting(1.102)into(1.101)yieldsxe−x2N1(x)=eex−x21(1.103)x−xeN(x)=12eex−x21Consequently,theelectricpotentialatanypointinsideanelementcanbewrittenasxe−xx−xeV(x)=Ve2+Ve1(1.104)1xe−xe2xe−xe2121whereVeandVearethevaluesoftheelectricpotentialatthetwoendnodesoftheelement.12Theelectricfieldatanypointinsideanelementiscomputedbytakingthenegativegradientoftheelectricpotentialgivenby(1.104)E=−∇V(1.105)whichisequivalenttodV(x)E=−aˆx(1.106)dxsincetheelectricpotentialisonlyafunctionofthex-coordinate.Applying(1.106)on(1.104)yieldsVeVeE=−aˆ−1+2xxe−xexe−xe2121(1.107)Ve−VeVe−Ve=aˆ12=aˆ12xxe−xexle21wherele=xe−xeisthelengthoftheelement.Unliketheelectricpotential,whichiscontin-21uousacrosselementboundaries,theelectricfieldisdiscontinuous.Inaddition,duetotheuseoflinearinterpolationfunctions,theelectricfield(whichisproportionaltothegradientofthe P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS271.0FEMExact0.80.60.40.20.0Electricpotential(V)–0.2–0.40.00.010.020.030.040.050.060.070.08Distance(m)FIGURE1.7:ComparisonbetweentheelectricpotentialproducedbytheFEMusinglinearelementsandtheexactanalyticalsolutionelectricpotential)isconstantoveranelement.Useofhigherorderinterpolationfunctionswillresultinabetterrepresentationoftheelectricfieldinthefiniteelementdomain.Usingafour-elementmesh,theelectricpotentialovertheentiredomainiscomputedandplotted.ThefiniteelementsolutioniscomparedagainsttheanalyticalsolutionobtainedinSection1.2.ThecomparisonisshowninFigure1.7.Asillustrated,theelectricpotentialatthenodesofthefiniteelementmeshmatchesperfectlytheanalyticalsolution,whereasatintermediateevaluationpointsthereisadeviationbetweenthetwosolutions.Thereasonforthisdeviationstemsfromthefactthatthenumericalsolutionatintermediatepointsisaninterpolationofthenodalvaluesusinglinearshapefunctions.Anacceptablerepresentationofthenumericalerrorbetweenthefiniteelementsolutionandtheexactanalyticalsolutionisdefinedastheareaboundedbythetwocurves,whicharedepictedinFigure1.7,ascomparedtothetotalareaunderthecurvedescribedbytheexactsolution.Inequationform,1Ne(e)(e)error(%)=A−A100%(1.108)|A|ExactNumericalExacte=1ThenumericalpercenterrorversusthenumberoflinearelementsinthefiniteelementdomainisdepictedinFigure1.8(a).Notethatthisplotwasgeneratedbyusinguniform P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1428INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS2–210101–31010-norm2Error(%)L0–41010–1–5101051015202530354045505101520253035404550NumberofelementsNumberofelements(a)(b)FIGURE1.8:(a)Numericalpercenterrorbasedonarea.(b)ErrorbasedontheL2-normdefinitiondiscretization,meaningthatallelements(linesegments)hadexactlythesamelength.Fromthisfigure,itisobservedthatbydoublingthenumberofelements,whichisequivalenttoreducingthelengthoftheelementstohalf,thepercenterrorinthenumericalsolution,ascomparedtotheexactanalyticalsolution,isreducedbyafactorof4.ThiscanbeclearlyseenfromTable1.2,whichshowsthenumericalvalueofthecomputedpercenterrorasafunctionofthenumberoflinearelements.TABLE1.2:NumericalerrorasafunctionofthenumberoflinearelementsNUMBEROFLINEARNUMERICALPERCENTERRORBASEDONTHEELEMENTSERRORL2NORM523.48627.5×10−3105.87161.9×10−3152.60968.3×10−4201.46794.7×10−4250.93943.0×10−4300.65242.1×10−4350.47931.5×10−4400.36701.2×10−4450.29009.2×10−5500.23497.5×10−5 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS29Anotherwaytoquantifythenumericaldiscrepancybetweentwosolutions,specificallythefiniteelementsolutionVfeandtheexactanalyticalsolutionVex,istocomputetheL2norm,whichrepresentsthedistancebetweenthetwosolutions,givenbyN"#1/2!!e2Lnorm=!V!=V(e)(e)dx(1.109)2ex−Vfe2ex−Vfee=1eThenumericalerror,basedontheL2norm,wasalsocomputedandisshownplottedinFigure1.8(b).ComparingFigure1.8(a)withFigure1.8(b),itisobservedthattheerrorforbothcasesdecreasesatthesamerate.Notethatbydoublingthenumberoflinearelementsinthefiniteelementdomain,theerrorbasedontheL2normdecreasesbyafactorof4,whichwasalsothecaseforthenumericalpercenterrorcalculatedbasedontheareaboundedbythetwosolutions.Eitherwayofcomputingthenumericalerrorisacceptable.Bothmethodsprovideagoodindicationoftheaccuracyofthefiniteelementsolutionascomparedtotheexactanalyticalsolution.Besidesplottingtheelectricpotentialinthefiniteelementdomain,onemaydecidetoplottheelectricfieldasafunctionofthex-coordinate.Theexactanalyticalexpressionoftheelectricfieldisgivenby(1.7)and,asshown,itisalinearfunctionofx.Thefiniteelementsolution,basedonlinearelements,isgivenby(1.107).Asshownfromthisexpression,theelectricfieldisconstantinsideanelementandcertainlynotnecessarilycontinuousacrosselementboundaries.Figure1.9showsacomparisonbetweenthefiniteelementsolutionusingfourlinearelements(uniformdiscretization)andtheexactanalyticalexpressionoftheelectricfield.Notethattheelectricfieldisavectorquantityand,forthisproblem,ithasadirectionalongthepositivex-axis.Asstatedbefore,theelectricfieldobtainedfromthenumericalapproachisshowntobeconstantovertheelementanddiscontinuousacrosselementboundaries.Asthefiniteele-mentmeshbecomesincreasinglydenser,thenumericalsolutionapproachestheexactanalyticalsolution.1.10ONE-DIMENSIONALHIGHERORDERINTERPOLATIONFUNCTIONSIntheprevioussections,linearinterpolationfunctionswereusedforthesolutionoftheelec-trostaticBVPathand.Inthissection,weareintroducinghigherorderinterpolationfunctions,specificallyquadraticandcubic,inordertomoreaccuratelyrepresentthefiniteelementsolu-tionwithinthediscretizeddomain.Itisexpectedthatthenumericalerrorwillbesubstantiallyreducedwiththeuseofhigherorderelementsasopposedtolinearelements. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1430INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS60FEM50Exact403020100Electricfield(V/m)–10–20–30–400.00.010.020.030.040.050.060.070.08Distance(m)FIGURE1.9:ComparisonbetweentheelectricfieldproducedbytheFEMusinglinearelementsandtheexactanalyticalsolution1.10.1QuadraticElementsQuadratic,insteadoflinear,shapefunctionswillbeusedtointerpolatethesolutionofaBVPoveranelement.Alinearrepresentationofthesolutionoveranelementrequiresthevaluesoftheprimaryunknownquantityatonlytwonodes,whichcoincidewiththeendnodesoftheelement.Ontheotherhand,aquadraticrepresentationofthesolutionoveranelementrequiresthevaluesoftheprimaryunknownquantityatthreenodesinsteadofjusttwo.Twoofthesenodescoincidewiththeendnodesoftheelementwhereasthethirdonemustbeaninteriornode.Althoughthethirdnodecouldbechosenatanyinteriorpoint,themostconvenientchoiceisthemidpointoftheelement.Thus,thegeometryofaquadraticfiniteelementalongthex-axis,aswellasthegeometryofaquadraticfiniteelementalongtheξ-axis(naturalcoordinatesystem),areshowninFigures1.10(a)and1.10(b),respectively.132132xxxeee1x3x2–10+1(a)(b)FIGURE1.10:Geometryofaquadraticfiniteelement:(a)alongthex-axisand(b)alongtheξ-axis P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS31Notethattheleftmostnodeisgiventhelocalnumber1,therightmostnodeisgiventhelocalnumber2,andthemiddlenodeisgiventhelocalnumber3.Thecoordinatetransformationusedtoconvertfromthex-coordinatesystemtothenaturalξ-coordinatesystemisgivenby2x−xe3ξ=(1.110)xe−xe21wherexe+xexe=12(1.111)32Inotherwords,thecoordinatex=xemapstoξ=−1,thecoordinatex=xemapstoξ=+1,12andthecoordinatex=xemapstoξ=0.Theunknownquantity(inourcase,theelectrostatic3potential)isrepresentedoveranelementbyV(ξ)=VeN(ξ)+VeN(ξ)+VeN(ξ)(1.112)112233whereNj(ξ)forj=1,2,3arequadraticshapefunctions.ThesearealsoknownasLagrangeshapefunctions.ThequadraticshapefunctionN1(ξ)mustbe1atnode1and0attheothertwonodes.Inotherwords,N1(−1)=1N1(0)=0(1.113)N1(1)=0Therefore,thisquadraticshapefunctionhastworoots:oneatξ=0andanotheratξ=1.Thus,itmustbeoftheformN1(ξ)=cξ(ξ−1)(1.114)wherecisaconstant.Todetermineconstantc,wemustimposetheconditionN1(−1)=1(1.115)Asaresult,1c(−1)(−2)=1⇒c=(1.116)2Therefore,1N1(ξ)=ξ(ξ−1)(1.117)2 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1432INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSUsingasimilarapproach,onecanderivetheothertwoquadraticshapefunctions:N(ξ)=1ξ(1+ξ)22(1.118)N3(ξ)=(1+ξ)(1−ξ)Thesethreequadraticshapefunctionsareshownplottedbetweenξ=−1andξ=+1inFigure1.11.Thecoordinatetransformationfromthex-coordinatesystemtothenaturalξ-coordinatesystemgivenby(1.110),canbederivedfromthemappingexpressionx=xeN(ξ)+xeN(ξ)+xeN(ξ)(1.119)112233Noticeherethatthesameinterpolationfunctionsthatareusedtorepresenttheprimaryunknownquantityinthenaturalcoordinatesystemarealsousedfortherepresentationofthex-coordinate.1.01.00.80.80.60.610.420.4NN0.20.200–0.2–0.2–1–0.8–0.6–0.4–0.200.20.40.60.81–1–0.8–0.6–0.4–0.200.20.40.60.81xx(a)(b)1.00.80.63N0.40.20–1–0.8–0.6–0.4–0.200.20.40.60.81x(c)FIGURE1.11:Quadraticshapefunctions:(a)N1,(b)N2,and(c)N3 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS33ThisiscalledisoparametricrepresentationandtheelementsusedintheFEMformulationarewidelyreferredtoasisoparametricelements.Substitutingtherespectivequadraticshapefunctionsinto(1.119)yieldsxexex=1ξ(ξ−1)+2ξ(1+ξ)+xe(1+ξ)(1−ξ)322xexexe+xe=1ξ(ξ−1)+2ξ(1+ξ)+12(1+ξ)(1−ξ)222xexe=1(1−ξ)+2(1+ξ)(1.120)22xe+xexe−xe1221=+ξ22xe−xe=xe+21ξ32Solvingintermsofξ,oneobtains2x−xe3ξ=(1.121)xe−xe21whichisthecoordinatetransformationgivenin(1.110).Theelementmatrixandright-hand-sidevectoroftheBVPathandwillbeevaluatedusingquadraticelementsinSection1.11.Exercise1.1.FollowingthesameapproachusedforthederivationofN1(ξ),provethattheremainingtwoquadraticshapefunctionsaregivenby(1.118).1.10.2CubicElementsCubicrepresentationofthesolutionoveranelementrequiresfournodes.Twoofthesecoincidewiththeendnodesoftheelementandtheothertwocorrespondtointeriorpoints.Thegeometryofacubicfiniteelement,inthex-coordinatesystemandtheξ-coordinatesystem,isshowninFigure1.12.Notethatthenodes,alongboththex-andξ-axes,areequidistant.Forisoparametric13421342xxxexexexe–1–1/3+1/3+11342(a)(b)FIGURE1.12:Geometryofacubicfiniteelement:(a)alongthex-axisand(b)alongtheξ-axis P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1434INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICScubicfiniteelements,boththeprimaryunknownquantityVandthex-coordinateareexpressedaccordingtoV(ξ)=VeN(ξ)+VeN(ξ)+VeN(ξ)+VeN(ξ)(1.122)11223344x(ξ)=xeN(ξ)+xeN(ξ)+xeN(ξ)+xeN(ξ)(1.123)11223344wherethecubicshapefunctionsNj(ξ)forj=1,2,3,4canbederivedbyforcingNj(ξ)tobe1atthejthnodeand0atallothernodes.Specifically,fornode1N1(−1)=1N−1=013(1.124)N+1=013N1(+1)=0Thus,therootsofN1areξ=13ξ=−1(1.125)3ξ=−1Consequently,thecubicshapefunctionN1(ξ)canbewrittenasN(ξ)=cξ−1ξ+1133(ξ−1)(1.126)wherecisaconstant.ConstantccanbedeterminedbyenforcingtheconditionN1(−1)=1(1.127)thushaving1=c−1−1−1+1(−1−1)33⇒1=c−4−2(−2)33(1.128)⇒1=c−169⇒c=−916Therefore,thecubicshapefunctionthatcorrespondstonode1isN(ξ)=−9ξ−1ξ+1(ξ−1)(1.129)11633Similarly,thecubicshapefunctionsthatcorrespondtotheremainingthreenodesareN(ξ)=9(ξ+1)ξ−1ξ+121633N(ξ)=27(ξ+1)ξ−1(ξ−1)(1.130)3163N(ξ)=−27(ξ+1)ξ+1(ξ−1)4163 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS35110.80.80.60.610.420.4NN0.20.200–0.2–0.2–1–0.8–0.6–0.4–0.200.20.40.60.81–1–0.8–0.6–0.4–0.200.20.40.60.81xx(a)(b)1.21.2110.80.80.60.630.440.4NN0.20.200–0.2–0.2–0.4–0.4–1–0.8–0.6–0.4–0.200.20.40.60.81–1–0.8–0.6–0.4–0.200.20.40.60.81xx(c)(d)FIGURE1.13:Cubicshapefunctions:(a)N1,(b)N2,(c)N3,and(d)N4TheseareshownplottedovertheelementdomaininFigure1.13.Theelementmatrixandright-hand-sidevectorfortheBVPathandwillbeevaluatedusingcubicelementsinSection1.12.Thecoordinatetransformationfromthex-coordinatesystemtotheξ-coordinatesystemcanbederivedbysubstitutingtheexpressionsforthecubicshapefunctionsinto(1.123)andnoticingthat2xe+xexe=1233ee(1.131)2x+xxe=2143Incorporating(1.131)into(1.123)andmanipulatingtheterms,onecanshowthat2x−xe2x−xeccξ==(1.132)xe−xele21 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1436INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwherexeisthex-coordinateofthemidpointoftheelementgivenbycxe+xexe=12(1.133)c2Exercise1.2.FollowingthesameapproachusedforthederivationofN1(ξ),provethattheremainingthreecubicshapefunctionsaregivenby(1.130).Exercise1.3.Usingisoparametriccubicelements,showthatthecoordinatetransformationfromthex-coordinatesystemtotheξ-coordinatesystemisgovernedby(1.132).1.11ELEMENTMATRIXANDRIGHT-HAND-SIDEVECTORUSINGQUADRATICELEMENTSAccordingtotheweakformulationpresentedinSection1.6,theentriesoftheelementcoeffi-cientmatrixandright-hand-sidevectoraregivenbyxe2dNidNjKe=εedxfori,j=1,2,3(1.134)ijxedxdx1andxe2fe=Nρdxfori=1,2,3(1.135)iivxe1Notethattheglobalright-hand-sidevectorddoesnotcontributetotheglobalmatrixsystembecausetheonlynonzeroentriesofthisvectorbelongtothefirstandthelastrows.ThesewilleventuallybediscardedafterimposingthetwoDirichletboundaryconditionsattheendnodesofthedomain.ToobtaintheentriesoftheelementcoefficientmatrixKe,theintegralin(1.134)mustbeevaluatedeitheranalyticallyornumerically.Usingquadraticelements,itwasshownthattheexpressionthatmapsthex-coordinatetotheξ-coordinatesystemisgivenby2x−xe3ξ=(1.136)xe−xe21Takingthederivativewithrespecttoxyields22dξ=dx=dx(1.137)xe−xele21 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS37Thus,ledx=dξ(1.138)2Inaddition,usingthechainruleofdifferentiationdNidNidξ2dNi==(1.139)dxdξdxledξAsaresult,theintegralin(1.134)canbeequivalentlyexpressedintermsofthenaturalcoordinateinsteadofthex-coordinate.Doingthis,thelimitsofintegrationremainalwaysthesameforallelementsinthedomain,i.e.,from−1to+1.Thus,afterthecoordinatetransformationhastakenplace,theintegralthatisusedintheevaluationoftheentriesoftheelementcoefficientmatrixtakestheform+1ee2dNie2dNjlK=εdξijledξledξ2−1(1.140)e+12εdNidNj=dξfori,j=1,2,3le−1dξdξwhereN(ξ)=1ξ(ξ−1)12N(ξ)=1ξ(1+ξ)(1.141)22N3(ξ)=(1+ξ)(1−ξ)DifferentiatingthesethreequadraticshapefunctionswithrespecttothenaturalcoordinateξyieldsdN11=ξ−dξ2dN21=ξ+(1.142)dξ2dN3=−2ξdξSubstitutingin(1.140)fori=j=1,theentryKebecomes11e+12e2ε17εKe=ξ−dξ=(1.143)11le23le−1 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1438INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSSimilarly,εeKe=Ke=12213le2εeKe=Ke=−13313le7εeKe=(1.144)223le8εeKe=Ke=−23323le16εeKe=333leThus,theelementcoefficientmatrixisa3×3symmetricsquarematrixgivenby⎡⎤e71−8eε⎢⎥K=⎣17−8⎦(1.145)3le−8−816Theelementright-hand-sidevectorisa3×1vectorwhoseentriesareobtainedbyevaluatingthefollowingintegral:e+1elρ0fi=−Ni(ξ)dξfori=1,2,3(1.146)2−1whereρvwasreplacedby−ρ0,whichistheuniformelectronchargedensitybetweenthetwoparallelplates.SubstitutingtheshapefunctionsNi(ξ)fori=1,2,3in(1.146)andevaluatingthecorrespondingintegral,onecanobtainanalyticallytheentriesoftheelementvectorfe:e+1eelρ01lρ0f=−ξ(ξ−1)dξ=−1226−1e+1eelρ01lρ0f=−ξ(1+ξ)dξ=−(1.147)2226−1e+1eelρ02lρ0f3=−(1+ξ)(1−ξ)dξ=−2−13Asaresult,theelementright-hand-sidevectorforquadraticnodalfiniteelementsisgivenby⎧⎫⎪⎨1⎪⎬leρe0f=−1(1.148)6⎪⎩⎪⎭4Exercise1.4.Usingquadraticelements,showthattheelementcoefficientmatrixisgivenby(1.145). P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS391.12ELEMENTMATRIXANDRIGHT-HAND-SIDEVECTORUSINGCUBICELEMENTSTheprocedureusedtoderivetheelementcoefficientmatrixandright-hand-sidevectorofthedifferentialequationathandisidenticaltotheoneusedforquadraticelementsbutcertainlyabitmoreinvolved.Symbolicmathpackages,suchasMaple,canbeusedtoanalyticallyobtaintheentriesoftheelementcoefficientmatrixandright-hand-sidevector,thusavoidingunnecessarymathematicalcomplexitiesbyhand.Itcanbeshownthatthegoverningelementcoefficientmatrixandright-hand-sidevectorforcubicelementsaregivenby⎡⎤148−13−18954e⎢⎥eε⎢−1314854−189⎥K=⎢⎥(1.149)40le⎣−18954432−297⎦54−189−297432⎧⎫⎪⎪1⎪⎪⎪⎨⎪⎬leρ1e0f=−(1.150)8⎪⎪3⎪⎪⎪⎩⎪⎭3Exercise1.5.FollowthesameformulationpresentedinSection1.11toshowthattheelementcoefficientmatrixandright-hand-sidevectorforcubicelementsaregivenby(1.149)and(1.150),respectively.Toavoidtheevaluationofcomplicatedintegralsbyhand,itisrecommendedthatthesymbolicmathpackagecalledMaplebeutilized.1.13POSTPROCESSINGOFTHESOLUTION:QUADRATICELEMENTSInthecontextofFEM,aBVPisrepresentedbyasetofindependentlinearequationsthatcanbesolvednumerically,usinglinearalgebratechniques,toobtainthevaluesoftheunknownquantityatthenodesofthefiniteelementmesh.Forabetterrepresentationoftheunknownquantitiesoverthecomputationaldomain,itisinstructivethatthesequantitiesbeevaluatedatpointsotherthanthenodesofthemesh.ForthespecificBVPconsideredinthischapter,theprimaryunknownquantityistheelectrostaticpotentialwhereasthesecondaryunknownquantityistheelectricfield.Thetypeofinterpolationfunctionsusedoverasingleelementisquantraticand,therefore,theelectrostaticpotentialatanypointinsideanelementcanbeexpressedasV=VeN+VeN+VeN1122331(1.151)=Veξ(ξ−1)+Veξ(ξ+1)+Ve123(1+ξ)(1−ξ)2 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1440INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS1.0FEM(Quadraticelements)Exact0.80.60.40.20.0Electricpotential(V)–0.2–0.40.00.010.020.030.040.050.060.070.08Distance(m)FIGURE1.14:ComparisonbetweentheelectrostaticpotentialobtainedusingtheFEMandtheexactanalyticalexpression.TheFEMsolutionwasobtainedusingfourquadraticelementsandplottedbasedontenevaluationpointsperelementwhereξisgivenby(1.136).Therefore,theelectrostaticpotentialcanbeeasilyplottedintermsofthex-coordinatebyloopingthroughalltheelementsone-by-oneandevaluatingitsvalue,using(1.151),atapredefinedsetofpoints.Figure1.14illustratestheelectrostaticpotentialasafunctionofx,whichwasobtainedaftertheBVPunderconsiderationwassolvedusingfourquadraticelementsandtenevaluationpointsperelement.Itisinterestingtoobservethatthenumericalsolutionisidenticaltotheexactanalyticalsolution,andtherefore,thetwoplotsareindistinguishable.Thus,thenumericalerroriseffectivelyzerosincequadraticelementswereusedtointerpolateasolutionwhichisquadraticinnature[seeEq.(1.5)].Thiswouldnothavebeenthecaseiftheexactanalyticalsolutionwereofhigherorder.Theelectricfieldiscomputedbytakingthefirstderivativeoftheelectrostaticpotentialin(1.151)andevaluatingtheresultingexpressionatanumberofpointsalongeachelement.Specifically,theelectricfieldalongthex-directioniswrittenasdVdVdξ2dVEx=−=−=−(1.152)dxdξdxledξ P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS41Substituting(1.151)into(1.152),thex-directedelectricfieldintheregionbetweenthetwoparallelplatesusingquadraticshapefunctionscanbeexpressedas211E=−Veξ−+Veξ++Ve(1.153)xe123(−2ξ)l22whereξisgivenby(1.136).Toevaluatetheelectricfieldinthedomain,wemustloopthroughalltheelementsand,foreachelement,evaluatetheexpression(1.153)ataselectednumberofpoints.Usingtenevaluationpointsperelement,thex-directedelectricfieldbetweenthetwoparallelplatesisshownplottedinFigure1.15.Aswasthecasewiththeelectrostaticpotential,thefiniteelementsolutionisidenticaltotheexactanalyticalsolution,thusresultinginzeronumericalerror.Thereasonstemsfromthefactthattheexactexpressionoftheelectricfieldwithinthetwoparallelplatesislinearinnature[seeEq.(1.7)]and,asaresult,quadraticshapefunctionscanadequatelyrepresentlinearvariationofaquantitythatisproportionaltothefirstderivativeoftheprimaryunknownquantity.Itisworthemphasizingherethatthenumericalerrorisnotidenticaltozerobutclosetothemachineerrorofthecomputerusedtoperformthecomputations.60FEM(Quadraticelements)50Exact403020100Electricfield(V/m)–10–20–30–400.00.010.020.030.040.050.060.070.08Distance(m)FIGURE1.15:ComparisonbetweentheelectricfieldobtainedusingtheFEMandtheexactanalyticalexpression.TheFEMsolutionwasobtainedusingfourquadraticelementsandplottedbasedontenevaluationpointsperelement P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1442INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSTheeffectivenessandaccuracyofquadraticshapefunctionsinthecontextofthenodalFEM,whichisusedheretosolvea1-DBVP,canbebetterevaluatedif,forthesametypeofproblem,weconsideranonuniforminsteadofauniformchargedistributionbetweenthetwoparallelplates.Theseparationbetweentheplatesandtheboundaryconditionsaremaintainedthesameasbefore.Theonlymodificationtotheoriginalproblemistheprofileofthechargedistributionwhichisnowgivenbyx2ρv=−ρ01−(1.154)dInotherwords,insteadofhavingauniformchargedistributionintheregionbetweentheplates,wenowhaveachargedistributionthatvariesinaparabolicmannerhavinganegativemaximumvalueof−ρ0attheleftmostplateandaminimumvalueof0attherightmostplate.Thisproblemcanbesolved—asbefore—analyticallytoobtainacloseformsolutiongivenby242dρ0xdρ0V0dρ0V(x)=1−+−x+V0−(1.155)12εd12εd12εwhichisafourth-orderequation,andtherefore,cannotbeperfectlyinterpolatedbyquadraticshapefunctions.Thefiniteelementformulationforthistypeofproblemisexactlythesameastheoneoutlinedforauniformchargedistributionwiththeonlydifferencebeingtheelementright-hand-sidevector.Usingthenonuniformchargedistributiondescribedby(1.154)andafiniteelementformulationusingquadraticshapefunctions,theelementright-hand-sidevectorbecomes⎧⎫⎪⎨fe⎪⎬1fe=fe(1.156)⎪⎩2⎪⎭fe3whereρ22α2α2e0f=−++1410d23d3ρ22α2α2e0(1.157)f=−−+2410d23d3ρ24α2e0f=−+3215d23Noticethatisthelengthoftheelementgivenby=xe−xe(1.158)21 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS43andxeα=1−3(1.159)dUsingthisformulationoftheright-hand-sidevector,theBVPathandwassolvedforanonuniformchargedistributionusingquadraticshapefunctions.Thedistributionofthegov-erningelectrostaticpotentialisplottedinFigures1.16(a)and1.16(b)usingatwo-elementmeshandafour-elementmesh,respectively.Inbothfigures,theexactanalyticalsolution,givenby(1.155),isprovidedforthepurposeofcomparison.Itcanbeclearlyseenthatthenumericalsolutioncloselyapproachestheexactsolutionasthenumberofquadraticelementsinthemeshisincreased.Thesameconclusiveremarkappliestotheelectricfield.However,theelectricfieldislinearlyinterpolatedbecauseitiscomputedfromthegradientoftheelectrostaticpo-tential.Inaddition,thereisadiscontinuityoftheelectricfieldatelementboundaries.Thiswasthecasewithlinearelementsaswell.Thisdiscontinuityofthesecondaryunknownquantitytendstobecomeincreasinglysmallerasthenumberofelementsincreases.Notethattheun-derlinedweakformulationandassociatedshapefunctionsguaranteecontinuityoftheprimaryunknownquantityacrosselementsbutnotnecessarilyofthesecondaryunknownquantity.Thisisdemonstratedbyplottingtheelectricfielddistributionbetweenthetwoparallelplatesforatwo-elementmeshandafour-elementmeshusingquadraticshapefunctions.Thecorre-spondinggraphs,togetherwiththeexactanalyticalsolution,areillustratedinFigures1.17(a)and1.17(b),respectively.Onceagain,thefiniteelementsolution,forboththeelectrostaticpotentialandtheelectricfield,becomesmoreaccurateasthenumberofquadraticelementsincreases.Theaccuracyofthenumericalsolutioncanbeevaluatedbycomputingthenumericalerror,basedeitherontheareaboundedbetweenthetwocurvesorthedefinitionofL2norm,asafunctionofthenumberofquadraticelements.Thiswasdonehereusingthefirstdefinitionofnumericalerror,i.e.,basedontheareaboundedbetweenthefiniteelementsolutionandtheexactanalyticalsolution.Thisisreferredtoasthenumericalpercenterror.Aplotofthenumericalpercenterror,asfarasthecomputationoftheelectrostaticpotentialisconcerned,asafunctionofthenumberofquadraticelementsisdepictedinFigure1.18.ComparingthisgraphwiththecorrespondinggraphofFigure1.8(a),whichdepictsthepercenterrorbetweenthetwosolutionsbutusinglinearinsteadofquadraticfunctions,itisclearthattherateofconvergencehasincreasedsignificantlywiththeuseofquadraticshapefunctions.Specifically,withtheuseof20quadraticelements,thepercenterrorhasdroppeddownto10−3,whereasinthecaseofusing20linearshapefunctionsthepercenterroriscloseto1.TherateofconvergenceusingquadraticshapefunctionscanbemoreclearlyseenfromTable1.3,whichtabulatesthedatathatcorrespondstothegraphillustratedinFigure1.18.Fromthistable,onecandeducethat P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1444INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS1.0FEM(2Quadraticelements)Exact0.80.60.4Electricpotential(V)0.20.00.00.010.020.030.040.050.060.070.08Distance(m)(a)1.0FEM(4Quadraticelements)Exact0.80.60.4Electricpotential(V)0.20.00.00.010.020.030.040.050.060.070.08Distance(m)(b)FIGURE1.16:Electrostaticpotentialforanonuniformchargedistributionintheregionbetweentheplates.TheFEMusesquadraticshapefunctions:(a)two-elementmeshand(b)four-elementmesh P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS4540FEM(2Quadraticelements)Exact3530252015Electricfield(V/m)10500.00.010.020.030.040.050.060.070.08Distance(m)(a)40FEM(4Quadraticelements)Exact3530252015Electricfield(V/m)10500.00.010.020.030.040.050.060.070.08Distance(m)(b)FIGURE1.17:Electricfieldforanonuniformchargedistributionintheregionbetweentheplates.TheFEMusesquadraticshapefunctions:(a)two-elementmeshand(b)four-elementmesh P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1446INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS10Quadraticelements52152–11052Percenterror(%)–21052–31002468101214161820NumberofelementsFIGURE1.18:Numericalpercenterrorbetweenthefiniteelementsolutionandtheexactanalyticalsolution.ThecomputationswereperformedbasedontheelectrostaticpotentialTABLE1.3:NumericalpercenterrorasafunctionofthenumberofquadraticelementsNUMBEROFQUADRATICELEMENTSNUMERICALPERCENTERROR21.474640.184460.054680.0230100.0118120.0068140.0043160.0029180.0020200.0015 P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARY-VALUEPROBLEMS47thepercenterrordropsbyafactorof8ifthenumberofquadraticelementsisdoubled,whereasinthecaseoflinearshapefunctions,thisfactorwas4.Exercise1.6.Showthat(1.155)isindeedtheanalyticalsolutionoftheBVPathand.Rememberthattheelectronchargedistributionintheregionbetweenthetwoparallelplateshaschangedfromaconstantprofiletoaparabolicprofilegivenby(1.154).Exercise1.7.UsingquadraticelementsandtheformulationpresentedinSection1.11forthederivationoftheelementright-hand-sidevectorfe,showthatthelatterisgivenby(1.156)–(1.159).1.14POSTPROCESSINGOFTHESOLUTION:CUBICELEMENTSThemajorstepsinvolvedinthepostprocessingofthefiniteelementsolutionobtainedusingcubicelementsareidenticaltothestepsinvolvedinthepostprocessingofthefiniteelementsolutionobtainedusingquadraticelements.Theonlydifferenceistheorderofshapefunctionsusedtointerpolatetheprimaryunknownvariableoveranelement.Toavoidunnecessaryrepetitions,wewillrestrictourdiscussiononthedevelopmentofthegoverningexpansionequationsdescribingtheprimaryandsecondaryunknownvariablesoveranelementafterthefiniteelementsolutionhasbeenobtained.Itwillbeleftasanexerciseforthereadertorepeat,usingcubicelements,thesameprocedureaswasdoneintheprevioussection,usingquadraticelements,forthecomputationofthenumericalpercenterrorbetweentheexactandthefiniteelementsolutionsandtocomparetheresultwithFigure1.18orTable1.3.SolvingtheglobalmatrixsystemafterimposingthegoverningDirichletboundarycondi-tions,asetofvaluesfortheprimaryunknownvariable(i.e.,electrostaticpotential)isobtained.Thesevaluescorrespondtotheglobalnodesofthefiniteelementmesh.Forsomeonetoplottheprimaryunknownvariablewithgoodenoughresolutioninthediscretizeddomain,itisnec-essarythattheprimaryunknownvariablebeevaluatedatmultiplepointsalongeachelement.Whenusingcubicshapefunctions,theelectrostaticpotentialisgivenbyV=VeN(x)+VeN(x)+VeN(x)+VeN(x)(1.160)11223344whereVefori=1,2,3,4arethevaluesoftheelectrostaticpotentialatthefournodesoftheielement,andNi(x)fori=1,2,3,4arethecorrespondingcubicshapefunctionsintermsofx.Thecubicshapefunctionsexpressedintermsofxareobtainedbysubstituting(1.132)intothesetoffourequationsgivenby(1.129)and(1.130).Aswasdonewithquadraticelements,theelectricfieldoveracubicelementiscomputedbytakingthefirstderivativeof(1.160)andevaluatingtheresultingexpressionataselectionofpointsalongtheelement.Thisderivativecanbeevaluatedmoreconvenientlybyimplementing P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1448INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSthechainruleofdifferentiation4dVdVdξ2dV2edNiEx=−dx=−dξdx=−ledξ=−leVidξ(1.161)i=1whereNifori=1,2,3,4aregivenby(1.129)and(1.130).Fromthisdiscussion,itbecomesclearnowthattheelectricfield,whichisthesecondaryunknownvariable,isquadraticoveracubicelementwhereastheelectrostaticpotential,whichistheprimaryunknownvariable,iscubic.1.15SOFTWAREThreeMatlabcodeswerewritteninordertosolvetheelectrostaticBVPdiscussedinthischapter.Thefirstcode,FEM1DL,useslinearelementstosolvethePoisson’sequationinonedimensionassumingthatthechargedistributionbetweenthetwoparallelplatesisuniform.Dirichletboundaryconditionsareimposedonthetwoplates.Thecodecomputesboththeelectrostaticpotentialandtheelectricfieldwithinthedomainofinterest.Thesecondcode,FEM1DQ,usesquadraticelementstosolvethesameexactproblem.Thethirdcode,FEM1Dqnucd,usesquadraticelementstosolvethesameproblembutwithanonuniform,insteadofuniform,chargedistribution.Allthreecodesarecapableofcomputingthenumericalerrorassociatedwiththefiniteelementsolution.Certainparameterssuchasmeshsizeandseparationofplatescanbemodifiedbytheuser.AllthreeMatlabcodescanbedownloadedfromthepublisher’sURL:www.morganclaypool.com/page/polycarpou.Exercise1.8.WriteafiniteelementcodeinMatlabtosolvethesameBVPconsideredinthischapterbutusingcubicelementsinsteadoflinearorquadraticelements.Assumethattheelectronchargedistributionbetweentheplatesisgivenby(1.154).Computeandplotthenumericalerror,eitherasapercentageorusingtheL2-normdefinition,andcomparewithTables1.2and1.3.REFERENCES[1]O.C.Zienkiewicz,TheFiniteElementMethod.NewYork:McGraw-Hill,1977.[2]P.P.SilvesterandR.L.Ferrari,FiniteElementsforElectricalEngineers,2nded.London:CambridgeUniversityPress,1990.[3]J.Jin,TheFiniteElementMethodinElectromagnetics,2nded.NewYork:Wiley-IEEEPress,2002.[4]J.L.Volakis,A.Chatterjee,andL.C.Kempel,FiniteElementMethodElectromagnetics:Antennas,MicrowaveCircuits,andScatteringApplications.NewYork:Wiley-IEEEPress,1998. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:14ONE-DIMENSIONALBOUNDARYVALUEPROBLEMS49[5]G.Pelosi,R.Coccioli,andS.Selleri,QuickFiniteElementsForElectromagneticWaves.Boston:ArtechHouse,1998.[6]J.N.Reddy,AppliedFunctionalAnalysisandVariationalMethodsinEngineering.NewYork:McGraw-Hill,1986.[7]S.G.Mikhlin,VariationalMethodsinMathematicalPhysics.NewYork:Macmillan,1964.[8]M.N.O.Sadiku,NumericalTechniquesinElectromagnetics,2nded.BocaRaton:CRCPress,2001.[9]D.K.Cheng,FundamentalsofEngineeringElectromagnetics.NewYork:Addison-Wesley,1993.[10]C.H.Edwards,Jr.,andD.E.Penney,ElementaryLinearAlgebra.NewJersey:Prentice-Hall,1988. P1:IML/FFXP2:IML/FFXQC:IML/FFXT1:IMLMOBK021-01MOBK021-Polycarpou.clsApril29,200619:1450 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2851CHAPTER2Two-DimensionalBoundary-ValueProblems2.1INTRODUCTIONInthischapter,thenodalFEMwillbeappliedtoageneric2-DBVPinelectromagnetics.Suchproblemsusuallyinvolveasecond-orderdifferentialequationofasingledependentvariablethatissubjecttoasetofboundaryconditions.TheseboundaryconditionscouldbeoftheDirichlettype,theNeumanntype,orthemixedtype.Thedomainoftheproblemisa2-Dgeometrywithanarbitraryshape.Thus,anaccuraterepresentationofthedomaininthecontextoftheFEMpresumesdiscretizationofthedomainusingthemostappropriateshapeofbasicelementscalledthefiniteelements.Themostcommonlyusedfiniteelementsintwodimensionsarethetriangularandquadrilateralelements.Bothtypesofelementswillbeusedinthischaptertosolve2-DBVPsinelectromagnetics.Furthermore,higherorderelementswillbeusedtoillustratetheimprovementinaccuracyofthenumericalsolutionwithoutnecessarilyincreasingthenumberofelementsinthefiniteelementmesh.ThemajorstepsinvolvedintheapplicationoftheFEMforthesolutionof2-DBVPsareidenticaltothoseinvolvedinthesolutionof1-Dproblems.Specifically,aproperapplicationoftheFEMforthesolutionof2-DBVPsmustinvolvethefollowingmajorsteps:a)Discretizationofthe2-Ddomain.b)Derivationoftheweakformulationofthegoverningdifferentialequation.c)Properchoiceofinterpolationfunctions.d)Derivationoftheelementmatricesandvectors.e)Assemblyoftheglobalmatrixsystem.f)Impositionofboundaryconditions.g)Solutionoftheglobalmatrixsystem.h)Postprocessingoftheresults.Problemdefinition:ABVPcharacterizedbyagenericformofasecond-orderpartialdifferentialequationwillbeconsideredinthischaptertoillustratethemajorstepsinvolvedina2-Dnodal P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2852INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFEM.Thisgenerictypepartialdifferentialequationcanbeexpressedas∂∂u∂∂uαx+αy+βu=g(2.1)∂x∂x∂y∂ywhereαx,αy,β,andgareconstantstobedefinedbythespecificapplicationanduistheprimaryunknownquantity.Poisson’sequation,foralinearandisotropicmedium,isgivenby∇(ε∇V)=−ρv(2.2)Ina2-Dspace,(2.2)isoftenwrittenas∂∂V∂∂Vε+ε=−ρv(2.3)∂x∂x∂y∂yEquation(2.3)isaspecialcaseofthegenericformgivenby(2.1).Comparing(2.1)with(2.3),itcanbeeasilyrealizedthatthesetwopartialdifferentialequationswouldbeidenticalifu=Vαx=αy=ε(2.4)β=0g=−ρvConsequently,the2-DPoisson’sequation,whichiswidelyusedtosolveelectrostaticproblems,isaspecialcaseof(2.1).ThesetofboundaryconditionscouldbeeitherofDirichlettypeu=u0on1(2.5)ormixedtype∂u∂uαxaˆx+αyaˆy·aˆn+γu=qon2(2.6)∂x∂ywhereaˆnistheunitvectornormaltotheboundary2andγ,qareconstantstobedefined.2.2DOMAINDISCRETIZATIONThedomainofa2-DBVPusuallyhasanirregularshape,asshowninFigure2.1(a).UsingtheFEM,thefirststepistoaccuratelyrepresentthephysicaldomainoftheproblembyasetofbasicshapescalledthefiniteelements.Theuseofarectangle,forexample,asabasicfiniteelementtodiscretizeanirregulardomainiscertainlythesimplestbutnotthemostsuitablechoicebecauseanassemblyofrectanglescannotaccuratelyrepresentthearbitrarygeometricalshapeofthedomain.Insuchacase,thediscretizationerrorissignificant,asshowninFigure2.1(b),althoughittendstodecreaseasthesizeofrectanglesinthedomainbecomessmaller.Ontheotherhand,ifatriangleisusedinsteadofarectangleasthebasicelementforthemeshingof P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS53Discretizationerror(a)(b)FIGURE2.1:(a)Irregular2-Ddomain.(b)Finiteelementmeshusingrectangularelementsthe2-Ddomain,thediscretizationerrorwouldbeeffectivelymuchsmaller.ThisisillustratedgraphicallyinFigure2.2(a).Thequadrilateralisanotherbasicelementthatiscommonlyusedin2-Dfiniteelementanalysis.AcoarsemeshoftheirregulardomainusingquadrilateralelementsisshowninFig-ure2.2(b).Aswasthecasewiththetriangularelement,thequadrilateralelementresultsinasmallerdiscretizationerrorthantheonecausedbytheuseoftherectangularelement.Notethattherearecertainadvantagesinusingtriangularelementsascomparedtoquadrilateralel-ements,andtheseadvantagesbecomeincreasinglyimportantwhenusingvectorelementstosolveelectromagneticproblems.SincethescopeofthisbookistohelpthereaderunderstandthebasicsoftheFEMandnottoapplythemethodtoadvancedtopicsinelectromagnetics,ourDiscretizationerror(a)(b)FIGURE2.2:(a)Finiteelementmeshusingtriangularelements.(b)Finiteelementmeshusingquadri-lateralelements P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2854INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSdiscussionwillberestrictedtonodalelements,only.Forvectorelementsandtheirapplicationtoelectromagneticproblems,thereaderisreferredtotheliterature[1–5].Meshgenerationhascertainrulesthatmustbefollowedatalltimes:•Foratriangularmesh,theshapeoftrianglesmustbeclosetoequilateral.•Foraquadrilateralmesh,theshapeofquadrilateralsmustbeclosetosquare.•Nodesmustappearatsourcepoints.•Thefiniteelementmeshmustaccuratelyrepresentthegeometricaldomainoftheproblem.•Inregionswherethesolutionisexpectedtohavelargevariations,theelementsmustbesufficientlysmall.•Avoidelementswithverylargeaspectratios,i.e.,theratioofthelargestsidetothesmallestside.•Numberthenodesinascendingorderstartingfrom1.Thenumberingofthenodesdirectlyaffectsthebandwidthoftheglobalmatrix.•Theremustbenooverlapofelements.•Neighboringelementsmustshareacommonedge.•Aninteriornode(nonboundarynode)mustbelongtoatleastthreeelements.2.3INTERPOLATIONFUNCTIONSProperinterpolationfunctionsmustbedevelopedfortriangularandquadrilateralelementssincetheyarebothwidelyusedinthediscretizationof2-Ddomains.AswasindicatedinChapter1,theseinterpolationfunctionsmustsatisfycertainkeyrequirements.First,theymustguaranteecontinuityoftheprimaryunknownquantityacrossinterelementboundaries.Second,theymustbeatleastoncedifferentiablesincethegoverningdifferentialequationathandisofsecondorderand,third,theymustbecompletepolynomialstoprovidesufficientrepresentationofthesolution’sbehaviorinthefiniteelementdomain.Initially,wewillconcentrateonthedevelopmentofinterpolationfunctionsbasedonlinear/bilinearelementsand,then,moveontohigherorderelementsbasedontheLagrangepolynomials.2.3.1LinearTriangularElementAlineartriangularelementinthexy-planeisillustratedinFigure2.3(a).Thetrianglecon-sistsofthreeverticeswhichcorrespondtothethreenodesoftheelement.Thenodesarelocallynumberedinacounter-clockwisedirectiontoavoidhavinganegativeareausingthe P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS55yh33(0,1)1212xx(0,0)(1,0)(a)(b)FIGURE2.3:(a)Lineartriangularelementinthexy-plane.(b)Lineartriangularelement(masterelement)intheξη-planeJacobian1definition.Alinearinterpolationfunctionspanningatrianglemustbelinearintwoorthogonaldirections.Thesecouldbetheorthogonalaxesdefinedbythenaturalcoordinatesξandη.Thus,atriangleofarbitraryshape,suchastheoneshowninFigure2.3(a),couldbemappedtothemastertriangle,showninFigure2.3(b),whichliesonthenaturalcoordinatesystem.Eachlinearinterpolationfunctioncorrespondstoatrianglenode.DenotingthethreeinterpolationfunctionsbyN1(ξ,η),N2(ξ,η),andN3(ξ,η),theseareassignedtonodes1,2,and3,respectively.Usingaparallelargumenttothe1-Dcase,theshapefunctionN1(ξ,η)mustbe1atnode1and0attheothertwonodes,namelynodes2and3.Thus,startingfromalinearrepresentationofshapefunctionN1(ξ,η)N1(ξ,η)=c1+c2ξ+c3η(2.7)andusingtheconditionsthatAtnode1:ξ=0,η=0⇒N1(0,0)=c1=1Atnode2:ξ=1,η=0⇒N1(1,0)=1+c2+0=0⇒c2=−1(2.8)Atnode3:ξ=0,η=1⇒N1(0,1)=1+0+c3=0⇒c3=−1onecandeducethattheinterpolationfunctionthatcorrespondstonode1isgivenbyN1(ξ,η)=1−ξ−η(2.9)1ItwillbeshownlaterinthischapterthattheJacobianisdirectlyrelatedtotheareaofthetrianglewhichisapositivequantity.Ifthelocalnodesofthetrianglearenumberedinacounter-clockwisedirection,thentheJacobiancomestobepositiveand,therefore,itcanbederiveddirectlyfromtheareaofthetriangle. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2856INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSSimilarly,N2(ξ,η)mustbe1atnode2and0atnodes1and3.AlinearrepresentationofN2(ξ,η)isN2(ξ,η)=c1+c2ξ+c3η(2.10)Imposingtheaboveconditions,wehaveAtnode1:ξ=0,η=0⇒N2(0,0)=c1=0Atnode2:ξ=1,η=0⇒N2(1,0)=0+c2+0=1⇒c2=1(2.11)Atnode3:ξ=0,η=1⇒N2(0,1)=0+0+c3=0⇒c3=0Thus,N2(ξ,η)=ξ(2.12)Finally,N3(ξ,η)mustbe1atnode3and0attheothertwonodesofthemastertriangle.Ingeneral,N3(ξ,η)canbeexpressedasN3(ξ,η)=c1+c2ξ+c3η(2.13)Toobtainconstantsc1,c2,andc3,theaboveconditionsmustbeimposed,i.e.,Atnode1:ξ=0,η=0⇒N3(0,0)=c1=0Atnode2:ξ=1,η=0⇒N3(1,0)=c2=0(2.14)Atnode3:ξ=0,η=1⇒N3(0,1)=c3=1ThefinalformofN3(ξ,η)isthereforegivenbyN3(ξ,η)=η(2.15)ThesethreeinterpolationfunctionsareplottedinFigure2.4.Itisimportanttoempha-sizeherethattheseinterpolationfunctionsarenotlinearlyindependent.OnlyN2andN3areindependent;N1isalinearcombinationofN2andN3;i.e.,N1=1−N2−N3(2.16)whichcanalsobewrittenasN1+N2+N3=1(2.17)Thesetriangle-basedlinearinterpolationfunctionscanalsobewrittenintermsofareacoordinates.Consideranarbitrarypoint(ξ,η)insidethemastertriangleshowninFigure2.5.Connectingallthreeverticestotheinteriorpoint(ξ,η),threesubtrianglesareformedwithrespectiveareasA1,A2,andA3.NotethatA1correspondstothesubtriangleoppositetolocal P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS57hhN133N2111122xxh1N3312xFIGURE2.4:Lineartriangularinterpolationfunctionsh3(x,)hA1A221xA3FIGURE2.5:Areacoordinates P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2858INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSnode1.Asimilarnotationappliestotheothertwosubtriangles.Basedonthisnotation,thelinearinterpolationfunctionsN1,N2,andN3canbeexpressedasA1N1=AA2N2=(2.18)AA3N3=AwhereAistheareaofthemastertriangle.NoticethatA1A2A3N1+N2+N3=++AAAA1+A2+A3=(2.19)AA==1AThissetoftriangularbasisfunctionsareusedintheFEMtointerpolatetheprimaryunknownquantityintheinteriorofanelement.Incaseofusinglineartriangularelementstodiscretizetheproblemdomain,theprimaryunknownquantity—letussayu—insideanelementcanbeexpressedas3u=ueN+ueN+ueN=ueN(2.20)112233iii=1whereue,ue,anduearethenodalvaluesoftheprimaryunknownquantityatthethreevertices123ofthetriangle.Forisoparametricelements,thesameshapefunctionsusedtointerpolatetheprimaryunknownquantityinsideanelementarealsousedtointerpolatethespacecoordinatesxandy.Inotherwords,3x=xeN+xeN+xeN=xeN112233iii=13y=yeN+yeN+yeN=yeN(2.21)112233iii=1SubstitutingN1,N2,andN3into(2.21)yieldsx=xe+x¯ξ+x¯η12131y=ye+y¯ξ+y¯η(2.22)12131 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS59wherex¯=xe−xe2121x¯=xe−xe3131y¯=ye−ye(2.23)2121y¯=ye−ye31312.3.2BilinearQuadrilateralElementAbilinearquadrilateralelementinthexy-planeisshowninFigure2.6(a).Thequadrilateralhasfourlocalnodeswhicharenumberedinacounter-clockwisedirection.Knowingthesolutionatthefournodesoftheelement,theprimaryunknownquantitycanbeevaluatedatanypointinsidetheelementbyusingtheappropriateinterpolationfunctions.Inthissection,wewillconstructbilinearinterpolationfunctionsforthemasterquadrilateralelement.Anisoparametricrepresentationwillbeusedtotransformafunctionfromthenaturalcoordinatesystemtothexy-coordinatesystemandviceversa.Themasterelement,whichisdefinedintheξη-coordinatesystem(naturalcoordinatesystem)hasasquareshapeandisdepictedinFigure2.6(b).Agenericbilinearinterpolationfunctionforlocalnode1spanningthegeometricaldomainofthemasterquadrilateralelementhastheformN1(ξ,η)=c1+c2ξ+c3η+c4ξη(2.24)AccordingtothepropertiesofLagrangepolynomials,1atnode1N1=(2.25)0atallothernodesh3443(–1,1)(1,1)yx1(–1,–1)(1,–1)212xFIGURE2.6:(a)Quadrilateralelementinthexy-plane.(b)Quadrilateralmasterelementintheξη-plane P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2860INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSApplyingtheseconditionsatthefournodesofthemasterquadrilateralelement,theresultisasystemoffourequationswithfourunknowns,theunknownsbeingthefourconstantsof(2.24).N1(−1,−1)=c1−c2−c3+c4=1N1(1,−1)=c1+c2−c3−c4=0N1(1,1)=c1+c2+c3+c4=0(2.26)N1(−1,1)=c1−c2+c3−c4=0Thissystemofequationscanbesolvedinastraightforwardmannertoobtaintheconstants.Itcanbeshownthatc=1,c=−1,c=−1,c=1(2.27)14243444Thus,theinterpolationfunctionthatcorrespondstonode1isgivenbyN(ξ,η)=1(1−ξ−η+ξη)(2.28)14whichcanalsobewrittenasN(ξ,η)=1(1−ξ)(1−η)(2.29)14Usingasimilarapproach,onecanconstructallfourbilinearinterpolationfunctionsspanningaquadrilateralelement:N(ξ,η)=1(1−ξ)(1−η)14N(ξ,η)=1(1+ξ)(1−η)24N(ξ,η)=1(1+ξ)(1+η)(2.30)34N(ξ,η)=1(1−ξ)(1+η)44Assuminganisoparametricquadrilateralelement,theprimaryunknownquantityandthexandyspacecoordinatesinsideanelementcanbeexpressedintermsofthesefourbasisfunctionsgivenby(2.30).Inotherwords,4u=ueN+ueN+ueN+ueN=ueN(2.31)11223344iii=1and4x=xeN+xeN+xeN+xeN=xeN11223344iii=14y=yeN+yeN+yeN+yeN=yeN(2.32)11223344iii=1 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS61Exercise2.1.Showthattheshapefunctionsgoverningabilinearquadrilateralelementaregivenby(2.30).UsethesameapproachfollowedforthederivationofN1(ξ,η).2.4THEMETHODOFWEIGHTEDRESIDUAL:THEGALERKINAPPROACHThegeneric2-DBVPconsideredatthebeginningofthischapterischaracterizedbyasecond-orderpartialdifferentialequationgivenby∂∂u∂∂uαx+αy+βu=g(2.33)∂x∂x∂y∂ywhereαx,αy,β,andgareconstants.Theweakformulationofthisproblemcanbeobtainedbyfirstconstructingtheweightedresidualof(2.33)forasingleelementwithdomaine.Theelementresidualisformedbymovingtheright-handsideof(2.33)totheleft-handside:∂∂u∂∂ure=α+α+βu−g(2.34)xy∂x∂x∂y∂yThiselementresidualisideallyzero,providedthatthenumericalsolutionutobeobtainedisidenticaltotheexactsolution.However,thisisnotthecase,andtherefore,theelementresidualreis,ingeneral,nonzero.Ourobjectiveistominimizethiselementresidualinaweightedsense.Toachievethis,wemustfirstmultiplyrewithaweightfunctionw,thenintegratetheresultovertheareaoftheelement,andfinally,settheintegraltozero.∂∂u∂∂uwαx+αy+βu−gdxdy=0(2.35)e∂x∂x∂y∂yIntroducingtheidentity∂∂u∂w∂u∂∂uwαx=αx+wαx(2.36)∂x∂x∂x∂x∂x∂xwhichcanberearrangedas∂∂u∂∂u∂w∂uwαx=wαx−αx∂x∂x∂x∂x∂x∂x∂∂u∂w∂u=wαx−αx(2.37)∂x∂x∂x∂xandsubstitutingthelatterinto(2.35)yields∂∂u∂∂u∂w∂u∂w∂uwαx+wαydxdy−αx+αydxdye∂x∂x∂y∂ye∂x∂x∂y∂y+βwudxdy=wgdxdy(2.38)ee P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2862INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSThen,weimplementtheGreen’stheoremwhichstatesthattheareaintegralofthedivergenceofavectorquantityequalstothetotaloutwardfluxofthevectorquantitythroughthecontourthatboundsthearea.Inequationform,∇t·AdA=A·aˆnd(2.39)eeorsimply∂Ax∂Ay+dxdy=aˆxAx+aˆyAy·aˆnd(2.40)e∂x∂yewhereAisthevectorquantityofinterestandaˆnistheoutwardunitvectorthatisnormaltotheboundaryoftheelement.Thecontourintegralin(2.40)mustbeevaluatedalongtheperipheryoftheelementinacounter-clockwisedirection.Comparingthefirstintegralof(2.40)withthefirstintegralof(2.38),itbecomesquiteclearthat∂uAx=wαx(2.41)∂xand∂uAy=wαy(2.42)∂yBydefiningthenormalunitvectorasaˆn=aˆxnx+aˆyny(2.43)andapplyingtheGreen’stheoremtothefirstintegralof(2.38),thelatterbecomes∂∂u∂∂u∂u∂uwαx+wαydxdy=wαxnx+αynyd(2.44)e∂x∂x∂y∂ye∂x∂ySubstitutingthisresultinto(2.38),theweakformofthedifferentialequationreducesto∂w∂u∂w∂u−αx+αydxdy+βwudxdy=wgdxdye∂x∂x∂y∂yee∂u∂u−wαxnx+αynyd(2.45)e∂x∂yAccordingtotheGalerkinapproach,theweightfunctionwmustbelongtothesamesetofshapefunctionsthatareusedtointerpolatetheprimaryunknownquantitywhich,inourcase,isu.Intheprevioussection,itwasshownthattheprimaryunknownquantityisinterpolated P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS63usingasetofLagrangepolynomials.Thus,nu=ueN(2.46)jjj=1whereNj’sarethecorrespondingshapefunctionsbasedonLagrangepolynomialsandnisthenumberoflocalnodesperelement.Substituting(2.46)into(2.45),andsettingw=Nifori=1,2,...,n(2.47)theweakformofthegoverningdifferentialequationisdiscretized:⎡⎤nn∂ueN∂ueN⎢jjjj⎥⎢⎢∂Nij=1∂Nij=1⎥⎥−⎢αx+αy⎥dxdye⎣∂x∂x∂y∂y⎦n∂u∂u+βNueNdxdy=Ngdxdy−Nαn+αnd,ijjiixxyyej=1ee∂x∂yfori=1,2,...,n(2.48)Noticethattheprimaryunknownquantityuinthecontourintegralof(2.48)hasnotbeenreplacedbythesetofinterpolationfunctionsgivenby(2.46).Thisintegralwillbetreatedseparatelyabitlater.Equation(2.48)canalsobewritteninthefollowingform:nn∂Nie∂Nj∂Nie∂Nj−αxuj+αyujdxdye∂xj=1∂x∂yj=1∂yn∂u∂u+βNueNdxdy=Ngdxdy−Nαn+αnd,ijjiixxyyej=1ee∂x∂yfori=1,2,...,n(2.49)Thisequationcanbeconvenientlyexpressedinamatrixformgivenby⎡⎤⎧⎫⎡⎤⎧⎫⎧⎫⎧⎫MeMe···Me⎪⎪ue⎪⎪TeTe···Te⎪⎪ue⎪⎪⎪⎪fe⎪⎪⎪⎪pe⎪⎪⎢11121n⎥⎪⎪1⎪⎪⎢11121n⎥⎪⎪1⎪⎪⎪⎪1⎪⎪⎪⎪1⎪⎪⎢MeMe···Me⎥⎨ue⎬⎢TeTe···Te⎥⎨ue⎬⎨fe⎬⎨pe⎬⎢21222n⎥2+⎢21222n⎥2=2+2⎢⎣............⎥⎦⎪⎪...⎪⎪⎢⎣............⎥⎦⎪⎪...⎪⎪⎪⎪...⎪⎪⎪⎪...⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪MeMe···Mnne⎩une⎭TeTe···Tnne⎩une⎭⎩fne⎭⎩pne⎭n1n2n1n2(2.50) P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2864INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwheree∂Ni∂Nj∂Ni∂NjMij=−αx+αydxdy(2.51)e∂x∂x∂y∂yTe=βNNdxdy(2.52)ijijefe=Ngdxdy(2.53)iieand∂u∂upe=−Nαn+αnd(2.54)iixxyye∂x∂yInamorecompactform,thematrixsystemin(2.50)canbeexpressedas⎡⎤⎧⎫⎧⎫KeKe···Ke⎪⎪ue⎪⎪⎪⎪be⎪⎪⎢11121n⎥⎪⎪1⎪⎪⎪⎪1⎪⎪⎢KeKe···Ke⎥⎨ue⎬⎨be⎬⎢21222n⎥2=2(2.55)⎢⎣............⎥⎦⎪⎪...⎪⎪⎪⎪...⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪KeKe···Knne⎩une⎭⎩bne⎭n1n2whereKe=Me+Teijijijbe=fe+pe(2.56)iiiNoticethatthecontourintegralin(2.54)mustbeevaluatedalongtheclosedboundaryofeverysingleelementinthedomain.Forexample,ifthefiniteelementmeshconsistsoflineartriangularelements,thiscontourintegralmustbeevaluatedalongthethreeedgesofeachtriangleinacounter-clockwisedirection.However,itisimportanttorealizethatanonboundaryedgebelongstotwoneighboringtriangles,asshowninFigure2.7(a).Asaresultofthisobservation,evaluatingthelineintegralin(2.54)forelemente1[seeFigure2.7(b)]alongtheedgefromnode1tonode2givesexactlythesameresult,butoppositesign,withevaluatingthesamelineintegralforelemente2alongtheedgefromnode3tonode1.Thereasonfortheoppositesignstemsfromthefactthatthetwooutwardunitvectorsnormaltothecommonedgeofthetwoneighboringtriangles,asshowninFigure2.7(b),pointinoppositedirections.Inotherwords,aˆn1=−aˆn2(2.57) P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS6523e2a2n2e2e11ae1n113(a)(b)FIGURE2.7:(a)Interioredgesharedbytwoneighboringtriangles.(b)Theoutwardunitvectorsnormaltothecommonedgepointinoppositedirectionswhereaˆ=aˆn(e1)+aˆn(e1)n1xxyyaˆ=aˆn(e2)+aˆn(e2)(2.58)n2xxyyThus,n(e1)=−n(e2)(2.59)xxandn(e1)=−n(e2)(2.60)yyTogiveanexample,letuscomputethecontributiontotheglobalentryoftheright-hand-sidevectorpbylocalnode1ofelemente1andlocalnode1ofelemente2astheintegrationin(2.54)isevaluatedalongtheedgecommontothetwotriangles.Denotingthiscontributiontothespecificglobalentrybyp1−1,wehavep=−N(e1)α∂un(e1)+α∂un(e1)d−N(e2)α∂un(e2)+α∂un(e2)d1−11xxyy1xxyy1→2∂x∂y3→1∂x∂y(2.61)Substituting(2.59)and(2.60)intothefirstintegralof(2.61),andusingthefactthatN(e1)=N(e2)(2.62)11 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2866INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSalongthepathofintegration,whichisthecommonedge,itbecomesevidentthatthetwointegralsareequalinmagnitudebutoppositeinsign,thuscancelingeachotherout.Noticethatthelineintegralsin(2.61)areevaluatedalongthecommonedgeand,therefore,theremainingtermsinvolvedintheintegrandof(2.54)areequalforbothtriangles.Similarly,thecontributiontotheglobalentryoftheright-hand-sidevectorpbylocalnode2ofelemente1andlocalnode3ofelemente2astheintegrationisevaluatedalongthecommonedgeisgivenbyp=−N(e1)α∂un(e1)+α∂un(e1)d−N(e2)α∂un(e2)+α∂un(e2)d2−32xxyy3xxyy1→2∂x∂y3→1∂x∂y(2.63)Again,substituting(2.59)and(2.60)intothefirstintegralof(2.63),andusingthefactthatN(e1)=N(e2)(2.64)23alongthepathofintegration,itisevidentthatthetwointegralscanceleachotherout.Conse-quently,thecontributionofthelineintegralin(2.54)totheglobalright-hand-sidevectorpiszeroforallinterioredges.Itisnonzeroonlyforedgesthatbelongtothedomainboundary,where=1∪2(2.65)Forboundaryedgesthatbelongto1,whereaDirichletboundaryconditionistobeimposed,thecontributionofthelineintegralin(2.54)willbediscarded.Thus,theonlycontributionofthelineintegralin(2.54)isattributedonlytoboundaryedgesthatresideon2.Rememberthat2ischaracterizedbyamixedboundaryconditionoftheform∂u∂uαxnx+αyny+γu=qon2(2.66)∂x∂yRearrangingthetermsin(2.66)andsubstitutinginto(2.54),thelineintegralbecomespe=−Nq−γu)d(2.67)ii(2Thelineintegralin(2.67)existsonlyforboundaryelements,i.e.,elementsthathaveoneormoreedgesontheouterboundaryofthefiniteelementdomain.Specifically,itmustbeevaluatedonlyalongboundaryedgesthatresideon2.Forinterioredges,aswasstatedbefore,itscontributioniszero.2.5EVALUATIONOFELEMENTMATRICESANDVECTORSThemainpurposeofthissectionistoderiveanalyticallytheexpressionsfortheentriesofallelementmatricesandvectorsthatarepresentinthelinearsystemofequationsgivenby(2.50). P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS67TheseentryvaluesaredependentonthetypeandorderofinterpolationfunctionsusedintheFEM.Inthissection,wearegoingtoconsidertwotypesofinterpolationfunctions:oneforthelineartriangularelement,andanotherforthebilinearquadrilateralelement.HigherorderelementswillbeconsideredseparatelyinSection2.11.2.5.1LinearTriangularElementsWebegintheevaluationofelementmatricesandvectorswithmatrixMewhoseentriesaregiven,accordingto(2.51),bye∂Ni∂Nj∂Ni∂NjMij=−αx+αydxdy(2.68)e∂x∂x∂y∂ywhereαxandαyareconstants.Atthispoint,itisimportanttoremindthereaderthatthegoverninginterpolationfunctionsforlineartriangularelementsaregivenbyN1=1−ξ−ηN2=ξN3=η(2.69)andthatthexandycoordinatesofanypointinsideanelementcanbeexpressedasx=xe+xξ+xη12131y=ye+yξ+yη(2.70)12131wherethenotationx=xe−xe(2.71)ijijhasbeenusedin(2.70).Usingthechainruleofdifferentiation,onecanwritethat∂N∂N∂x∂N∂y=+∂ξ∂x∂ξ∂y∂ξ∂N∂N∂x∂N∂y=+(2.72)∂η∂x∂η∂y∂ηInmatrixform,⎧⎫⎡⎤⎧⎫⎪⎪∂N⎪⎪∂x∂y⎪⎪∂N⎪⎪⎨∂ξ⎬⎢∂ξ∂ξ⎥⎨∂x⎬=⎢⎥(2.73)⎪⎪∂N⎪⎪⎣∂x∂y⎦⎪⎪∂N⎪⎪⎩⎭⎩⎭∂η∂η∂η∂y P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2868INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwherethe2×2squarematrixiscalledtheJacobianmatrix,denotedbyJ,andcanbeevaluatedusingtheexpressionsin(2.70).Specifically,theJacobianmatrixisgivenbyx21y21J=(2.74)x31y31Thecoordinatetransformationin(2.73)canberearrangedbyinvertingtheJacobianmatrixandexpressingthematrixsysteminthefollowingform:⎧⎫⎧⎫⎪⎪∂N⎪⎪⎪⎪∂N⎪⎪⎨⎬⎨⎬∂x−1∂ξ=J(2.75)⎪⎪∂N⎪⎪⎪⎪∂N⎪⎪⎩⎭⎩⎭∂y∂ηwhereJ−1,whichdenotestheinverseoftheJacobianmatrix,isgivenby−11y31−y21J=(2.76)|J|−x31x21Notethat|J|isthedeterminantoftheJacobianmatrixandisgivenby|J|=x21y31−x31y21=2Ae(2.77)whereAedenotestheareaofthetriangle.ThedeterminantoftheJacobianmatrixisequaltotwicetheareaofthetriangularelementprovidedthatthelocalnodenumbersofthetrianglefollowacounter-clockwisesenseofnumbering.Thus,informingtheconnectivityinformationarrayofthefiniteelementmesh,itisinstructivethatthelocalnodesofeachtrianglebenumberedinacounter-clockwisedirection.Using(2.75)–(2.77)inconjunctionwith(2.69),itfollowsthat⎧⎫⎧⎫⎪⎪∂N1⎪⎪⎪⎪∂N1⎪⎪⎨⎬⎨⎬∂x1y31−y21∂ξ=⎪⎪∂N1⎪⎪2Ae−x31x21⎪⎪∂N1⎪⎪⎩⎭⎩⎭∂y∂η!1y31−y21−1=2Ae−x31x21−1!1y21−y31=2Aex31−x21!1y23=(2.78)2Aex32 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS69Inotherwords,∂N1y23=∂x2Ae∂N1x32=(2.79)∂y2AeSimilarly,itcanbeshownthat∂N2y31=∂x2Ae∂N2x13=(2.80)∂y2Aeand∂N3y12=∂x2Ae∂N3x21=(2.81)∂y2AeToevaluatethedoubleintegralin(2.68),itisnecessarytochangethevariablesofin-tegrationfromxandytoξandη.Inotherwords,insteadofintegratingoverthetriangularelementontheregularcoordinatesystem,itismoreconvenientthattheintegrationbecarriedoutonthemastertrianglewhichliesonthenaturalcoordinatesystem.Thetransformationofadoubleintegralfromtheregularcoordinatesystemtothenaturalcoordinatesystemisgivenby[6]11−ηf(x,y)dxdy=f(x(ξ,η),y(ξ,η))|J|dξdη(2.82)e00whichisattributedtotheGermanmathematicianCarlGustavJacobJacobi(1804–1851).Using(2.79)–(2.81)andtheJacobitransformationin(2.82),theentriesofelementmatrixMecanbeevaluatedinastraightforwardmanner.Specifically,11−ηyyxxe23233232M11=−αx+αy2Aedξdη002Ae2Ae2Ae2Ae22(y23)(x32)=−αx+αy(2.83)4Ae4AeSimilarly,eey23y31x32x13M12=M21=−αx+αy(2.84)4Ae4Ae P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2870INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSeey23y12x32x21M13=M31=−αx+αy(2.85)4Ae4Ae22e(y31)(x13)M22=−αx+αy(2.86)4Ae4Aeeey31y12x13x21M23=M32=−αx+αy(2.87)4Ae4Ae22e(y12)(x21)M33=−αx+αy(2.88)4Ae4AeNoticethatthematrixissymmetric,i.e.,Me=Me(2.89)ijjiThus,someoftheentriesdonothavetobeexplicitlyevaluated.AnotherelementmatrixthatispartofthegoverninglinearsystemofequationsismatrixTegivenby(2.52),i.e.,Te=βNNdxdy(2.90)ijijewhereβisaconstant.UsingtheJacobitransformationin(2.82),theintegralin(2.90)canbeconvenientlyexpressedoverthemastertriangularelementas11−ηTe=βNN|J|dξdηijij0011−η=β2AeNiNjdξdη(2.91)00Specifically,thefirstdiagonalentryofmatrixTeisgivenby11−ηTe=β2A2dξdη11e(N1)0011−η2=β2Ae(1−ξ−η)dξdη00βAe=(2.92)6 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS71TheremainingentriesofmatrixTecanbeevaluatedinasimilarway.However,inordertosavetimeintheevaluationofmatrixTe,thereisasimplegenericformulathatcanbeusedinstead[7]:mn!m!n!(N1)(N2)(N3)dxdy=2Ae(2.93)e(+m+n+2)!Tocheckthevalidityoftheformulain(2.93),letususeittoreevaluateentryTe.Thus,11Te=β(N2dxdy111)e2!0!0!=β2Ae(2.94)(2+0+0+2)!βAe=6whichisthesameresultobtainedin(2.92).TheremainingentriesofmatrixTearegivenbyeeβAeT12=T21=(2.95)12eeβAeT13=T31=(2.96)12eβAeT22=(2.97)6eeβAeT23=T32=(2.98)12eβAeT33=(2.99)6Next,wemustevaluatetheright-hand-sidevectorfewhoseentriesaregivenby(2.53),i.e.,fe=Ngdxdy(2.100)iieThisintegralcanbeexpressedinamoreconvenientformusingtheJacobitransformationin(2.82):11−ηfe=2AN(ξ,η)gdξdη(2.101)iei00Ifgintheargumentoftheintegralisconstant,thenitcanbetakenoutoftheintegralandproceedwithintegratingonlytheshapefunctionoverthemasterelement.If,however,gisafunctionofthespacecoordinatesxandy,thenithastobemappedtothenaturalcoordinatesystemusing(2.70)beforeproceedingwiththeintegrationoverthemasterelement.Incase P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2872INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICStheintegrationisdifficulttobeevaluatedanalytically,itmaybemoreconvenienttoevaluateitnumerically.Here,itisassumedthatgisconstantinsidetheintegraland,thus,theentriesoftheright-hand-sidevectorfearegiven,accordingto(2.93),bye1!0!0!gAef1=gN1dxdy=g2Ae=e(1+0+0+2)!3e0!1!0!gAef2=gN2dxdy=g2Ae=(2.102)e(0+1+0+2)!3e0!0!1!gAef3=gN3dxdy=g2Ae=e(0+0+1+2)!3Atlast,wemustevaluatetheright-hand-sidevectorpewhoseentriesaregivenby(2.54).Notethat(2.54)reducesto(2.67)afterimposingthemixedboundaryconditionin(2.66).Thus,theithentryoftheelementvectorpeisgivenbype=−Nq−γu)d(2.103)ii(2whereqandγareconstants.Theintegralin(2.103)isnonzeroonlyforboundaryelementsthathaveatleastoneedgecoincidingwith2andwhereeitherqorγisnonzero;ifbothqandγarezero,thenautomaticallytheintegralin(2.103)becomeszero.Forexample,toimposeaNeumannboundaryconditionon2,qandγmustbesettozero.Thisisequivalenttohavingaright-hand-sidevectorpeequaltothezerovector0.Now,foragenericmixedboundarycondition,asgivenby(2.66),withnonzeroqandγ,theintegralin(2.103)becomes3pe=−Nq−γueNdiijjLebj=1(2.104)=−Nqd+NγueN+ueN+ueNdii112233LeLebbwhereLedenotestheboundaryedgeoftheelementthatcoincideswith.b2ConsiderthetriangularelementillustratedinFigure2.8(a)withoneedgelyingon2.Numberingthelocalnodesoftheelementasindicatedinthefigure,theintegralin(2.104)mustbeevaluatedalongtheedgefromnode1tonode2(1→2).Inotherwords,pe=−Nqd+NγueN+ueN+ueNd(2.105)iii1122331→21→2Toevaluatetheintegralin(2.105),triangleemustbemappedontothemastertriangleinthenaturalcoordinatesystem,asshowninFigure2.8(b).Thus,integratingalongedge1→2on P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS73FIGURE2.8:(a)Atriangularelementwithanedgeonboundary2.(b)Mastertriangularelement.Thepathofintegrationisindicatedwitharrowstheregulartriangleisequivalenttointegratingfrom0to1alongtheξ-axisofthenaturalcoordinatesystemmultipliedbythelengthoftheedge.Inotherwords,d=12dξ(2.106)Theexpressionin(2.106)canbeprovedveryeasily.Atanypointalongedge1→2,thexandycoordinatesaregivenbyx=xeN(ξ,0)+xeN(ξ,0)1122=xe(1−ξ)+xeξ12=xe+xξ121(2.107)y=yeN(ξ,0)+yeN(ξ,0)1122=ye(1−ξ)+yeξ12=ye+yξ121Thedifferentialdisdefinedas"22d=(dx)+(dy)(2.108)wheredx=x21dξ(2.109)dy=y21dξSubstituting(2.109)into(2.108)yields"22d=(x21)+(y21)dξ(2.110)=12dξ P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2874INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSThus,usingtheresultof(2.106),theintegralin(2.105),evaluatedforlocalnode1oftheelement,becomes1pe=−N(ξ,0)qdξ111201#$+N(ξ,0)γueN(ξ,0)+ueN(ξ,0)+ueN(ξ,0)dξ1112233120(2.111)11#$=−ueeξ+ue0dξ(1−ξ)q12dξ+(1−ξ)γ1(1−ξ)+u231200q12γ12eγ12ee=−+u1+u2+0u3236Thetermsenclosedinthesquarebracketsmustbetransferredtotheleft-handsideoftheelementmatrixsystemin(2.55).Thisisequivalenttosubtractingthecoefficientsofue,ue,and12uefromthematrixentriesKe,Ke,andKe,respectively,i.e.,3111213eeγ12K11=K11−3eeγ12K12=K12−(2.112)6Ke=Ke−01313Similarly,thevectorentrythatcorrespondstolocalnode2isgiven,accordingtotheintegralin(2.105),by1pe=−N(ξ,0)qdξ221201#$+N(ξ,0)γueN(ξ,0)+ueN(ξ,0)+ueN(ξ,0)dξ2112233120(2.113)11#$=−ξqdξ+ξγueeξ+ue0dξ121(1−ξ)+u231200q12γ12eγ12ee=−+u1+u2+0u3263Again,thetermswithinthesquarebracketsmustbetransferredtotheleft-handsideoftheelementmatrixsystemandbesubtractedfromtheentriesofmatrixKe;i.e.,eeγ12K21=K21−6eeγ12K22=K22−(2.114)3Ke=Ke−02333 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS75Finally,thevectorentrythatcorrespondstolocalnode3is1pe=−N(ξ,0)qdξ331201#$(2.115)+N(ξ,0)γueN(ξ,0)+ueN(ξ,0)+ueN(ξ,0)dξ3112233120=0duetothefactthatNξ,0)=0alongtheξ-axis.Theelementright-hand-sidevectorpefor3(aboundaryelementlocallynumberedasshowninFigure2.8(a)isthereforegivenby⎧⎫⎪⎨1⎪⎬eq12p=−1(2.116)2⎪⎩⎪⎭0Exercise2.2.Givenalineartriangularelementonthexy-planewhosenodesarenumberedinacounter-clockwisedirection,provethatthedeterminantoftheJacobianmatrixistwicetheareaofthetriangle.Exercise2.3.UsingtheJacobitransformation,evaluatethedoubleintegralin(2.68)toshowthattheentriesofmatrixMearegivenby(2.83)–(2.88).FollowthesameprocedureaswasdoneforMe.RepeattheexerciseformatrixTetoshowthatitsentriesaregivenby(2.94)–(2.99).112.5.2BilinearQuadrilateralElementsTheevaluationofelementmatricesandvectorsusingbilinearquadrilateralelementsfollowsthesamestep-by-stepprocedureastheoneusedforlineartriangularelements.Thegoverninginterpolationfunctionsforabilinearquadrilateralelementaregivenby1N1=(1−ξ)(1−η)41N2=(1+ξ)(1−η)4(2.117)1N3=(1+ξ)(1+η)41N4=(1−ξ)(1+η)4Usinganisoparametricrepresentation,thex,yspacecoordinatesofapointinsideaquadrilateralandtheprimaryunknownquantityareexpandedintermsofthesameinterpolationfunctions,i.e.,x=xeN+xeN+xeN+xeN11223344(2.118)y=yeN+yeN+yeN+yeN11223344 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2876INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSandu=ueN+ueN+ueN+ueN(2.119)11223344wherexe,yefori=1,2,3,4arethenodecoordinatesofthequadrilateralelement,andueforiiii=1,2,3,4arethevaluesoftheprimaryunknownquantityatthefournodes.Usingthechainruleofdifferentiation,thepartialderivativesoftheinterpolationfunctionsin(2.117)withrespecttoξandηcanbeexpressedas∂Ni∂Ni∂x∂Ni∂y=+∂ξ∂x∂ξ∂y∂ξ(2.120)∂Ni∂Ni∂x∂Ni∂y=+∂η∂x∂η∂y∂ηorinamoreconvenientmatrixform⎧⎫⎡⎤⎧⎫⎪⎪∂Ni⎪⎪∂x∂y⎪⎪∂Ni⎪⎪⎨∂ξ⎬⎢∂ξ∂ξ⎥⎨∂x⎬=⎢⎥⎪⎪∂Ni⎪⎪⎣∂x∂y⎦⎪⎪∂Ni⎪⎪⎩⎭⎩⎭∂η∂η∂η∂y(2.121)⎧⎫∂Ni⎪⎨⎪⎬∂x=J⎪⎩∂Ni⎪⎭∂ywhereJ,theJacobianmatrix,isgivenbyJ11J12J=(2.122)J21J22with1#$J=−(1−η)xe+(1−η)xe+(1+η)xe−(1+η)xe11123441#$J=−(1−η)ye+(1−η)ye+(1+η)ye−(1+η)ye1212344(2.123)1#$J=−(1−ξ)xe−(1+ξ)xe+(1+ξ)xe+(1−ξ)xe21123441#$J=−(1−ξ)ye−(1+ξ)ye+(1+ξ)ye+(1−ξ)ye2212344 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS77Toexpressthepartialderivativeswithrespecttoxandyintermsofthepartialderivativeswithrespecttoξandη,theJacobianmatrixmustbeinverted,i.e.,⎧⎫⎧∂N⎫⎪⎨∂Ni⎪⎬⎪⎪i⎪⎪⎨⎬∂x−1∂ξ=J⎪⎩∂Ni⎪⎭⎪⎪∂Ni⎪⎪⎩⎭∂y∂η⎧⎫(2.124)⎪⎪∂Ni⎪⎪⎨⎬1J22−J12∂ξ=|J|−J21J11⎪⎪∂Ni⎪⎪⎩⎭∂ηwhere|J|isthedeterminantoftheJacobianmatrixgivenby|J|=J11J22−J12J21(2.125)Asanexample,letusconsidertheinterpolationfunctionthatcorrespondstonode1:∂N11=−(1−η)∂ξ4(2.126)∂N11=−(1−ξ)∂η4Substituting(2.126)into(2.124),wehave⎧⎫⎪⎨∂N1⎪⎬1!∂x1J22−J12−4(1−η)=(2.127)⎪⎩∂N1⎪⎭|J|−J21J11−1(1−ξ)4∂ywhichcanalsobewrittenas∂N11=[−J22(1−η)+J12(1−ξ)]∂x4|J|(2.128)∂N11=[J21(1−η)−J11(1−ξ)]∂y4|J|Similarly,itcanbeshownthat∂N21=[J22(1−η)+J12(1+ξ)]∂x4|J|(2.129)∂N21=[−J21(1−η)−J11(1+ξ)]∂y4|J| P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2878INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS∂N31=[J22(1+η)−J12(1+ξ)]∂x4|J|(2.130)∂N31=[−J21(1+η)+J11(1+ξ)]∂y4|J|∂N41=[−J22(1+η)−J12(1−ξ)]∂x4|J|(2.131)∂N41=[J21(1+η)+J11(1−ξ)]∂y4|J|Thefirstmatrixtobeevaluated,onanelementbasis,isMe.Theentriesofthismatrixaregivenbye∂Ni∂Nj∂Ni∂NjMij=−αx+αydxdy(2.132)e∂x∂x∂y∂ywhereeisthedomainofthequadrilateralelement.UsingtheJacobitransformationintroducedfortriangularelementsin(2.82),modifiedhoweverforquadrilateralelements,theintegralin(2.132)canbeexpressedinamoreconvenientformgivenby11∂N∂N∂N∂NeijijMij=−αx+αy|J|dξdη(2.133)−1−1∂x∂x∂y∂yAsseen,thedoubleintegralisnowevaluatedoverthedomainofthemasterelementwhichcorrespondstoasquare[seeFigure2.6(b)].NoticethatthepartialderivativesofthegoverninginterpolationfunctionswithrespecttoxandywereconvenientlyexpressedusingtheJacobianmatrixintermsofthenaturalcoordinatesξandη.Theseexpressionsaregivenby(2.128)–(2.131).Asanexample,letusevaluateentryMe:112211∂N∂Ne11M11=−αx+αy|J|dξdη−1−1∂x∂y!(2.134)11112αx[−J22(1−η)+J12(1−ξ)]+=−2dξdη16−1−1|J|αy[J21(1−η)−J11(1−ξ)]Unlikethecaseofusinglineartriangularelements,theintegrandof(2.134)isafairlycomplexfunctionofξandη,whichmakestheintegraldifficulttoevaluateanalytically.Thisintegralmustbeevaluatednumerically.Althoughtherearemanynumericalmethods[8,9]thatcanbeusedtoevaluateadoubleintegration,someofwhichincludemidpointruleandtrapezoidrule,themostappropriateandwidelyusedapproachintheFEMisGaussquadrature. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS79Tointroducethisapproach,letusfirstconsiderthesingleintegral1I=f(ξ)dξ(2.135)−1Usingann-pointGaussquadratureapproximation,theintegralin(2.135)canbeexpressedasfollows:nI≈wif(ξi)(2.136)i=1wherewi’saretheGaussweightsandξi’saretheGausspoints.Table2.1tabulatesthevaluesofGaussweightsandpointsstartingfromaone-point(i.e.,n=1)quadratureallthewaytoaseven-point(i.e.,n=7)quadrature.ForacompletetableofGaussweightsandpointsuptoandincludinga13-pointquadrature,thereaderisreferredto[10].Ann-pointGaussquadraturebecomesexactforpolynomialsofdegree(2n−1)orless.NoticethatinordertopreservetheTABLE2.1:WeightsandpointsforGaussquadraturenξiwi10.00000000002.00000000002±0.57735026921.000000000030.00000000000.8888888889±0.77459666920.55555555564±0.33998104360.6521451549±0.86113631160.347854845150.00000000000.5688888889±0.53846931010.4786286705±0.90617984590.23692688516±0.23861918610.4679139346±0.66120938650.3607615730±0.93246951420.171324492470.00000000000.4179591837±0.40584515140.3818300505±0.74153118560.2797053915±0.94910791230.1294849662 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2880INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSaccuracyofGaussquadraturerule,alargenumberofsignificantfiguresmustbeconsideredwhenstoringtheGaussweightsandpointsinacomputer.ItisalsorecommendedthatdoubleprecisionarithmeticbeusedwhenimplementingGaussquadratureinacomputercode.Gaussquadraturecanalsobeappliedtodoubleintegralsoftheform11I=f(ξ,η)dξdη(2.137)−1−1Employingann-pointGaussquadratureapproximation,onecanwritethat1nI≈wif(ξi,η)dη−1i=1nn≈wjwifξi,ηj(2.138)j=1i=1nn≈wiwjfξi,ηjj=1i=1Similarly,Gaussquadraturecanbeextendedtovolumeintegralsoveracube.Itisimportant,however,toemphasizeherethatGaussweightsandpointsfortrianglesaredifferentfromtheonestabulatedinTable2.1.FormoreinformationonGaussquadraturefortriangulardomains,referto[10,11].NowthatthereaderhasbecomefamiliarwithGaussquadrature,theintegralin(2.134)canbeevaluatednumericallyusing(2.138)andthedatapostedinTable2.1.Thesecondelementma-trixinvolvedinthe2-DnodalfiniteelementformulationismatrixTewhoseentriesaregivenbyTe=βNNdxdyijije11(2.139)=βNiNj|J|dξdη−1−1AswithmatrixMe,theintegrandinvolvedinmatrixTeisalsocomplexandthatmakesitextremelydifficulttointegrateanalytically.Thus,choosingtheappropriaten-pointGaussquadrature,dependingonthedegreeofthepolynomialinvolved,theintegralin(2.139)canbeevaluatednumerically.Asanexample,entryTebecomes11nn# $%%Te≈wwβNξ,η2%Jξ,η%(2.140)11ij1ijijj=1i=1wherenmustbesuchthatthedegreeofthepolynomialinvolvedintheintegrandis(2n−1)orless.TheremainingentriesofmatrixTeareevaluatedinasimilarway.Inaddition,the P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS81right-hand-sidevectorfe,whoseentriesaregivenbyfe=Ngdxdyiie11(2.141)=Nig|J|dξdη−1−1canbeevaluatedusingthesameapproach.Asanexample,considerentryfe:1nn%%fe≈wwNξ,ηgξ,η%Jξ,η%(2.142)1ij1ijijijj=1i=1TocompletetheevaluationofallelementmatricesinvolvedinthenodalfiniteelementformulationofthegenericBVPathandusingbilinearquadrilateralelements,wemustalsoevaluatetheright-hand-sidevectorpewhoseentriesaregiven,accordingto(2.104),by#$pe=−Nqd+NγueN+ueN+ueN+ueNd(2.143)iii11223344LeLebbwhereLedenotesaboundaryedgethatcoincideswith.Ifaquadrilateralelementhasnob2edgeonboundary2,thenpe=0fori=1,2,3,4(2.144)iIfaquadrilateralhasoneormoreedgesonboundary2,thetwointegralsinvolvedin(2.143)mustbeevaluatedalongeachoftheseboundaryedges.Toillustratehowthesetwointegralsareevaluated,considertheboundaryquadrilateralelementshowninFigure2.9(a)withoneedgelyingon2.Basedonthisfigureandtheassociatednodalnumberingscheme,theintegralin(2.143)becomes#$pe=−Nqd+NγueN+ueN+ueN+ueNd(2.145)iii112233441→21→2h143(–1,1)(1,1)2Γ2exy(–1,–1)(1,–1)12x34(a)(b)FIGURE2.9:(a)Aquadrilateralelementwithanedgeonboundary2.(b)Masterquadrilateralelement.Thepathofintegrationisindicatedwitharrows P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2882INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSEvaluationoftheintegralin(2.145)requiresthatquadrilateraleinFigure2.9(a)bemappedontothemasterelementshowninFigure2.9(b).Doingthat,theintegralin(2.145)canbeequivalentlyexpressedas11ueN(ξ,−1)+ueN(ξ,−1)+pe=−12N(ξ,−1)qdξ+12N(ξ,−1)γ1122dξiiiueN(ξ,−1)+ueN(ξ,−1)2−12−13344(2.146)Noticethatintegratingalongedge1→2isequivalenttointegratingonthenaturalcoordinatesystemalongtheξ-axisfromξ=−1toξ=1withdreplacedby12d=dξ(2.147)2Theexpressionin(2.147)canbederivedfairlyeasilyoncewerealizethatalongedge1→2thenaturalcoordinateη=−1;i.e.,x=xeN(ξ,−1)+xeN(ξ,−1)+xeN(ξ,−1)+xeN(ξ,−1)1122334411=xe(1−ξ)+xe(1+ξ)+xe0+xe0122234(2.148)xe+xexe−xe1221=+ξ22Similarly,itcanbeshownthatye+yeye−ye1221y=+ξ(2.149)22Thedifferentialdisgivenby"22d=(dx)+(dy)(2.150)wherexe−xex2121dx=dξ=dξ22ee(2.151)y2−y1y21dy=dξ=dξ22Substituting(2.151)into(2.150)yields"122d=(x21)+(y21)dξ2(2.152)12=dξ2 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS83Now,letusevaluate(2.146)fori=1,2,3,and4assumingqandγareconstants.Specifically,fori=1,11ueN(ξ,−1)+ueN(ξ,−1)+pe=−12N(ξ,−1)qdξ+12N(ξ,−1)γ1122dξ111ueN(ξ,−1)+ueN(ξ,−1)2−12−1334411#$12q12qeeee=−(1−ξ)dξ+(1−ξ)u1(1−ξ)+u2(1+ξ)+u30+u40dξ4−18−112qγ12eγ12eee=−+u1+u2+0u3+0u4(2.153)2312Aswasdoneinthecaseoftriangularelements,thetermsenclosedinthesquarebracketsmustbetransferredtotheleft-handsideoftheelementmatrixin(2.55).Thisisequivalenttosubtractingthecoefficientsofue,ue,ue,anduefromthematrixentriesKe,Ke,Ke,and1234111213Ke,respectively,i.e.,14eeγ12K11=K11−3eeγ12K12=K12−(2.154)12Ke=Ke−01313Ke=Ke−01414Similarly,e12qγ12eγ12eeep2=−+u1+u2+0u3+0u4(2.155)2123Again,thetermsinsidethesquarebracketsmustbetransferredtotheleft-handsideoftheelementmatrixsystemandbesubtractedfromtheentriesofmatrixKe;i.e.,eeγ12K21=K21−12eeγ12K22=K22−(2.156)3Ke=Ke−02323Ke=Ke−02424TheothertwoentriesofvectorpearezerosinceN3(ξ,−1)=0(2.157)N4(ξ,−1)=0 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2884INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSThus,forthequadrilateralelementonboundary2depictedinFigure2.9(a),theright-hand-sidevectorpeisgivenby⎧⎫⎪⎪1⎪⎪⎪⎨⎪⎬e12q1p=−(2.158)2⎪⎪0⎪⎪⎪⎩⎪⎭0Itisimportanttopointoutherethatforeveryelementwithoneedgeonboundary2,thecorrespondingentriesofmatrixKemustbeupdatedaccordingto(2.154)and(2.156).ThisisequivalenttoupdatingtheentriesofmatrixKeonceforeveryedgeonboundary.Inother2words,foredge1→2,therelevantentriesofmatrixKemustbeupdatedaccordingtoeeγ12K=K−11113eeγ12K=K−121212(2.159)eeγ12K=K−212112eeγ12K=K−22223whichresultsaftercombining(2.154)and(2.156).Iftheedgeonboundary2isdefinedbylocalnodenumbers2and3,thematrixupdatemustoccuronentriesKe,Ke,Ke,andKe22233233insteadofKe,Ke,Ke,andKe.Furthermore,theright-hand-sidevectorpewouldbegiven11122122by⎧⎫⎪⎪0⎪⎪⎪⎨⎪⎬e23q1p=−(2.160)2⎪⎪1⎪⎪⎪⎩⎪⎭0Now,iftheelementhastwoedgesonboundary2,asshowninFigure2.10,thecorre-spondingvectorpewilltaketheform⎧⎫⎧⎫⎪⎪1⎪⎪⎪⎪0⎪⎪⎪⎨⎪⎬⎪⎨⎪⎬e12q123q1p=−−(2.161)2⎪⎪0⎪⎪2⎪⎪1⎪⎪⎪⎩⎪⎭⎪⎩⎪⎭00 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS8532e14FIGURE2.10:Aquadrilateralelementwithtwoedgesonboundary2.ThepathofintegrationisindicatedwitharrowswhereastheupdateofthecorrespondingentriesofmatrixKebecomeseeγ12K=K−11113eeγ12K=K−121212eeγ12K=K−212112eeγ12γ23(2.162)K=K−−222233eeγ23K=K−232312eeγ23K=K−323212eeγ23K=K−33333Exercise2.4.ShowthattheJacobianmatrixforabilinearquadrilateralelementisgivenby(2.122)–(2.123).Exercise2.5.Usingone-point,two-point,andthree-pointGaussquadrature,evaluatethelineintegral13x2+e2xdx−1Usetheexactanswertocalculatethecorrespondingerrorforeachcase. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2886INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSExercise2.6.Usingone-pointandtwo-pointGaussquadrature,evaluatethedoubleintegral112√2xy−xydxdy−1−1Usetheexactanswertocalculatethecorrespondingerrorforeachcase.2.6ASSEMBLYOFTHEGLOBALMATRIXSYSTEMIntheprevioussection,thegoverningelementmatricesandvectorsfortriangularandquadri-lateralelementsweregenerated.Foreachelementinthedomain,thereexistacoefficientmatrixKeandaright-hand-sidevectorbethatmustbemapped,accordingtotheconnectivityinfor-mation,andaddedtotheglobalcoefficientmatrixKandtheglobalright-hand-sidevectorb,respectively.Theprocessofmappingandaddingtheentriesofeachelementcoefficientmatrixandright-hand-sidevectortotheentriesoftheglobalcoefficientmatrixandright-hand-sidevector,respectively,iscalledassemblyprocess.Thedimensionoftheelementcoefficientmatrixisequaltothenumberofnodesoftheelement.Forexample,alineartrianglecorrespondstoa3×3elementcoefficientmatrix,whereasabilinearquadrilateralcorrespondstoa4×4elementcoefficientmatrix.Thedimensionoftheglobalcoefficientmatrixisequaltothetotalnumberofnodesinthefiniteelementdomain.Asanexample,considerthe2-DdomainillustratedinFigure2.11,whichisdiscretizedusingquadrilateralelements.Thetotalnumberofnodesinthefiniteelementdomainiseightand,therefore,thesizeoftheglobalcoefficientmatrixis8×8.Thenodesareuniquely32412518367FIGURE2.11:Discretizationofa2-Ddomainusingthreequadrilateralelements P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS87TABLE2.2:ElementconnectivityinformationLOCALNODENUMBERSELEMENTNUMBER,en1n2n3n4118322854337658numberedfrom1to8,asshowninFigure2.11.Thesearereferredtoastheglobalnodenumbersasopposedtothelocalnodenumbers(1–4forquadrilaterals)whicharereferredtotheelementitself.Usingasimilarargument,theelementright-hand-sidevectoris4×1(forbilinearquadrilateralelements)andtheglobalright-hand-sidevectoris8×1(forthefiniteelementdomainillustratedinFigure2.11).Theassemblyprocessbeginswiththeformationoftheelementconnectivityarray,whichholdstheglobalnodenumbersofeachelementinacounter-clockwisedirection,asshowninTable2.2.Notethatitdoesnotmatterwhichparticularnodegoesfirstorlast.Inacomputercode,thisinformationisstoredinsidea2-Darrayn(e,i)=nifori=1,2,3,4(2.163)whereedenotestheelementnumber,andnidenotestheglobalnodenumberthatcorrespondstotheithlocalnodenumber.Forexample,theelementconnectivityarrayforelement1inFigure2.11isgivenbyn(1,1)=1n(1,2)=8(2.164)n(1,3)=3n(1,4)=2Similarly,forelement2wehaven(2,1)=8n(2,2)=5(2.165)n(2,3)=4n(2,4)=3andsoon. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2888INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSOncetheelementconnectivityarrayisformed,an8×8globalmatrixiscreatedwithitsentriesinitializedtozero.Thus,initiallytheglobalcoefficientmatrixisanullmatrix⎡⎤00000000⎢⎥⎢00000000⎥⎢⎥⎢00000000⎥⎢⎥⎢00000000⎥K=⎢⎥(2.166)⎢00000000⎥⎢⎥⎢⎥⎢00000000⎥⎢⎥⎣00000000⎦00000000Theprocessofassemblybeginsbyloopingthroughalltheelementsone-by-oneandupdatingtheentriesoftheglobalcoefficientmatrixaccordingtothefollowingMatlabalgorithm:fore=1:3%loopthroughtheelementsinthedomainfori=1:4%loopthroughthelocalnodes(1stloop)ofelementeforj=1:4%loopthroughthelocalnodes(2ndloop)ofelementeK(n(e,i),n(e,j))=K(n(e,i),n(e,j))+ke(i,j);endendendwhereKdenotestheglobalcoefficientmatrixandkedenotestheelementcoefficientmatrix.Bycompletingtheassemblyofelement1,theglobalcoefficientmatrixhasthefollowingform:⎡(1)(1)(1)(1)⎤KKK0000K11141312⎢(1)(1)(1)(1)⎥⎢KKK0000K⎥⎢41444342⎥⎢(1)(1)(1)(1)⎥⎢⎢K31K34K330000K32⎥⎥⎢⎥⎢00000000⎥K=⎢⎥(2.167)⎢00000000⎥⎢⎥⎢⎥⎢00000000⎥⎢⎥⎢⎥⎣00000000⎦(1)(1)(1)(1)KKK0000K21242322Inotherwords,accordingtotheelementconnectivityinformationinTable2.2,elemententry(1)(1)K23mapstoglobalentryK83,elemententryK42mapstoglobalentryK28,andsoon.Itisalsoworthnotingherethatthematrixin(2.167)issymmetricaboutthemaindiagonal.Thiscanbededuceddirectlyfromtheintegralsdescribingeachentryintheweakformulationoftheproblem. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS89Oncethesecondquadrilateralelementisassembled,theglobalcoefficientmatrixisaug-mentedbytheentriesofthecorrespondingelementmatrix.Specifically,theglobalcoefficientmatrixbecomes⎡⎤(1)(1)(1)(1)KKK0000K11141312⎢⎥⎢K(1)K(1)K(1)0000K(1)⎥⎢41444342⎥⎢⎥⎢K(1)K(1)K(1)+K(2)K(2)K(2)00K(1)+K(2)⎥⎢3134334443423241⎥⎢(2)(2)(2)(2)⎥⎢00KKK00K⎥K=⎢34333231⎥(2.168)⎢(2)(2)(2)(2)⎥⎢00KKK00K⎥⎢24232221⎥⎢⎥⎢00000000⎥⎢⎥⎢00000000⎥⎣⎦(1)(1)(1)(2)(2)(2)(1)(2)KKK+KKK00K+K2124231413122211TocompletetheassemblyprocessforthefiniteelementmeshinFigure2.11,wemustalsoassembleelement3.Followingthesamealgorithmdescribedabove,thefinalformoftheglobalcoefficientmatrix,afterallthreequadrilateralelementshavebeenassembled,becomes⎡⎤(1)(1)(1)(1)KKK0000K11141312⎢⎥⎢K(1)K(1)K(1)0000K(1)⎥⎢41444342⎥⎢⎥⎢K(1)K(1)K(1)+K(2)K(2)K(2)00K(1)+K(2)⎥⎢3134334443423241⎥⎢(2)(2)(2)(2)⎥⎢00KKK00K⎥K=⎢34333231⎥⎢(2)(2)(2)(3)(3)(3)(2)(3)⎥⎢00K24K23K22+K33K32K31K21+K34⎥⎢⎥⎢0000K(3)K(3)K(3)K(3)⎥⎢23222124⎥⎢⎥⎢0000K(3)K(3)K(3)K(3)⎥⎣13121114⎦(1)(1)(1)(2)(2)(2)(3)(3)(3)(1)(2)(3)KKK+KKK+KKKK+K+K212423141312434241221144(2.169)Toassembletheglobalright-hand-sidevector,theprocessisverysimilar.Wealwaysstartwithaninitialized-to-zero8×1globalvectorbwhoseentriesareupdatedaccordingtothefollowingMatlabalgorithm:fore=1:3%loopthroughtheelementsinthedomainfori=1:4%loopthroughthelocalnodesofelementeb(n(e,i))=b(n(e,i))+be(i);endend P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2890INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSwherebistheglobalright-hand-sidevectorandbeistheelementright-hand-sidevector.Basedontheabovealgorithm,oncetheassemblyprocessisover,theglobalright-hand-sidevectorwillhavethefollowingform:⎧⎫⎪⎪(1)⎪⎪b⎪⎪1⎪⎪⎪⎪(1)⎪⎪⎪⎪b4⎪⎪⎪⎪⎪⎪⎪⎪(1)(2)⎪⎪⎪⎪b3+b4⎪⎪⎪⎪⎪⎪⎨b(2)⎬3b=(2.170)⎪⎪b(2)+b(3)⎪⎪⎪⎪23⎪⎪⎪⎪(3)⎪⎪⎪⎪b⎪⎪⎪⎪2⎪⎪⎪⎪(3)⎪⎪⎪⎪b⎪⎪⎪⎪1⎪⎪⎩(1)(2)(3)⎭b+b+b2142.7IMPOSITIONOFBOUNDARYCONDITIONSTherearetwotypesofboundaryconditionsthatmustbeimposedfortheBVPathand:theDirichlettypeandthemixedtype.RememberthattheNeumannboundaryconditionisaspecialcaseofthemixedboundarycondition.TheprocedurethathastobefollowedinordertoproperlyimposetheDirichletboundaryconditiononalltheglobalnodesthatlieonboundary1isexactlythesameastheprocedureoutlinedinSection1.8forthe1-Dcase.Asfarasthemixedboundaryconditionisconcerned,thishasalreadybeenincorporatedintheweakformulation,whichwasdevelopedinSection2.4andevaluatedforspecificinterpolationfunctionsinSection2.5.Re-memberthattheexistenceofamixedboundaryconditionon2resultsinanadditionallinein-tegral,givenby(2.67),whichmustbeevaluatedalongalledgesthatcoincidewiththisboundary.2.8SOLUTIONOFTHEGLOBALMATRIXSYSTEMTheFEMwhenappliedtoaBVPalwaysresultsinasetoflinearequationsthatisusuallypresentedinamatrixform:Ku=b(2.171)whereKistheglobalcoefficientmatrix,uistheglobalvectorrepresentingtheprimaryunknownquantityatthenodesofthedomain,andbistheglobalright-hand-sidevector.Thesizeoftheglobalcoefficientmatrixisequaltothetotalnumberofnodesinthefiniteelementdomain.Animportantcharacteristicoftheglobalcoefficientmatrixisitssparsitymeaningthatmostofitsentriesarezero.Astraightforwardapproachtosolvethelinearsystemin(2.171)isbyusingLUdecomposition,wherematrixKisdecomposedintoalowertriangularmatrixLmultipliedwithanuppertriangularmatrixU.However,thisapproachwillcertainlydestroythesparsityof P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS91thematrixandpopulatethememoryofyourcomputer.Forverylargeproblems,inotherwordsproblemswithaverylargenumberofnodesandthereforeunknowns,itwillbealmostimpossibletosolvesuchalinearsystemusingLUdecompositionbecauseoflimitationsincomputermemoryandanextensivecode-executiontime.Suchanapproachcanonlybeimplementedforrelativelysmallfiniteelementproblems,i.e.,problemswithasmallnumberofnodes.UsingMatlab,onecanmakeuseofthecommandsparsetotakeadvantageofthesparsityoftheglobalcoefficientmatrixwithoutallocatingmemoryforthezeroentries.Inprofessionalcommercialfiniteelementpackages,theresultinglinearsystemofequationsissolvedusingiterativetechniques.Thistypeoftechniquesdoesnotdestroythesparsityofthematrix,andcomputermemorydoesnotgrowoutofproportions.Suchtechniquesincludeconjugategradient(CG)andbiconjugategradient(BiCG)methodscombinedwithavarietyofpreconditionersandaccelerators.Otherpopulariterativemethodsthatareusedinthesolutionoftheglobalfiniteelementmatrixsystemincludethegeneralizedminimalresidual(GMRES)methodandthequasiminimalresidual(QMR)method.WewillnotdiscussthesemethodsheresincetheycanbefoundinavarietyoflinearalgebrabooksundertheumbrellaofIterativeMethodsinLinearAlgebra[12–15].2.9POSTPROCESSINGOFTHERESULTSBysolvingthematrixsystemin(2.171),avectorcontainingthevaluesoftheprimaryunknownquantityattheglobalnodesofthefiniteelementdomainisobtained.Toevaluatetheprimaryunknownquantityatapointinsideanelement,wehavetomakeuseofthegoverningsetofinterpolationfunctions,i.e.,nu(x,y)=ueN(x,y)(2.172)iii=1NoticethattheinterpolationfunctionsNimustbeexpressedintermsofxandy.InSection2.3,theinterpolationfunctionsforbothtriangularandquadrilateralelementsweregivenintermsofthenaturalcoordinates.TobeabletoexpressNiintermsofxandy,itisimportantthatthenaturalcoordinates,ξandη,beexpressedintermsofxandy,first.Thiswillbedoneherefortriangularelements,butthesameexactprocedurecanbeappliedforquadrilateralelementsaswell.From(2.22),thexandycoordinatesofapointinsideatrianglearegivenbyx=xe+xξ+xη12131(2.173)y=ye+yξ+yη12131whichcanalsobewrittenas!!x−xexxξ1=2131(2.174)y−yeyyη12131 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2892INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSTosolveforξandη,weneedtoinvertthe2×2squarematrixgiving!!ξ1y−xx−xe=31311(2.175)ηxy−xy−yxy−ye2131312121211SubstitutingξandηintheinterpolationfunctionsNi(ξ,η),thelatterbecomefunctionsofxandy;i.e.,Ni(x,y).2.10APPLICATIONPROBLEMSInthissection,thenodalFEMwillbeusedtosolve2-Delectromagneticproblems.Specif-ically,theunderlinedmethodwillbeappliedseparatelytosolveanelectrostaticproblemandascatteringproblem.TheelectrostaticproblemisdescribedbytheLaplace’sequationandasetofDirichletboundaryconditionsimposedontheoutersurfaceofarectangulardomain.Thescatteringproblemisatime-harmonicproblemdescribedbythehomogeneousscalarwaveequation(Helmholtzequation)inconjunctionwithanABCthatisusedtoeffectivelytruncatetheunboundeddomain.Theeffectivenessandaccuracyofthemethodforthesolutionofstaticaswellastime-harmonicproblemswillbeevaluatedbycomparingthenumericalsolutiontotheexactanalyticalsolution.Bothproblemsconsideredinthissectionhaveanexactanalyticalsolution.ItisextremelyimportanttoemphasizeherethatthecapabilityoftheFEMextendsbeyondthescopeofthesetworepresentativeelectromagneticproblems.Duringthelastfewdecades,theFEMhasbeensuccessfullyappliedtoaplethoraofelectromagneticproblemsincludingscatteringby2-Dand3-Dstructures[16],antennaradiationandwavepropagation[17],anisotropicandfrequency-dependentmaterial(e.g.,plasma,ferrites,etc.)[18],phasedarrays[19],eigenvalueproblems[20],andmanymore.SpecialemphasiswasalsoplacedonhybridizingtheFEMwithboundaryintegral(BI)formulations[21]andthegeometricaltheoryofdiffraction(GTD)[22].2.10.1ElectrostaticBoundary-ValueProblemProblemdefinition:ConsidertheinfinitelylongrectangularboxwithmetallicwallsshowninFigure2.12.Theverticalandbottomwallsaremaintainedatazeroelectricpotentialwhereasthetopwall,whichisseparatedbytinygapsfromtheverticalsidewalls,hasafixedelectricpotentialofV0.Theregioninsidetheboxisfreeofcharge.UsetheFEMtosolvetheLaplace’sequationsubjecttothegivensetofboundaryconditionsinordertofindandplottheelectricpotentialdistributionintheinteriorofthebox. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS93yVV=0hV=0V=00x0wV=0FIGURE2.12:InfinitelylongrectangularboxwithmetallicwallsSolution:TheLaplace’sequation∂2V∂2V+=0(2.176)∂x2∂y2subjecttoasetofDirichletboundaryconditionsV(x,0)=0V(x,h)=V0(2.177)V(0,y)=0V(w,y)=0canbesolvedanalyticallyusingthemethodofseparationofvariables.Thecloseformsolutionoftheelectricpotentialasafunctionofxandyintheinteriorofthemetallicboxisgivenby[23]4V∞sin(2k−1)πxsinh(2k−1)πy0wwV(x,y)=(2.178)π(2k−1)sinh(2k−1)πhk=1wAcontourplotoftheelectricpotentialdistributionbasedonthecloseformexpressionin(2.178)isdepictedinFigure2.13.Thedimensionsoftherectangularboxarespecifiedtobe1×1mandtheconstantpotentialimposedatthetopsurfaceoftheboxisequalto1V.ThisBVPwassolvedusingtheFEMbasedonthe2-Dnodalanalysisoutlinedinthischapterfortriangularandquadrilateralelements.OneofthemainadvantagesoftheFEMoverothernumericaltechniquesisitsabilitytotreatdomainsofarbitraryshapeandnotnecessarilyofacanonicalshape,asisthedomainoftheproblemathand.Thisspecificdomainwaschosenbecausethereisanexactanalyticalsolutionforitwithwhichwecancompareinordertoevaluate P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2894INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFIGURE2.13:Contourplotoftheelectricpotentialdistributioninsidethebox(analyticalsolution)theaccuracyofthenumericalmethod.Tobemorespecific,thedomainoftheproblemathandisdiscretizedusinglineartriangularelements,althoughquadrilateralelementscouldalsobeused.AcoarsemeshofthemetallicsquareboxisshowninFigure2.14.Eachelementofthemeshisgivenauniquenumber(inourcasefrom1to32)knownastheelementnumber.Eachelementhasthreelocalnodeswhicharenumberedinacounter-clockwisedirectionfrom1to212223242525272931262830321617181920171921231820222411121314159111315101214166789101357246812345FIGURE2.14:Finiteelementmeshusingtriangularelements.Elementnumbersareshowninbold P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS953.Inaddition,eachnodeisgivenauniquenumberthatiswidelyreferredtoastheglobalnodenumber.FortheparticularmeshshowninFigure2.14,thereare25globalnodeswhicharenumberedsequentiallyfrom1to25.Thus,theglobalcoefficientmatrixwillhavedimensions25×25,whereastheglobalright-hand-sidevectorwillhavedimensions25×1.Theglobalmatrixandright-hand-sidevectorareformedthroughtheprocessofassembly,asoutlinedinSection2.6.Now,toconstructtheelementmatricesandvectors,whichareusedduringtheassemblyoftheglobalmatrixandvector,itisimportanttorealizethattheLaplace’sequationin(2.176)isaspecialcaseofthegenericsecond-orderpartialdifferentialequationin(2.1),i.e.,∂∂u∂∂uαx+αy+βu=g(2.179)∂x∂x∂y∂yBycomparing(2.179)with(2.176),itcanbededucedthatu=Vαx=αy=1(2.180)β=0g=0Inaddition,forthespecificproblemathand,thereexistonlyDirichletboundarycondi-tionsandnomixedboundaryconditions.Inotherwords,thereisnoboundary2.ConcerningtheDirichletboundaryconditions,notethatthreeofthemarehomogeneouswhereasoneofthemisnonhomogeneous.A2-DnodalfiniteelementMatlabcode,namedFEM2DLBox,waswrittentosolvetheLaplace’sequationwiththeassociatedDirichletboundaryconditions.Thecodeutilizeslineartriangularelements.Contourplotsoftheelectricpotentialintheinte-riorofthemetallicboxaredepictedinFigures2.15(a)and2.15(b)foracoarseandafinemesh,respectively.ByvisuallycomparingthecontourplotsinFigures2.15(a)and2.15(b)withthecorrespondingplotinFigure2.13,itisobservedthatagoodrepresentationoftheexactsolutioncanbeobtainedusingtheFEMprovidedthediscretizationofthedomainisfineenough.TheaccuracyofthenumericalsolutioncanbeevaluatedbycomputingtheerrorasafunctionofmeshsizeusingtheL2-normdefinition&''Np(1en2error=u−u(2.181)iiNpi=1whereuerepresentstheexactanalyticalsolutionataspecificgridpoint,unrepresentstheiinumericalsolutionatthesamepoint,andNprepresentsthetotalnumberofpointswherethetwosolutionsareevaluated.AplotofthenumericalerrorasafunctionofdiscretizationsizeisillustratedinFigure2.16.Thehorizontalaxiscorrespondstothelengthofanedgealongthex P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28(a)(b)FIGURE2.15:(a)Finiteelementsolutionusingacoarsemesh.(b)Finiteelementsolutionusingafinemesh96 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS970.070.060.05norm)2L0.040.030.02Numericalerror(0.010.00.020.040.060.080.10.120.140.160.180.2hFIGURE2.16:NumericalerrorbasedontheL2normasafunctionofdiscretizationsizeorydirection.Notethatthetotalwidthandheightoftheboxareequaltounity.FromFigure2.16,itisobservedthatbydecreasingthesizeofelements,whichimpliesincreasingthenumberofelementsinthemesh,thenumericalerroraccordingto(2.181)ismonotonicallyreduced.2.10.2Two-DimensionalScatteringProblemProblemdefinition:ATMzuniformplanewaveisnormallyincidentuponaperfectlyconduct-ingcircularcylinderofradiusa,asshowninFigure2.17.TheincidentelectricfieldcanbewrittenasEinc=aˆEe−jk0x=aˆEe−jk0ρcosφ(2.182)z0z0whereE0istheamplitudeoftheplanewave,k0isthepropagationconstantinfreespace,ρistheradialdistancefromthecenterofthecylindertotheobservationpoint,andφisthecorrespondinganglemeasuredfromthepositivex-axis.Assumingthata=0.5λandE0=1V/m,calculatethetotalelectricfield(incidentfield+scatteredfield)asafunctionofangleφataradialdistanceρ=λawayfromthecenterofthecylinder.Notethattheexactanalyticalsolutiontotheproblemisgivenby[24]+∞Jka)E=Ej−nJkρ)−n(0H(2)ejnφ(2.183)z0n(0(2)n(k0ρ)n=−∞Hn(k0a)whereEzcorrespondstothetotalelectricfield. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:2898INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSyriHafxrizEσ=∞FIGURE2.17:ATMzplanewaveisincidentuponaperfectlyconductingcircularcylinderSolution:Thisscatteringproblem,likeanyotherelectromagneticproblem,isgovernedbythewell-knowntime-harmonicMaxwellequationsgivenby∇×E=−jωμ0μrH(2.184)∇×H=J+jωε0εrE(2.185)ρv∇·E=(2.186)ε0εr∇·H=0(2.187)However,theproblemathandisa2-Dscatteringproblemwithanincidentelectricfieldpolarizedalongthez-direction.Asaresult,thescatteredelectricfieldfromtheinfinitelylong(alongthez-direction)circularcylinderwillalsobepolarizedalongthez-direction.Inotherwords,thetotalelectricfieldwillhaveonlyaz-directedcomponent.Inaddition,therewillbenofieldvariationsalongthez-direction,i.e.,∂Ez=0(2.188)∂zForthiscase,thetime-harmonicMaxwellequationsin(2.184)–(2.187)canbesimplifiedcon-siderably.Theprocessofsimplificationbeginsbywriting(2.184)as1H=−∇×E(2.189)jωμ0μrandsubstitutingthelatterinto(2.185)toyield11−∇×∇×E=J+jωε0εrE(2.190)jωμ0μr P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS99Therelativepermeabilityμrcannotbetakenoutofthecurl–curlexpressionunlessitisaconstant.Toretaingenerality,itwasdecidedthattherelativepermeabilityofthemediumbetreatedasaspacedependentquantity.Now,sincethetotalelectricfield(incidentplusscattered)isz-directed,wecanwritethat1∇×∇×(aˆE=−jωμaˆJ+k2εaˆE(2.191)zz)0(zz)0r(zz)μrwherek2=ω2με(2.192)000UsingthedefinitionofthecurlofavectorinCartesiancoordinates∂Ez∂Ey∂Ex∂Ez∂Ey∂Ex∇×E=aˆx−+aˆy−+aˆz−(2.193)∂y∂z∂z∂x∂x∂yandthefactthatfieldvariationalongthez-directioniszero,see(2.188),thecurl–curlexpressionin(2.191)canbewrittenas1∂1∂Ez∂1∂Ez∇×∇×(aˆzEz)=aˆz−−(2.194)μr∂xμr∂x∂yμr∂ySubstituting(2.194)into(2.191),eliminatingtheunitvectoralongthez-direction,andrear-rangingtheterms,wehave∂1∂Ez∂1∂Ez2++k0εrEz=jωμ0Jz(2.195)∂xμr∂x∂yμr∂ywhichisknownastheinhomogeneousscalarwaveequation.Incasewehaveasource-freeregion,asistheproblemathand,theimpressedsourcecurrentJzmustbesettozero,thereforehaving∂1∂Ez∂1∂Ez2++k0εrEz=0(2.196)∂xμr∂x∂yμr∂ywhichisknownasthehomogeneousscalarwaveequation.Thelatteristhegoverningpartialdifferentialequationofthescatteringproblemunderconsideration.Afiniteelementformulationwasdevelopedearlyinthischaptertosolvethegenericpartialdifferentialequationin(2.1).Comparing(2.196)with(2.1),weobservethattheformerisaspecialcaseofthelatter.Inotherwords,u=Ez1αx=αy=μr(2.197)β=k2ε0rg=0 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28100INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSConsequently,thesameexpressionsderivedinthischapterfortheevaluationoftheelementmatricesandvectorscanbeusedforthefiniteelementsolutionofthehomogeneousscalarwaveequationin(2.196).Asfarasboundaryconditionsareconcerned,thetotaltangentialelectricfieldonperfectlyconductingsurfaces,denotedby1,mustvanish.Inotherwords,onthesurfaceofthecircularcylinderthez-directedtotalelectricfieldiszero,i.e.,Ez=0on1(2.198)Furthermore,thedomainoftheproblemisunboundedsincethisisascatteringproblemand,therefore,thefieldscatteredbythecylindermustcontinuepropagatingtowardinfinitywithoutdisturbance.Thequestionishowdowesimulatenumericallythisundisturbedwavepropagationtoinfinity?Inanelectromagneticanechoicchamber,absorbingconesareplacedonthewalls,ceilingandfloorinordertoreducetheunwantedreflectionsbythesurroundings.Thisineffectsimulatesoutwardandundisturbedwavepropagation.InthecontextoftheFEM,aswellasothernumericalmethodssuchasthefinite-differencetime-domain(FDTD)method[25,26],thenumericaldomainmustbetruncatedatacertaindistanceawayfromtheobject,andontheouterboundaryofthedomain,anABC[27–31]mustbeimposed.ThepurposeofthisABC,whichcorrespondstoapartialdifferentialequationinvolvingtheprimaryunknownvariableandderivativesofit,istoallowthewavetopropagateintheoutwarddirectionwithoutcausinganynumericalreflectionsbacktotheobject.DependingonthehighestdegreeofthederivativesinvolvedinanABC,thereexistdifferentorderssuchasfirst-orderABC,second-orderABC,etc.AstheorderoftheABCincreases,thecomplexityofitincreasesaswell;however,theeffectivenessoftheABCimprovesaccordingly.Otherwaysofproperlyterminatingtheunboundeddomainofascatteringorradiationproblemalsoexistincludingthepopularandwidelyusedperfectlymatchedlayer(PML)[32–35],whichwillnotbediscussedhere.TheadvantageofanABCoverthePMListhatitdoesnotincreasethesizeoftheglobalmatrixandretainssparsitywhereasthePML,althoughitretainssparsity,itincreasesthenumberofunknownsandworsenstheconditionnumberofthematrixsystem.ThemainadvantagesofPMLoverABCincludestraightforwardimplementationandimprovedaccuracy.Forthepresentscatteringproblem,itwasdecidedthatafirst-orderABCbeusedinordertoillustratetheprocessofimposingsuchacondition.Asecond-orderABCcanalsobeusedbut,ofcourse,withabitmorecomplexityinvolvedduringitsimplementation.Thiswillbeleftasanexercisetothereader.AccordingtoBaylissetal.[28,29],afirst-orderABCtobeimposedonacircularboundaryatacertaindistanceawayfromthescattererisgivenby∂Fsca1sca+jk0+F=0(2.199)∂ρ2ρ P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS101whereFscadenotesthescatteredfieldwhetherthatiselectricormagnetic,andρistheradialdistancefromthecenterofthecirculardomaintotheouterboundary.Thefartherawaythiscircularboundaryisplacedfromthescatterer,themoreaccuratethefiniteelementsolutionwillbe[29].Agoodruleofthumbistoplacethefirst-orderABCatleasthalfawavelengthawayfromtheoutersurfaceofthescatterer.Substitutingthescatteredelectricfieldin(2.199)andrealizingthatthetotalelectricfieldcanbewrittenasasumoftheincidentandscatteredfields,scaincscaincEz=Ez+Ez⇒Ez=Ez−Ez(2.200)wecanexpress(2.199)as∂E1∂Einc1zzinc+jk0+Ez=+jk0+Ez(2.201)∂ρ2ρ∂ρ2ρByrealizingthattheρ-directionisnormaltotheouterboundarywheretheABCisbeingimposed,thefirstderivativewithrespecttoρcanbewritteninamoreconvenientformgivenby∂Ez∂Ez∂Ez=aˆx+aˆy·aˆn(2.202)∂ρ∂x∂ywhereaˆnisthenormalunitvectortotheouterboundary.Thus,substituting(2.202)into(2.201),yields∂E∂E1∂Einc∂Einc1zaˆzaˆzaˆzaˆincx+y·aˆn+jk0+Ez=x+y·aˆn+jk0+Ez∂x∂y2ρ∂x∂y2ρ(2.203)Theincidentelectricfieldisgiven,initsgeneralform,byEinc=Ee−jk0(xcosφi+ysinφi)(2.204)z0whereφiistheincidentanglemeasuredwithrespecttothepositivex-axis.Note,however,thataccordingtotheproblemdefinitiontheincidentanglewassettozero,thusEinc=Ee−jk0x(2.205)z0Substituting(2.205)intotheright-handsideof(2.203)andcarryingoutthederivatives,weobtain∂Ez∂Ez1inc1incaˆx+aˆy·aˆn+jk0+Ez=−jk0Ez(aˆx·aˆn)+jk0+Ez∂x∂y2ρ2ρ(2.206) P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28102INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSComparing(2.206)withthemixedboundaryconditionin(2.6),aftersetting1αx=αy==1(2.207)μrsincetheouterboundaryofthedomain(i.e.,2)belongstofreespace—itisconcludedthat1γ=jk0+(2.208)2ρ1incincq=−jk0Ez(aˆx·aˆn)+jk0+Ez(2.209)2ρAsaresult,afirst-orderABConboundary2canbeimposedthroughtheuseofthemixedboundaryconditionin(2.6)andtheassignmentofγandqaccordingto(2.208)and(2.209).TheprocessofimposingamixedboundaryconditionwasexplainedindetailinSection2.5.However,intheimplementationofthemixedboundaryconditioninSection2.5,itwasassumedthatbothγandqareconstantalong2.Asseen,however,from(2.208)and(2.209)onlyγisconstantwhereasqisafunctionofspacecoordinatex.Thus,ithastobetakenintoaccountwhenevaluatingtheintegralof(2.67).Inotherwords,qcannotbetakenoutoftheintegralaswasdoneinSection2.5.Specifically,assumingedge1→2liesontheABCboundary,thecorrespondingentriesofelementvectorpearegivenby1−jkx−e−jk0x21e−jkxe021p=Eqe01100122(k0x21)−1+(jkx+1)e−jk0x21e−jkxe021(2.210)p=Eqe01200122(k0x21)pe=03wherejkye−ye021q0=γ−(2.211)12Ifitwereedge3→1lyingontheABCboundary,insteadofedge1→2,theentriesofelementvectorpewouldtakethefollowingform:−1+(jkx+1)e−jk0x13e−jkxe013p=Eqe03100132(k0x13)pe=0(2.212)21−jkx−e−jk0x13e−jkxe013p=Eqe03300132(k0x13)wherejkye−ye013q0=γ−(2.213)13 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS103Andatlast,ifitwereedge2→3lyingontheABCboundary,theentriesofelementvectorpewouldbepe=011−jkx−e−jk0x32e−jkxe032p=Eqe02200232(2.214)(k0x32)−1+(jkx+1)e−jk0x32e−jkxe032p=Eqe02300232(k0x32)wherejkye−ye032q0=γ−(2.215)23Itisimportanttoreemphasizeatthispointthattheglobalnodesofeachtriangleinthemeshbestoredintheconnectivityinformationarrayinacounter-clockwisedirection,otherwisetheaboveexpressionswillnotbevalid.A2-Dnodalfiniteelementcode,namedFEM2DL--Cyl,waswritteninMatlabtosolvethescatteringproblemdescribedinthissection.Themajorpartsofthecodeincludediscretizationofthedomainusinglineartriangularelements,generationofpropermeshdata,constructionofelementmatricesandvectorsforeachelementinthedomain,assemblyoftheseelementmatricesandvectorsintotheglobalmatrixandright-hand-sidevector,impositionofDirichletboundaryconditionsonthesurfaceofthecircularcylinder,andfinallysolutionofthematrixsystemtoobtainthetotalelectricfieldatthenodesofthedomain.Theelectricfieldatanypointinthedomaincanbeevaluated,asoutlinedinSection2.9,byusingthegoverningsetofinterpolationfunctions.ArelativelycoarsemeshofthedomainisshowninFigure2.18.Theradiusofthecircularconductingcylinder,identifiedinthefigurebytheinnercircularboundary,isλ/2whereastheradiusoftheABCboundary,identifiedinthefigurebytheoutercircularboundary,is3λ/2.Inotherwords,thefirst-orderABCisimposedadistanceofonewavelengthawayfromtheoutersurfaceoftheobject.Figure2.19showsacontourplotofthenormalized(totheamplitudeoftheincidentfield)magnitudeofthetotalelectricfieldinthediscretizeddomain.Thediscretizationsizewasapproximatelyλ/25.Foracceptableresults,thediscretizationsizeshouldnotbelargerthanλ/10whereastheABCboundarymustnotbeplacedcloserthanadistanceofλ/2fromtheoutersurfaceoftheobject.AcomparisonbetweenthefiniteelementsolutionandtheexactanalyticalsolutionisshowninFigure2.20.TheevaluationofthefieldoccurshalfwaybetweentheoutersurfaceoftheconductingcylinderandtheABCboundary.Notethatthediscretizationsizeisroughlyλ/25,andtheABCboundaryisofcirculartypewithradius3λ/2. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28104INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSFIGURE2.18:AcoarsetriangularmeshofthefiniteelementdomainFIGURE2.19:Contourplotofthetotalelectricfieldbasedonafiniteelementsolution P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS1051.6ExactFEM1.41.21.0oE|/z0.8|E0.60.40.20.0060120180240300360(degrees)FIGURE2.20:ComparisonbetweenthetotalelectricfieldcalculatedusingtheexactanalyticalsolutionandtheelectricfieldcomputedusingtheFEMwithfirst-orderABCExercise2.7.StartingfromMaxwellequationsandtheassumptionthat∂Hz=0(2.216)∂zshowthattheinhomogeneousscalarwaveequationfortheTEzpolarizationisgivenby∂1∂Hz∂1∂Hz2∂1∂1++k0μr=−Jy+Jx(2.217)∂xεr∂x∂yεr∂y∂xεr∂yεrExercise2.8.SolvethesamescatteringproblemastheonepresentedinthissectionbutusingaTEz,insteadofaTMz,incidentplanewave.FirstderivetheweakformulationoftheproblemusingtheGalerkinapproachand,then,modifytheexistingMatlabcodetoobtainafiniteelementsolution.Useafirst-orderABCandlineartriangularelements.Exercise2.9.Giventhedefinitionofγandqin(2.208)and(2.209),respectively,showthattheentriesoftheelementright-hand-sidevectorpearegivenby(2.210)–(2.215).2.11HIGHERORDERELEMENTSTheaccuracyoftheFEMcanbeimprovedbyeitherusingafinermeshorhigherorderelements[36,37].Uptothispointinthetext,wehaveusedonlylineartriangularelementsandbilinearquadrilateralelements.Inthissection,wewillconcentrateonderivinghigherorderelements P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28106INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSforquadrilateralsandtriangles.Duetotheirsimplicity,higherorderquadrilateralswillbedevelopedfirst.Specifically,wewilldeveloptheexpressionsfortheshapefunctionsofanine-nodequadraticquadrilateralelement,andtheexpressionsfortheshapefunctionsofasix-nodequadratictriangularelementandaten-nodecubictriangularelement.Theprocedurefollowedtoconstructthesehigherorderelementswillbeillustratedineverydetailincasethereaderwantstousetheunderlinedapproachtodevelopadifferent-orderelementbutofthesametype.2.11.1ANine-NodeQuadraticQuadrilateralElementRememberthatabilinearquadrilateralelementhasfournodes;oneateachofthefourcorners.Aquadraticquadrilateralelementhasanadditionalnodeinthemiddleofeachofthefoursidesandanothernodeinthecenteroftheelement.Inotherwords,thetotalnumberofnodesisnine.Suchanelementisshownplottedinthexy-planeinFigure2.21(a),whereasthecorrespondingsquaremasterelementisshownplottedintheξη-planeinFigure2.21(b).ThemathematicalexpressionsofthegoverningshapefunctionsareobtainedbyimposingtheLagrangeanconditionthattheshapefunctionfornodeimustbeunityatnodeiandzeroatallothernodes,i.e.,1atnodeiNi(ξ,η)=(2.218)0atallothernodesh747334(–1,1)(1,1)6989x68(–1,–1)(1,–1)y515212x(a)(b)FIGURE2.21:(a)Anine-nodequadraticquadrilateralelementplottedinthexy-plane.(b)Squaremasterelementplottedintheξη-plane P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS107Toillustratetheapproachusedinderivingtheseshapefunctions,letusfirstconsidernode1:1atnode1;i.e.,ξ=−1andη=−1N1(ξ,η)=(2.219)0atallothernodesReferringtothesquaremasterelementinFigure2.21(b),itisobservedthatallthenodes,besidesnode1,resideonthehorizontallinesη=0(2.220)η−1=0andtheverticallinesξ=0(2.221)ξ−1=0Therefore,toimposetheconditionthatN1bezeroatallnodesbesidesnode1,itsexpressionmusthavetheformN1(ξ,η)=cξη(ξ−1)(η−1)(2.222)ConstantcisdeterminedbyimposingthesecondconditionwhichsaysthatN1beunityatnode1.Doingso,wehave1=c(−1)(−1)(−2)(−2)1(2.223)⇒c=4Substitutingconstantcinto(2.222),thefinalexpressionforN1isgivenby1N1(ξ,η)=ξη(ξ−1)(η−1)(2.224)4Thesameexactprocedurecanberepeatedfortheremainingeightnodestodeterminethegoverningexpressionsforthecorrespondingshapefunctions.Thesearegivenbelowforthesakeofcompleteness:1N2(ξ,η)=ξη(ξ+1)(η−1)(2.225)41N3(ξ,η)=ξη(ξ+1)(η+1)(2.226)41N4(ξ,η)=ξη(ξ−1)(η+1)(2.227)4 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28108INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS1N5(ξ,η)=η(ξ+1)(ξ−1)(η−1)(2.228)21N6(ξ,η)=−ξ(ξ+1)(η+1)(η−1)(2.229)21N7(ξ,η)=−η(ξ+1)(ξ−1)(η+1)(2.230)21N8(ξ,η)=−ξ(ξ−1)(η+1)(η−1)(2.231)2N9(ξ,η)=(ξ+1)(ξ−1)(η+1)(η−1)(2.232)Exercise2.10.UsingthesameapproachillustratedforshapefunctionN1,showthattheremainingshapefunctionsgoverninganine-nodequadraticquadrilateralelementaregivenby(2.225)–(2.232).2.11.2ASix-NodeQuadraticTriangularElementAquadratictriangularelementcanbeconstructedbyintroducinganadditionalnodeatthemidpointofeachofthethreeedgesofalineartriangle.Inotherwords,aquadratictriangle,shownplottedinthexy-planeinFigure2.22(a),hasatotalofsixnodes:oneateachvertexandoneatthemidpointofeachofthethreeedges.Thecorrespondingmastertriangularelement,plottedintheξη-plane,isdepictedinFigure2.22(b).x=02x-1=0235352h-1=0664h=0y1421x2x+2h−1=0x+h−1=0(a)(b)FIGURE2.22:(a)Asix-nodequadratictriangularelementplottedinthexy-plane.(b)Masterelementplottedintheξη-plane P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS109Assumingisoparametricrepresentation,theprimaryunknownquantityandthespacecoordinatesxandycanbeexpressedaccordingto6u=ueNiii=16x=xeN(2.233)iii=16y=yeNiii=1whereNiaretheLagrangeshapefunctionsdefinedas1atnodeiNi=(2.234)0atallothernodesFornode1,thecorrespondingshapefunction,referringtothemasterelementinFigure2.22(b),musthavetheformN1=c(2ξ+2η−1)(ξ+η−1)(2.235)Evaluating(2.235)atnode1whereξ=0andη=0,constantcisfoundtobec=1(2.236)Thus,thefinalformofN1isN1=(2ξ+2η−1)(ξ+η−1)(2.237)Similarly,1atnode2N2=(2.238)0atallothernodesThesecondrequirementin(2.238)issatisfiedprovidedthattheshapefunctionN2hastheformN2=cξ(2ξ−1)(2.239)TherequirementthatN2=1atnode2mustbeimposedinordertodeterminethevalueofconstantc.Doingso,itisfoundthatc=1(2.240)Thus,thefinalexpressionofN2isgivenbyN2=ξ(2ξ−1)(2.241) P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28110INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSUsingthesamemethodology,itcanbeshownthatN3=η(2η−1)(2.242)N4=−4ξ(ξ+η−1)(2.243)N5=4ξη(2.244)N6=−4η(ξ+η−1)(2.245)Exercise2.11.ShowthattheshapefunctionsN3−N6governingthequadratictriangularelementillustratedinFigure2.22aregivenby(2.242)–(2.245).2.11.3ATen-NodeCubicTriangularElementAten-nodecubictriangularelementcanbeconstructedbyintroducingtwoadditionalnodesateachofthethreeedgesofthetriangleandasinglenodeatthecenteroftheelement,asshowninFigure2.23(a).ThecorrespondingmastertriangleisshowninFigure2.23(b).Theequationsoflinesonwhichthesenodesresideareindicatedinthefigure.Thesewillbeusedtoconstructtheproperexpressionsoftheshapefunctionsspanningtheelement.ReferringtothemastertriangleinFigure2.23(b),itcanbeseenthattheshapefunctionfornode1musthavetheformN1=c(3ξ+3η−1)(3ξ+3η−2)(ξ+η−1)(2.246)x=03310xx--=320=762732h-=0310858106310h-=9492h=0yh1451xx3310xh+−=3320xh+−=xh+−10=(a)(b)FIGURE2.23:(a)Aten-nodecubictriangularelementplottedinthexy-plane.(b)Masterelementplottedintheξη-plane P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS111Allnodes,besidesnode1,resideonthethreeinclinedlineswhoseequationsareindicatedinthefigure.Thecorrespondingshapefunctionfornode1mustbezeroatallthosenodesandthatiswhytheproperformofN1isgivenby(2.246).Todetermineconstantc,itisnecessarytoimposetherequirementthattheshapefunctionbeunityatnode1.Doingso,itcanbeshownthat1c=−(2.247)2Thus,theshapefunctionfornode1isgivenby1N1=−(3ξ+3η−1)(3ξ+3η−2)(ξ+η−1)(2.248)2Inasimilarway,wecanderivetheremainingshapefunctions.Forthesakeofcompleteness,thesearegivenbelow:1N2=ξ(3ξ−1)(3ξ−2)(2.249)21N3=η(3η−1)(3η−2)(2.250)29N4=ξ(3ξ+3η−2)(ξ+η−1)(2.251)29N5=−ξ(3ξ−1)(ξ+η−1)(2.252)29N6=ξη(3ξ−1)(2.253)29N7=−ξη(3ξ−2)(2.254)29N8=−η(3η−1)(ξ+η−1)(2.255)29N9=η(3ξ+3η−2)(ξ+η−1)(2.256)2N10=−27ξη(ξ+η−1)(2.257)Exercise2.12.UsingthesameapproachillustratedforthederivationofN1,showthattheshapefunctionsN2−N10forthecubictriangularelementinFigure2.23aregivenby(2.249)–(2.257).2.12SOFTWARETwoMatlabcodeswerewrittentosolvethetwoapplicationproblemsdiscussedinthischapter.Thefirstcode(FEM2DL--Box),whichuseslineartriangularelements,waswrittentosolvetheelectrostaticBVP,presentedinthischapter,forarectangulardomain.Threeofthewallsaresettoazeropotentialwhereasthetophorizontalwallwassettoanonzeropotential.Theuser,amongotherparameters,hasthefreedomtochangethedensityofthemesh.Thecodecomputestheprimaryunknownquantity,whichistheelectrostaticpotential,andthenumerical P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28112INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICSerrorbasedontheL2norm.Thesecondcode(FEM2DL--Cyl),whichalsouseslinearelements,waswrittentosolvethescatteringproblem,presentedinthischapter,foraconductingcircularcylinder.Theprogramgeneratesatriangularmesh,basedontheuser’schoiceofdiscretizationsize,andcomputesthetotalfieldeverywhereinthedomain.Notethatthetotalfieldisasuperpositionofthescatteredfieldandtheincidentfield.ATMzpolarizationwasconsidered.ThereaderisencouragedtoexecutethesecodesinMatlab,modifycertainparameterssuchasdiscretizationsize,objectdimensions,distancetotheABC,etc.,andvisualizetheresults.WhoeverisinterestedcaneventrytoslightlymodifytheseMatlabcodestosolveothertypesofelectromagneticproblems.BothMatlabcodescanbedownloadedfromthepublisher’sURL:www.morganclaypool.com/page/polycarpou.References[1]J.C.Nedelec,“AnewfamilyofmixedfiniteelementsinR3,”Numer.Methods,vol.30,pp.57–81,1986.[2]M.L.BartonandZ.J.Cendes,“Newvectorfiniteelementsforthree-dimensionalmagneticfieldcomputation,”J.Appl.Phys.,vol.61,pp.3919–3921,1987.doi:10.1063/1.338584[3]Z.J.Cendes,“Vectorfiniteelementsforelectromagneticfieldcomputations,”IEEETrans.Magn.,vol.MAG-27,no.5,pp.3958–3966,Sept.1991.doi:10.1109/20.104970[4]J.P.Webb,“Edgeelementsandwhattheycandoforyou,”IEEETrans.Magn.,vol.MAG-29,pp.1460–1465,Mar.1993.doi:10.1109/20.250678[5]G.Mur,“Thefinite-elementmodelingofthree-dimensionalelectromagneticfieldsusingedgeandnodalelements,”IEEETrans.AntennasPropag.,vol.41,no.7,pp.948–953,July1993.doi:10.1109/8.237627[6]H.Anton,I.Bivens,andS.Davis,Calculaus,7thed.NewYork:Wiley,2002,pp.1075–1090.[7]O.C.ZienkiewiczandR.L.Taylor,TheFiniteElementMethod,4thed.,vol.1:BasicFormulationofLinearProblems.NewYork:McGraw-Hill,1989.[8]D.Zwillinger,HandbookofIntegration.Boston:JonesandBarlett,1992.[9]E.W.CheneyandD.R.Kincaid,NumericalMathematicsandComputing,5thed.Monterey:BrooksCole,2003.[10]M.AbramowitzandI.A.Stegun,HandbookofMathematicalFunctions.NewYork:DoverPublications,1972.[11]G.R.Cowper,“Gaussianquadratureformulasfortriangles,”Int.J.Numer.MethodsEng.,vol.7,pp.405–408,1973.doi:10.1002/nme.1620070316 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28TWO-DIMENSIONALBOUNDARY-VALUEPROBLEMS113[12]R.S.Varga,MatrixIterativeAnalysis,Prentice-Hall,EnglewoodCliffs,NewJersey,1962.[13]D.M.Young,IterativeSolutionofLargeLinearSystems.NewYork:DoverPublications,2003.[14]Y.Saads,IterativeMethodsforSparseLinearSystems.Boston:PWSPublishing,1996.[15]A.M.Bruaset,ASurveyofPreconditionedIterativeMethods.London,UK:Chapman&Hall,1995.[16]A.Chatterjee,J.M.Jin,andJ.L.Volakis,“Edge-basedfiniteelementsandvectorABC’sappliedto3-Dscattering,”IEEETrans.AntennasPropag.,vol.41,no.2,pp.221–226,Feb.1993.doi:10.1109/8.214614[17]J.M.JinandJ.L.Volakis,“Ahybridfiniteelementmethodforscatteringandradiationbymicrostrippatchantennasandarraysresidinginacavity,”IEEETrans.AntennasPropag.,vol.39,no.11,pp.1598–1604,Nov.1991.doi:10.1109/8.102775[18]A.C.Polycarpou,C.A.Balanis,J.T.Aberle,andC.Birtcher,“Radiationandscatteringfromferrite-tunedcavity-backedslotantennas:Theoryandexperiment,”IEEETrans.AntennasPropag.,vol.46,no.9,pp.1297–1306,1998.doi:10.1109/8.719973[19]D.T.McGrathandV.P.Pyati,“Phasedarrayantennaanalysiswiththehybridfiniteelementmethod,”IEEETrans.AntennasPropag.,vol.42,no.12,pp.1625–1630,Dec.1994.doi:10.1109/8.362811[20]A.Chatterjee,J.M.Jin,andJ.L.Volakis,“Computationofcavityresonancesusingedge-basedfiniteelements,”IEEETrans.Microw.TheoryTech.,vol.MTT-40,pp.2106–2108,Nov.1992.doi:10.1109/22.168771[21]J.M.JinandV.V.Liepa,“Applicationofhybridfiniteelementmethodtoelectromag-neticscatteringfromcoatedcylinders,”IEEETrans.AntennasPropag.,vol.36,no.1,pp.50–54,Jan.1988.doi:10.1109/8.1074[22]D.-H.Han,A.C.Polycarpou,andC.A.Balanis,“Hybridanalysisofreflectorantennasincludinghigher-orderinteractionsandblockageeffects,”IEEETrans.AntennasPropag.,vol.50,no.11,pp.1514–1524,2002.doi:10.1109/TAP.2002.803952[23]M.N.O.Sadiku,NumericalTechniquesinElectromagnetics,2nded.BocaRaton:CRCPress,2001.[24]C.A.Balanis,AdvancedEngineeringElectromagnetics.NewYork:Wiley,1989.[25]A.TafloveandS.C.Hagness,ComputationalElectrodynamics:TheFinite-DifferenceTimeDomainMethod,2nded.Boston:ArtechHouse,2000.[26]K.S.KunzandR.J.Luebbers,TheFiniteDifferenceTimeDomainMethodforElectro-magnetics.BocaRatorn:CRCPress,1993.[27]B.EngquistandA.Majda,“Absorbingboundaryconditionsforthenumericalsimulationofwaves,”Math.Comput.,vol.31,pp.329–351,1977. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28114INTRODUCTIONTOTHEFINITEELEMENTMETHODINELECTROMAGNETICS[28]A.BaylissandE.Turkel,“Radiationboundaryconditionsforwave-likeequations,”Com-mun.PureAppl.Math.,vol.33,pp.707–725,1980.[29]A.Bayliss,M.Gunzburger,andE.Turkel,“Boundaryconditionsforthenumericalsolutionofellipticequationsinexteriorregions,”SIAMJ.Appl.Math.,vol.42,pp.430–451,1982.doi:10.1137/0142032[30]A.F.Peterson,“Absorbingboundaryconditionsforthevectorwaveequation,”Microw.Opt.Technol.Lett.,vol.1,pp.62–64,1988.[31]J.P.WebbandV.N.Kanellopoulos,“Absorbingboundaryconditionsforthefiniteelementsolutionofthevectorwaveequation,”Microw.Opt.Technol.Lett.,vol.2,pp.370–372,1989.[32]J.P.Berenger,“Aperfectlymatchedlayerfortheabsorptionofelectromagneticwaves,”J.Comput.Phys.,vol.114,no.2,pp.185–200,Oct.1994.doi:10.1006/jcph.1994.1159[33]Z.S.Sacks,D.M.Kingsland,R.Lee,andJ.-F.Lee,“Aperfectlymatchedanisotropicabsorberforuseasanabsorbingboundarycondition,”IEEETrans.AntennasPropag.,vol.43,no.12,pp.1460–1463,1995.doi:10.1109/8.477075[34]U.PekelandR.Mittra,“Afiniteelementmethodfrequencydomainapplicationoftheperfectlymatchedlayer(PML)concept,”Microw.Opt.Technol.Lett.,vol.9,no.3,pp.117–122,June1995.[35]A.C.Polycarpou,M.R.Lyons,andC.A.Balanis,“Atwo-dimensionalfinite-elementformulationoftheperfectlymatchedlayer,”IEEEMicrow.GuidedWaveLett.,vol.8,pp.30–32,Jan.1997.[36]P.Silvester,“High-orderpolynomialtriangularfiniteelementsforpotentialproblems,”Int.J.Eng.Sci.,vol.7,pp.849–861,1969.doi:10.1016/0020-7225(69)90065-2[37]P.Silvester,“Ageneralhigh-orderfinite-elementwaveguideanalysisprogram,”IEEETrans.Microw.TheoryTech.,vol.MTT-17,pp.204–210,Apr.1969.doi:10.1109/TMTT.1969.1126932 P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28115AuthorBiographiesAnastasisC.PolycarpoureceivedtheB.S.(withsummacumlaude),M.S.,andPh.D.degreesinElectricalEngineeringfromArizonaStateUniversityin1992,1994,and1998,respectively.Duringhisgraduatestudies,hewaswiththeTelecommunicationsResearchCenter(TRC)ofASUwhereheworkedonvariousresearchprojectssponsoredbygovernmentorganizationsandprivatecompaniessuchastheUSNavy,USArmy,Boeing,Sikorsky,andafewmore.Inthesummerof1998,hejoinedtheDepartmentofElectricalEngineeringofASUasanAsso-ciateResearchFacultywhereheperformedresearchonavarietyofsubjectsinthebroadareaofelectromagnetics.WhilebeingatASU,heworkedonthedevelopmentandenhancementofnumericalmethods,inparticulartheFiniteElementMethod(FEM)andtheMethodofMoments(MoM),fortheanalysisofcomplexelectromagneticproblemssuchasmicrowavecircuits,interconnectsandelectronicpackaging,cavity-backedslotantennasinthepres-enceofmagnetizedferrites,andhelicopterelectromagnetics.Hewroteamultipurposethree-dimensionalfiniteelementcodeusingedgeelementstosolvegeometricallycomplexscatteringandradiationproblems.Thecodeutilizesadvancediterativetechniquesinlinearalgebraforthesolutionofextremelylargeindefinitematrixsystems.Dr.Polycarpouhaspublishedmorethan40journalsandconferenceproceedingsandtwochaptersinbooks.HeiscurrentlyanAssociateProfessoratIntercollegeinCyprus.HisresearchareasofinterestincludenumericalmethodsinelectromagneticsandspecificallytheFiniteElementMethodandtheMethodofMoments,antennaanalysisanddesign,electromagnetictheory,andferritematerials. P1:IML/FFXP2:IMLMOBK021-02MOBK021-Polycarpou.clsApril29,200620:28116