EE474 Computational Electromagnetics

Computer Assignment #3-2D FEM

Employ the finite element method for determining the capacitance per unit length of a cylindrical coaxial cable. Compute the capacitance from energy considerations (We =(CV^2)/2), or from charge considerations (C =QV). Normalize the calculations to e0. Let the inner conductor radius be a = 1 cm, and the outer conductor radius be b = 2 cm. Use the 2D finite element program given in Jin or Silvester and Ferrari as guidance in developing your finite element program. In particular, the subroutines MESHIN, ELMATR, and ELEMBD in Silvester and Ferrari will be very helpful. A 2D triangular patch mesh generator will be provided. Consider the following cases:

1. A cylindrical coaxial cable with an air dielectric filling. Investigate the errors associated with solution approximation (mesh density), and geometry approximation. Calculate the error between the exact solution and the finite element solution in all cases.

(a) CASE 1: Errors due to approximation with the basis or interpolation functions. Consider at least 5 meshing densities from very course to very fine. Strive for a relatively uniform mesh in each case. This is essentially a convergence type study.

The mesh density in my solution is based purely on the number of horizontal blocks (Elements (H)) used in my solution. For simplicity this number will be included in results (since it is what is ultimately used to get them) as well as the total number of elements (Total # of elements is proportional to mesh density.)

(b) CASE 2: Errors due to boundary approximation. Approximate the cylindrical boundaries with four, ten, and twenty segments. Use a fine mesh. Here several nodes may be put on each linear segment to get a sufficiently fine mesh that the capacitance value for a given boundary approximation has converged.

In my solution, I have no direct control over the boundary approximation. Increasing the variable that controls mesh density does, however, increase the number of segments along the boundary. (See program discussion.) Since I have no control over this, the variation in the number of segments on the boundary with the meshing density will have to suffice.

I will start with 4 “Elements(H)”, since it is the absolute minimum that produces a valid mesh. Here is the mesh and the result. There were 6 inner segments, 18 outer segments, 36 total elements. The solution yielded a capacitance of 82.47 (pF). For a charge of 1C the inner perimeter voltage was 1.2E10volts. (not shown on graph), while the outer was set to 0 volts.

Note: The background color (non meshed area) in the intensity plot has no meaning. It was set to be a nice contrasting shade.

The same exact problem with 25 “Elements (H)” has the results shown in the following graph. More elements are deemed unnecessary since this is very close to the exact solution. There were 82 inner segments, 166 outer segments, 1032 total elements. The solution yielded a capacitance of 80.40 (pF/m). For a charge of 1C the inner perimeter voltage was 1.243E10 volts, while the outer was set to 0 volts.

The exact solution from the mathcad document for this system is 80.261pF/m. It is used in the percentage error calculations. The variation with both the number of inner and outer segments, and meshing density (total elements) is summarized in the following table.

Elements(H)	Total Elements	Inner Segments	Outer Segments	V difference	Capacitance	% Error
4	36	6	18	1.2E10	82.47(pF/m)	2.752%
8	140	26	50	1.22E10	81.96	2.117%
16	466	50	106	1.241E10	80.58	0.397%
20	682	66	130	1.242E10	80.49	0.285%
25	1032	82	166	1.243E10	80.40	0.173%

2. Inhomogeneously filled cables. Compare the exact and computed solutions for the capacitance per unit length. Work toward a mesh that will give an error of less than 5% with the exact solution.

(a) Layered media.

The exact solution for this was C=202.556(pF/m) from the mathcad document.

Elements(H)	Total Elements	Inner Segments	Outer Segments	V difference	Capacitance	% Error
25	1284	82	166	4.927e9	202.91(pF)	0.175%

(b) Partially filled.

The exact solution for this was C=240.782(pF/m) from the mathcad document.

Elements(H)	Total Elements	Inner Segments	Outer Segments	V difference	Capacitance	% Error
25	1032	82	166	4.145e9	241.20(pF)	0.174%

(Note: The gray scale on the intensity graph is quantized so it doesn’t go quite to zero.)

Submit your program and results with a brief summary, including:

1. A summary of the formulation including development of element matrices etc., and the capacitance calculation.

The formulation of the A element matrix (3x3) is done exactly as indicated in the derivation. The incorporation of the A element matrix into the global matrix is done via simply matching global node reference to the local node reference and adding it in. I.E. If the element was A1,2 and the global nodes associated with local nodes 1 and 2 were 30 and 412 then the element would be added in at A30,412 where this A is the global matrix.

The D global matrix was determined directly via the manner discussed in the derivation. In general it is D[i]=D`(len[I]+len[I+1])/2. Where the only points to be determined where the inner boundary, since the outer boundary would be zero (see derivation) and the inside elements would cancel out. (This is one of the reasons for a requirement of counter clockwise numbered local nodes.)

The capacitance calculation is simply Q/V when all is said and done. (Q over the potential difference between the outer and inner conductor.)

2. Results and Conclusions

The main results for the middle split case is that for all combinations of er1 and er2 for the middle split the graph of the results is exactly the same except its simply scaled by some unknown (before calculation) constant. This is a direct result of there being no charge in the conductor and the voltage being constant (p.e.c) on the inside and outside of the coaxial cable.

The main result for the layered media case is exactly what we would expect. The potential tends to have a smooth variation in the region with lower Er and is more constant with higher Er. This is no great surprise since we expect less variation in voltage in better conductors.

The main overall conclusion is that, surprisingly, it doesn’t necessarily take that many elements to get reasonable solutions. Even the smallest mesh that my generation program could make without creating an error produced a result within 3%.

The main program design and implementation conclusion is that the hardest step tends to be mesh generation. There were many exceptions and what you might almost call fuzzy logic to cover irritating little cases where the mesh did not generate correct either due to finite precision or a node landing on a boundary and confusing the system.

Currently the circular boundaries approximated inscribe the circular boundaries created. If a way was designed to compensate so that the middle of the elements in the circular boundary created approximately align with the true boundary then a small increase in accuracy might be possible. It is questionable if this increase is worth the program complexity required.

Program Discussion

The program is designed around the idea of creating a regular 2D mesh using equilateral triangles. An example of such a mesh shown here.

Note the blue block. That is the segment of 4 triangles I chose to “fill space” with. After establishing this initial mesh the boundary conditions are imposed and the mesh is adjusted to match those boundary conditions. Where needed new triangles are created. The revised mesh is shown here.

In the interests of time and since the problem requirements assumed that I would be using a different mesh generation routine that was provided, no further discussion of the workings of this mesh generator will be provided. Suffice it to say it creates a properly ordered mesh with all the local nodes lined up counter clockwise etc, etc. It should, however, be noted that this mesh generation routine is directly a part of the solve program making it much more flexible for solving this or similar problems than the one provided which ran on a unix platform.

Program Design

1.0 Overall Design

The entire program from mesh formulation, variable inputs, boundary conditions, linear equation solving is all designed and implemented in native Labview code. The overall program design is based on the manipulation of data in a cluster array. This array contains all the node specific information as well as the information specific to an element. A few additional parameters that are not written to after initial setting are stored in global variables to simplify the code. The cluster array begins with default values and a series of sub VI’s (functions) are used to manipulate these variables. A summary of the important steps follows.

1. Mesh Construction—Fill entire domain with triangles and determine the region for every point. The program is designed around the idea of creating a regular 2D mesh using equilateral triangles. An example of such a mesh shown here. Note the blue block. That is the segment of 4 triangles I chose to “fill space” with.

2. Calculate Slope and intercept for all sides of all triangles in domain.
3. Calculate the points where the triangles intersect region boundaries. Using this info create new triangles as needed and recalculated the slopes and intercepts of those new triangles. After establishing this initial mesh the boundary conditions are imposed and the mesh is adjusted to match those boundary conditions. Where needed new triangles are created.
4. Remove regions that are not part of the computational domain. In this case, both the triangles inside the inner conductor boundary and those outside the outer conductor boundary. A typical revised mesh is shown here.

5. Assign global nodes Pass 1.
6. Using global nodes remove any non triangles. Non triangles result from the from finite precision placing two nodes so close together than they are mathematically indistinguishable to the computer. These elements are found by counting the number of unique nodes per element.
7. Calculate center points and angles and force CCW. The calculation of the center points of the triangles allows one to determine the relative angle with respect to the center point. This allows one to reorder any elements so that the local node numbering is CCW in the event it was assigned incorrectly in step 3. (No preference for node order was made in step 3. It was simpler to just double check them all and correct as necessary.)
8. Assign global nodes Pass 2 (final pass). This corrects for the node removal from step 6.
9. Calculate the A matrix using the formulas in the derivation. Store the alpha, betas, and gammas for later use.
10. Calculate the D vector using the formulas in the derivation.
11. Apply a fixed voltage to the inner boundary. Force the outer boundary to be a constant, but unknown voltage (See Derivation).
12. Send to the linear equations solver and get a result.
13. Plot the results. Plotting is accomplished by running points that fill the space of a square that is large enough to completely surround the element triangle. Determination of whether or not a point is inside or on the triangle is accomplished as follows. First a given side is selected and the y-intercept is calculated using that slope. Comparing this intercept with the intercept from the side chosen yields a relationship for points that are inside the triangle with respect to that side. The intercept of the point under test is calculated with the same slope to determine if it is inside the triangle with respect to the one side. The point is checked against the other two sides similarly and if it is inside with respect to all three points it is used to plot a point in the solution. The constants used in this (z=k1+k2*x+k3*y) are determined in the manner discussed in the derivation using the alpha, beta, and gamma information previously saved. The process repeats for all elements until finished. (The other triangle corner could have been used instead of the center point.)
14. Determination of the capacitance is based on the charge placed on the outer (or inner) conductor/divided by the potential difference between those two conductors. (Negative signs have no meaning and thus the absolute value is taken.)