同济大学有限元Chap06(S)

6.2 Finite Element Approximation
✧ We now consider the finite element approximation of the solution of our governing equation. We begin by placing Eq.6.8 in the following matrix form : ⎰⎰⎰=Φ-Φ-⎪⎪⎭
⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂Φ∂∂Φ∂⎥⎦⎤⎢⎣⎡⎥⎦⎥⎢⎣⎢∂Φ∂∂Φ∂V S y x V dS q QdV dV y x k k y x 000 *δδδδ (6.13) or, in abbreviated notation,
⎣⎦[]{}⎰⎰⎰=Φ-Φ-Φ∇Φ∇V S
V qdS QdV dV R 0δδδ (6.14) where
{}⎪⎪⎭
⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂Φ∂∂Φ∂=Φ∇y x (6.15)
✧ Recall from Chapter 5 that within a single element,
⎣⎦{}Φ=ΦN y x ),( (6.16)
and that {}⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧ΦΦΦ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂Φ∂∂Φ∂=Φ∇K J I K J
I
K J I y N y N y N x N x N x N y x (6.17)
abbreviated as
{}[]{}Φ'=Φ∇N (6.18)
✧ Also recall that
⎣⎦{}Φ=ΦδδN and {}[]{}Φ'=Φ∇δδN
✧ We now substitute the above FEM approximations into Eq. 6.14 and integrate. We must do so, of course, element by element . We begin with the volume integral associate with the stiffness term, where, for our generic element e , we have
⎣⎦[]{}⎣⎦[][][]{}⎣⎦[]{}
ΦΦ=Φ''Φ=Φ∇Φ∇⎰⎰e T V V K dV N R N dV R e
e δδδ (6.19)
✧ In expanded notation the element stiffness matrix is
[]dV y N y N y N x N x N x N k k y N x N y N x N y N x N K K J I K J I y x V K K
J J I I e e ⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂⎥⎦⎤⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂=⎰00 (6.20) a 3×3 matrix. The other volume integral is
⎣⎦⎣⎦⎣⎦{}e T V V Q dV y x Q N QdV e e
Φ=Φ=Φ⎰
⎰δδδ),( (6.21) where {}⎰⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧=e V K J I e dV y x Q N N N Q ),(
(6.22)
✧ The total contribution from element e to the volume integration in Eq.6.8 is
⎣⎦[]{}⎣⎦{}e e Q K Φ-ΦΦδδδ (6.23)
where the Φ terms are those associated with the given element.
✧ We now evaluate the surface integral in Eq. 6.14. Figure 6.2 illustrates a finite
element approximation of Φ along a boundary. Consider a typical segment of S , S s , between nodal point m and nodal point n . Along m -n , Φ is a linear function of S . We may therefore use the shape functions developed for our one-dimensional problems and write
Figure 6.2 Segment of boundary.
⎣⎦{}e S n m S S N l S l S S Φ=⎭
⎬⎫⎩⎨⎧ΦΦ⎥⎦⎥⎢⎣⎢⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛-=Φ 1)( (6.24) where l S is the length of the segment.
✧ The surface integral along S S can now be written as
⎣⎦{}dS S q N dS q S S
S S S )(*⎰⎰Φ=Φδδ (6.25)
which gives us ⎣⎦{}S S Q dS q S ⎰
Φ=Φδδ* (6.26)
where {}dS S q l S l S Q Q Q S S S S n m S )(1⎰⎪⎪⎭
⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛-=⎭⎬⎫⎩⎨⎧= (6.27) ✧ The governing interpretation for Q m and Q n is that they equal the first moment of
area of q (S ) about their respective nodal point, divided by the length of the segment. Hence, we place a point Q at each node that gives the same moment as the distributed q . Also, note that
⎰=+S
S n m qdS Q Q (6.28)
Figure 6.3 Integration of surface integral
✧ We now have the contributions to Eq.6.8 from a typical volume element e and a typical surface segment s . We write these as
⎣⎦[]{}⎣⎦{}e e e Q K J Φ+ΦΦ=δδδ (6.29)
and
⎣⎦{}S S Q J Φ-=δδ (6.30)
✧ To obtain the total variation, we need to sum all such contributions; hence,
∑∑+=S e J J J δδδ (6.31)
⎣⎦[]{}⎣⎦{}Q K J e Φ-ΦΦ=δδδ (6.32)
where [K ] is the matrix sum of all the [K ]e , and {Q} is the sum of all the {Q}e and {Q}s .
✧ We now seek values of {Φ} for which δJ = 0 for all {δΦ}. This requires
[]{}{}Q K e =Φ (6.33)
which is our final governing finite element equation .
6.3 Programming Preliminaries
✧Before we look at the program to create and solve Eq.6.33, several programming
details must be discussed. Other details will given with the listing of the code, and yet others you will be able to obtain from the example data and functions included with the test problem.
6.3.1 Element Identification
✧As in the codes for one-dimensional analyses, the core of our new code will be
the integration of the stiffness matrix, element by element. To accomplish this, we must identify the location and geometry of each element. This is done by identifying the three nodes associated with each element numbers are circled. The nodes associated with each of these elements can be identified using the following two-dimensional array:
Figure 6.4 A simple two-dimensional mesh.
✧The nodal point coordinates will be located in two arrays, XORD and YORD.
6.3.2 Integration of [K]e and {Q}e matrices
✧Because our shape functions ⎣⎦N are linear in x and y, the gradient terms in
N'are constants. If we also assume that the coefficients k x and k y, in [R] are ⎣⎦
constant, we can write
[][][][][][][][][][]V N R N dV N R N dV N R N K T
V T
T
V e e
e ''=''=''=⎰⎰ (6.34) ✧ Note again the notation we have chosen to use for our region o
f integration. In the above, V is actually the area of the element. The notation indicates that we are assumin
g it is an area associated wit
h a volume of unit thickness .
✧ The derivatives of the shape functions will be stored in the arrays DNDX and DNDY , where, for example;
x
N ∂∂=2DNDX(2) ✧ The element stiffness matrix, [K ], will be stored in the 3×3 array S(J, K), defined in terms of the DNDX and DNDY arrays as
VOL
*DNDY(K))*RY *DNDY(J) DNDX(K)*RX *(DNDX(J)K) S(J,+= (6.35) where
YJ))/2.0*(XK - YK)*((XJ VOL = (6.36)
✧ The expression for {Q }e is
{}{}QdV N Q e
V e ⎰= (6.37) ✧ If we assume Q constant within each element, it can be shown that all three components of {Q }e are equal and have the value
{}⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧=VOL/3.0*Q VOL/3.0*Q VOL/3.0*Q e Q (6.38)
6.3.3 The Global Stiffness Matrix
✧ Placement of an element S matrix into the global SK matrix is best explained by an example. Let
[]⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=KK S KJ S KI S JK S JJ
S JI
S IK S IJ S II
S S e e e e e e e e e e (6.39)
be a typical stiffness matrix. If, say, it is for element in Fig.6.4, the nodes {I , J , K } correspond to the nodal numbers {6, 5, 2} as obtained from the NP array. Thus, the nine terms of this matrix would appear in the global stiffness matrix as:
✧The matrix is symmetric and will also be banded; thus, appropriate storage will
be used.
6.3.4 User-Provided INCLUDE Codes
✧As our finite element codes become more complex, it is necessary to use more
user-written INCLUDE codes to define specific problems. In the code to follow, poisson.m, there are two such INCLUDE codes: INITIAL.m and COEF.m.
✧INITIAL.m
In codes wire.m and ode2.m, the initialization of variables is made by loading the MESH data file. As we now consider two-dimensional problems, it becomes more difficult to define the initial values of the variables this way. It is much more convenient to initialize variables using an INCLUDE code. This will appear in poisson.m after the mesh data has been loaded and after the general initialization has taken place. In INITIAL.m the user must define NPBC, PHI, and Q for all nodes that differ from the values given in the general initialization (all zeros). This will be accomplished through the use of a nodal point code defined in the MESH data file exactly as the NPBC array was defined in wire.m and ode.m. However, in poisson.m the NPcode array, as it will be called, is not used. Only in INITIAL.m will it be used to identify which nodes should have which boundary conditions. The NPcode values are assigned by the user and entered with the mesh data. These values can then be used in the user’s own INITIAL.m code to assign specific boundary conditions. An example of both a data file and an INITIAL.m code will be given with the test problem.
✧COEF.m
The three parameters, k x(x, y), k y(x, y), and Q(x, y), are defined for each element with the COEF.m INCLUDE. Available for the definitions will be the coordinates of the centroid of the element, XC and YC, as well as the element number, the NP array, etc. The parameters are assumed constant within the element for the integrations. This provides a good approximation for the integration of the stiffness matrix where all other terms in the integrand are constant. However, it is not as accurate a procedure for integration of the right-hand side. Care must be
taken to make elements small enough that Q(x, y) is approximately constant within each element, or at least approximately linear.
6.3.5 The MESHo, NODES, and NP Data files
✧Three files define the mesh used for the analysis. The first, MESHo, contains
three separate lines:
NUMNP = Number of nodal points
NUMEL = Number of elements
NNPE = Number of nodes per element
✧The number of nodes per element for poisson.m is 3. For codes to be considered
later, higher-order elements will be used for which the number of nodes per element will be greater than 3.
✧The second file, NODES, contains one line per node. Each line gives for its
particular node.
XORD YORD NPcode
That is, the x and y coordinates of the node and the node identification number. ✧The third file, NP, contains the connectivity array for each element explained
earlier. An example of such an array was given in Section 6.3.1.
✧One of the problems encountered in using finite element codes for
two-dimensional analyses is how to create the mesh data. It is not unusual for there to be several thousand nodes and elements. Obviously, the user will not want to create the corresponding input files by hand. For this purpose there are mesh-generating codes that create this data and require the user only to enter data related to boundary locations and number of elements desired. One such program, mesh.m, is presented in Appendix D. It creates all three of the data files just described. In addition, it provides the user with an easy method for assigning NPcode values to all boundary nodes.
6.3.6 Bandwidth and Bandwidth Reduction
✧ A problem encountered with two-dimensional analyses is the numbering of the
nodes to minimize the storage needed for the stiffness matrix. For our previous one-dimensional problems, the bandwidth was either 2 (for symmetric matrices) or 3 (for nonsymmetric matrices).
✧With two-dimensional problems, the needed bandwidth will depend on the way
the nodes have been numbered. Consider the two schemes shown in Figure 6.5.
For the first, the largest numerical difference between any two nodes associated with a single element is 3. This would require a bandwidth of 4 if the matrix were symmetric. For the second numbering scheme, the largest numerical difference is 8, which would require a bandwidth of 9.
Figure 6.5Two numbering schemes.
✧You might think that (1) the difference is not large and (2) it is easy to see how to
number the nodes to get the smaller bandwidth. However, for complex meshes, (1) the difference between a good numbering scheme and a poor numbering scheme can result in very large difference in bandwidth requirements, and (2) determining how to number the nodes to minimize the bandwidth can be difficult.
✧The optimal numbering scheme depends on the connectivity of the nodes, and
this is completely specified by the NP array. Hence, computer codes are written that use this array as input to determine a new numbering scheme that will provide a narrow bandwidth. One such code, newnum.m, is described and given in Appendix D.
✧The output from this code is the data file NWLD, which stands for new(old), that
is, the new node number in terms of the old number. This array of numbers can be used as input to any finite element code to create a new numbering system that will produce a stiffness matrix having a narrow bandwidth. Whether program newnum.m or another code is used, NWLD must be in the working directory. The last entry in this code must be the bandwidth that the new numbering creates (see the example problem).
6.4 Program poisson.m
✧Program poisson.m is the finite element code for solving Eq. 6.1, Poisson’s
equation. In addition to poisson.m, you will need the following files in your working directory:
✧MESHo.txt is the data file that defines NUMNP, NUMEL, and NNPE.
NODES.txt is the data file that defines, for each node, its XORD, YORD, and NPcode values. NP.txt is the connectivity array. NWLD.txt gives the new nodal numbering for bandwidth reduction. INITIAL.m is the user’s INCLUDE code for
the initialization of variables to define the specific problem being analyzed.
COEF.m is the user’s INCLUDE code to define the parameters (coefficients) k x, k y, and Q. sGAUSS.m is the supplied function for solving the banded, symmetric equations by Gauss elimination.
Flow Chart
6.4.2 Code
6.4.3 Test Problem
For our test problem, we again use the fact that whenever our finite element approximation can duplicate an exact solution, it will do so. The problem we have selected is defined as follows
where Φ satisfies
055=⎥⎦
⎤⎢⎣⎡∂Φ∂∂∂+⎥⎦⎤⎢⎣⎡∂Φ∂∂∂y y x x The exact solution to this problem is
x )8/6(9-=Φ
where x is measured from the left side of the rectangle. Clearly, our piecewise linear approximation can duplicate this function; thus, poisson.m should produce the exact solution.
✧ The first step in obtaining our finite element solution is to select a mesh. We
choose a mesh with 16 elements and 15 evenly spaced nodes. The mesh and the
corresponding data are as follows
✧ Here we have defined our NPcode values to indicate the side of the mesh on
which the node appears, and only on the sides where there is need to specify a boundary condition that differs from the default values. We have chosen to number the sides counterclockwise, starting with the lower boundary. Thus, nodes 11, 6, and 1 are given NPcode numbers equal to 4, and nodes 15, 10, and 5 are given NPcode numbers equal to 2. However, there is no need to indicate that nodes 2, 3, and 4 on side 1 or that nodes 12, 13, and 14 are on side 3, because there is no need to change any of the default values associated with these nodes. The important concept to understand is that the NPcode values are for your use and only your use in your INITIAL.m code. Use them there to identify the nodes for which you wish to assign boundary conditions that differ from the default
boundary conditions. See the following example for the INITIAL.m code. Note, however, if you wish to assign boundary values in INITIAL.m without reference to the NPcode values, that is your option; you can simply give all nodes an NPcode value of zero and use, for example, the nodal coordinates to determine what boundary conditions to assign.
✧The NWLD values have been assigned to make the numbering begin at the upper
left corner and increase download, similar to the example given earlier. Sketch the mesh with the new numbers and verify the NWLD data as well as the IB value shown. Keep in mind that only poisson.m uses the new numbers, not the user. In all data files, and in INITIAL.m, the original node numbers should be used.
✧The user INCLUDE codes for the test problem are:
✧In INITIAL.m, we search through all NPcode values to identify those equal to 2
(node is on side 2) and equal to 4 ( node is on side 4). We again emphasize that these numbers simply point out to the user wants them to point out. We could have used 3 and 9 to indicate the specified PHI values for these nodes and have accomplished the very same purpose.
✧The test problem has conastant coeffcients; hence, our COEF.m code is fairly
simple. In problems with more complicated coefficients, for example, those that depend on the x and y coordinates, functions can be written in this code. Keep in mind when you write this code that you have use of all variables that have been
calculated in poisson.m up to the line where COEF appears. This includes the coordinates of the element’s centroid, XC and YC, as well as the gradients of the shape functions and other quantities.
✧At the conclusion of the analysis, the resulting Φvalues are saved in your
working directory in a file named PHI. For the test problem, it will contain the exact answers, which are
9.0 7.5 6.0 4.5 3.0
9.0 7.5 6.0 4.5 3.0
9.0 7.5 6.0 4.5 3.0
shwon in the order corresponding to the mesh nodes.
✧Three auxiliary codes can be used in conjunction with poisson.m to help prepare
input data and to help analyze the solution. They can be found in Appendix D and are as follows
mesh.m Generates mesh files: MESHo.txt, NODES.txt, NP.txt
newnum.m Generates data file NWLD.txt
topo.m Generates plots of PHI values。