New version page

UK EE 699 - 2D FEM Solution of the Helmholtz Equation

This preview shows page 1-2 out of 6 pages.

View Full Document
View Full Document

End of preview. Want to read all 6 pages?

Upload your study docs or become a GradeBuddy member to access this document.

View Full Document
Unformatted text preview:

2D FEM Solution of the Helmholtz Equation (TMz)^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Description:2D solution of the Helmholtz wave equation for the TMz polarization. Arectangular domain is assumed. The side walls are zero Neumann boundaries.The bottom wall is a Dirichlet boundary. The top wall is a Neumann boundary.This version assumes a homogeneous media. The domain is discretized into Ntriangular elements. First order Nodal elements are assumed. The results areplotted and the error is estimated against the exact solution.^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^f300 106⋅16:= c0 2.998 108⋅:= λc0f:= λ 15.989=VDBC01:=k02 π⋅λ:=Read in the coordinates of the nodes from an external file:rNodenodeList.txt:= triNodeelemList.txt:=numNodes rows rNode():= numTris rows triNode():=numNodes 36=numTris 50=rN rNodeT:=Graph the Mesh:meshMapx index 0←node triNodeij,1−←tmpindexrNodenode 0,←index index 1+←j02..∈fori 0 numTris 1−..∈fortmp:= meshMapy index 0←node triNodeij,1−←tmpindexrNodenode 1,←index index 1+←j02..∈fori 0 numTris 1−..∈fortmp:=0 0.5 100.51meshMapymeshMapxRead in the boundary nodes:The node # (1st column) and the Dirichlet boundary ID (2nd Column) are read in.Note that a pre-processor has sorted all nodes so that nodes on Dirichlet boundaries are at theend of the list. boundaryNodesboundaryNodeList.txt:= numBoundaryNodes rows boundaryNodes():=numBoundaryNodes 12=DirchletboundaryNodes index 0←tmpindex 0,boundaryNodesi0,← boundaryNodesi1,1=iftmpindex 1,boundaryNodesi1,← boundaryNodesi1,1=ifindex index 1+← boundaryNodesi1,1=ifi 0 numBoundaryNodes 1−..∈fortmp:=numDBoundaryNodes rows DirchletboundaryNodes():= numDBoundaryNodes 6=numIntNodes numNodes numDBoundaryNodes−:= numIntNodes 30=nodeBctmpi0←i 0 numNodes 1−..∈fornode DirchletboundaryNodesi0,1−←bc DirchletboundaryNodesi1,←tmpnodebc←i 0 numDBoundaryNodes 1−..∈fortmp:=NeumannNodes index 0←tmpindexboundaryNodesi0,← boundaryNodesi1,2=ifindex index 1+← boundaryNodesi1,2=ifi 0 numBoundaryNodes 1−..∈fortmp:=Compute the unitary vectors for the triangles:ii 0 numTris..:=aSub1Vectors index 0←node0 triNodei0,1−←node1 triNodei1,1−←asub1 rNnode1〈〉rNnode0〈〉−←tmpji,asub1j←j02..∈fori 0 numTris 1−..∈fortmp:= aSub2Vectors index 0←node0 triNodei0,1−←node2 triNodei2,1−←asub2 rNnode2〈〉rNnode0〈〉−←tmpji,asub2j←j02..∈fori 0 numTris 1−..∈fortmp:=rootG2d index 0←a1 aSub1Vectorsi〈〉←a2 aSub2Vectorsi〈〉←rg2dia1 a2×←i 0 numTris 1−..∈forrg2d:=Compute the reciprocal unitary vectors for the triangles:zhat001⎛⎜⎜⎝⎞⎠:=aSup1Vectors index 0←asub2 aSub2Vectorsi〈〉←asup1asub2 zhat×rootG2di←tmpji,asup1j←j02..∈fori 0 numTris 1−..∈fortmp:= aSup2Vectors index 0←asub1 aSub1Vectorsi〈〉←asup2zhat asub1×rootG2di←tmpji,asup2j←j02..∈fori 0 numTris 1−..∈fortmp:=Compute the Global K-matrix:Ktmpij,0←i 0 numNodes 1−..∈forj 0 numNodes 1−..∈fornodejtriNodeij,1−←j02..∈forasup1 aSup1Vectorsi〈〉←asup2 aSup2Vectorsi〈〉←asup0 asup1− asup2−←asupj0,asup0j←asupj1,asup1j←asupj2,asup2j←j02..∈foraq asupq〈〉←as asups〈〉←Keaq as⋅2Teqs,−⎛⎜⎝⎞⎠rootG2di⋅←tmpnodeq()nodes,tmpnodeq()nodes,Ke+← nodesnumIntNodes<ifs02..∈for nodeqnumIntNodes<iftmpnodeq()nodeq,1← otherwiseq02..∈forblankspace 0←i 0 numTris 1−..∈fortmp:=K0123401234567891011121.998 -1.001 0 0 0-1.001 3.997 -1.001 0 00 -1.001 3.997 -1.001 00 0 -1.001 3.997 -1.0010 0 0 -1.001 3.9970 0 0 0 -1.001-0.5-4-5.147·10 0 0 00 -1.001-4-5.147·10 0 00 0 -1.001-4-5.147·10 00 0 0 -1.001-4-5.147·100 0 0 0 -1.0010000000000=Tek0212112121211212121⎛⎜⎜⎜⎜⎜⎜⎝⎞⎟⎟⎟⎟⎠⋅:=Compute the Right-Hand-Side:brhsj0←j 0 numNodes 1−..∈fornodejtriNodeij,1−←bcjnodeBcnodej1−←j02..∈forasup1 aSup1Vectorsi〈〉←asup2 aSup2Vectorsi〈〉←asup0 asup1− asup2−←asupj0,asup0j←asupj1,asup1j←asupj2,asup2j←j02..∈foraq asupq〈〉←as asups〈〉←Keaq as⋅2Teqs,−⎛⎜⎝⎞⎠rootG2di⋅←rhsnodeqrhsnodeqKe VDBCbcs⋅−← nodesnumIntNodes≥ifs02..∈for nodeqnumIntNodes<ifrhsnodeqVDBCbcq← otherwiseq02..∈forblankspace 0←i 0 numTris 1−..∈forrhs:=Compute the Solution:b001234567891011121314150.5011.0011.0011.0011.0010.50000000000=c K1−b⋅:=Compute the Error Relative to the Exact Solution:ymax 1:=ExactyrN1i,←tmpicos k0 y ymax−()⋅[]cos k0 ymax⋅()←i 0 numNodes 1−..∈fortmp:=Exact001234567891011121314151.0294577781951.0294577781951.0294577781951.0294577781951.0294577781951.0294577781951.0525601401491.0525601401491.0525601401491.0525601401491.0525601401491.0525601401491.0691644621151.0691644621151.0691644621151.069164462115= c001234567891011121314151.0294720981391.029470855861.0294532785351.0294305997461.0294130009841.0294117193181.0526214997981.0526058002961.0525595105681.0525041119881.0524577964351.0524420296861.0693277009421.0692763203671.0691801131091.069074146415=Error Exact c−:=meanError1numIntNodes0numIntNodes 1−iErrori∑=⋅:=meanError 1.968


View Full Document
Loading Unlocking...
Login

Join to view 2D FEM Solution of the Helmholtz Equation and access 3M+ class-specific study document.

or
We will never post anything without your permission.
Don't have an account?
Sign Up

Join to view 2D FEM Solution of the Helmholtz Equation and access 3M+ class-specific study document.

or

By creating an account you agree to our Privacy Policy and Terms Of Use

Already a member?