Question? Leave a message!




Finite element method for elliptic problems

a mixed finite element method for second order elliptic problems raviart thomas and a weak galerkin finite element method for second-order elliptic problems
Dr.NaveenBansal Profile Pic
Dr.NaveenBansal,India,Teacher
Published Date:25-10-2017
Website URL
Comment
3 The finite element method for elliptic problems 3.1 Difficulties associated with the application of weighted residual methods Although the weighted residual methods introduced in Chapter 2 have been used with success in many areas of physics and engineering, there are certain difficulties which prevent them being more widely used for the solution of practical problems. One obvious problem involves the choice of trial functions. It is clear that for an irregular-shaped boundary, such as in Fig. 3.1, it would in general be impossible to find one function, let alone a sequence of functions, which satisfies every essential boundary condition. Thus, immediately, the class of problems amenable to solution by this method is restricted to those problems with a ‘simple’ geometry. Even if the geometry is suitable and a sequence of functions satisying essen- tial boundary conditions is available, these functions are usually polynomials. It is not difficult to appreciate that, in general, very high-order polynomials would be required to approach the exact behaviour of the unknown over the whole region. A worse situation than this, however, concerns the case of discontinuous material properties. k (x; y) 2 n k (x; y) 1 y x Fig. 3.1 A general curved region in two dimensions in which there is an interface between differing media.72 The Finite Element Method Consider, for example, Poisson’s equation (3.1) −div(k grad u)= f(x,y), where k(x,y) is a function of position which is discontinuous across some interface in the region of interest; see Fig. 3.1. We shall refer to equations such as this as 2 Poisson’s equation. Strictly speaking, Poisson’s equation is−∇ u = f, i.e. k =1. Such a situation arises in electromagnetic theory at the junction of two dielectrics with different permittivities. In these cases ∂u/∂n, the normal derivative of u, would be discontinuous across the interface and polynomials, which are continuously differentiable, would not be suitable for accurate description of this situation. In weighted residual methods, all parts of the region are treated with the same degree of importance, no undue attention being paid to areas which may in fact be of more interest than the rest. Finally, weighted residual methods include coupling between points which are distant from one another, even though this coupling is weak. This yields dense matrices in the final analysis and is costly in terms of computer storage. In this chapter, the ideas behind the weighted residual method are extended so that the above-mentioned difficulties may be overcome. In fact, we shall restrict ourselves to the Galerkin approach; the use of a variational approach will be discussed in Chapter 5. Our examples in Chapter 2 suggest that the Galerkin approach is the most accurate in some sense. However, the importance of the method will be appreciated when we reach Chapters 5 and 7, where we shall see that it yields exactly the same equations as does the variational method and hence we can make specific statements about convergence and accuracy. 3.2 Piecewise application of the Galerkin method We consider an approach in which the region of interest is subdivided into a finite set of elements, connected together at a set of points called the nodes.In each of these elements, the function behaviour is considered individually and then an overall set of equations is assembled from the individual components. These individual components are found by a piecewise application of the Galerkin method. The distinction between element numbering and nodal numbering can some- times lead to confusion, and there is no generally accepted notation. We shall adopt the following: subscripts will refer to nodal numbers and superscripts to element numbers; where it is important to distinguish between them, we shall write i for ‘node i’and i for ‘element i’.The finite element method for elliptic problems 73 Actual boundary A e i Typical node Approximate boundary Typical element Fig. 3.2 Finite element idealization of a two-dimensional region using triangular ele- ments, with a fine mesh-grading in the region of node A. We can see that the difficulties encountered in the pure Galerkin approach are now overcome. Since the boundary geometry may be approximated as closely as required by a polygonal arc, or by polyhedra in three dimensions, even the most irregular boundaries are easily approximated; see Fig. 3.2. In Chapter 4, it will be shown how to use curved elements and get an even better boundary approximation. In the examples that follow, it will be seen that the enforcement of essential boundary conditions presents no problems. Because the trial func- tions are defined in a piecewise manner, discontinuities in normal derivatives over interfaces between different media are easily accounted for and relatively low-degree polynomials may be used to obtain suitable accuracy throughout the region; see Fig. 3.3. Also, grading of the finite element mesh may be used to concentrate on regions of specific interest; see Fig. 3.2 again. Finally, the piecewise nature of the application of the Galerkin method ensures that the influence of one element is limited to those elements actually connected to it, thus uncoupling remote regions. This yields sparse matrices when the overall system of equations is assembled, and hence the computational advantages of sparse matrices become available. 3.3 Terminology The solution of boundary-value problems such as that given by eqns (2.1) and (2.39) frequently represents a quantity associated with a scalar field such as a potential. Consequently, we often refer to such problems as field problems.74 The Finite Element Method u(x) Exact solution Finiteelement x Fig. 3.3 Hypothetical finite element approximation to the solution of a two-point boundary-value problem using linear trial functions. Because the finite element method was developed in its computational form by structural engineers (Argyris 1964, Zienkiewicz and Cheung 1965), the structural terminology has remained in the generalization to field problems. In Section 3.6 we shall develop equations of the form KU = F, where U is a vector of nodal variables, i.e. values of u,∂u/∂x,∂u/∂y etc., eval- uated at the nodes. The number of nodal variables associated with a particular e node is often called the number of degrees of freedom at that node. u is the vector of element nodal variables. K and F are called the overall stiffness matrix and the overall force vector, respectively, and are assembled from element matrices e e k and f . e f is a column vector of known quantities obtained from the non- homogeneous terms in the boundary-value problem under consideration. F is e the column vector of known quantities for the overall system. f is called the e element force vector, and k is called the element stiffness matrix. In some problems the field variables, q, are related to the potential function, u, by equations such as q =−gradu. These relationships then yield, after the finite element analysis, equations of the form q = SU.The finite element method for elliptic problems 75 S is called the overall stress matrix. In recent years, some of the structural terminology has been lost and the matrices K and F are referred to according to the physical properties from which they are derived; for example, in thermal problems, K is sometimes called the conductivity matrix and F is referred to simply as the source term. 3.4 Finite element idealization As in Chapter 2, two-dimensional problems will be considered for illustrative purposes; extensions of the method to three dimensions are, in principle, straight- forward. Consider the field problem (3.2) Lu = f in D with (3.3) Bu = g on C. L and B are differential operators. The finite element method seeks an approx- imation, u ˜(x,y), to the exact solution, u(x,y), in a piecewise manner, the approximation being sought in each of a total of E elements. Thus, in the general e element e, an approximation u ˜ (x,y) is sought in such a manner that outside e, e (3.4) u ˜ (x,y)=0,e=1,...,E. e e For example, in Fig. 3.4, u ˜ (x ,y ) is in general non-zero, but u ˜ (x ,y )=0. A A B B Using eqn (3.4), it follows that the approximate solution may be written as e (3.5) u ˜(x,y)= u ˜ (x,y), e B (x ,y ) B B e A (x ,y ) A A Fig. 3.4 A general element e containing a point A (x ,y ). The point B (x ,y )is A A B B not in e.76 The Finite Element Method where the summation is taken over all the elements. The reason for writing the approximate solution in this way is not immediately apparent, but it makes setting up the overall system of equations a little easier to understand. The first decision to be made is to choose a suitable subdivision of the region into finite elements which sufficiently approximates the boundary geometry. This is very much a matter of choice for the user. More will be said about this when we deal with individual elements in Sections 3.7 and 3.8 and with the isoparametric concept in Chapter 4. It also affects the accuracy of the solution, and this will be discussed in Chapter 6. Once the discretization has been decided on, the next choice is that of the e representation of the element approximation, u ˜ (x,y), in terms of the weighted residual parameters. The most common form of approximation in each element is polynomial approximation. This is probably due to the fact that polynomials are relatively easily manipulated, both algebraically and computationally. Polyno- mials are also attractive from the point of view of the Weierstrass approximation theorem (Wade 1995), which states that any continuous function may be approx- imated arbitrarily closely by a suitable polynomial. The choice of polynomial is a matter for the user, but some guidelines are useful. 1. The number of terms in the polynomial must be equal to the total number of degrees of freedom associated with the element, otherwise the polynomial many not be unique. Thus, for a triangular element with three nodes, one degree of freedom at each node, a three-term polynomial is used, giving e T u ˜ (x,y)= c + c x + c y =1xyc c c . 0 1 2 0 1 2 For an element with four nodes and two degrees of freedom at each node, an eight-term polynomial is used, giving e 2 2 2 2 u ˜ (x,y)= c + c x + c y + c x + c xy + c y + c x y + c xy 0 1 2 3 4 5 6 7 2 2 2 2 T =1xy x xy y x yxy c ,...,c . 0 7 The coefficients c are called generalized coordinates andineachcasethe i approximation is of the form e e e (3.6) u ˜ (x,y)= P (x,y)Δ , e e where P (x,y) is a row vector of linearly independent functions and Δ is a column vector of constants. 2. There should be no preference for either the x or the y direction. This is often referred to by saying that the approximation must have geometrical invariance.The finite element method for elliptic problems 77 3. There are two other requirements, which will be dealt with in Chapter 7. The structural description of these requirements is that the element must be able to exactly reproduce rigid-body motions and constant-strain deformations. Mathematically, these requirements say that to ensure convergence of the method the unknown must be continuous and must be allowed to assume any arbitrary linear form. Although it is not easy to give definite rules which are applicable in all cases, it is in general better not to retain high-order terms at the expense of lower ones. For this reason, complete polynomials are often favoured. Complete polynomials are those in which all possible terms up to any given degree are present. The necessary terms for all possible polynomials up to a complete quintic are shown in Table 3.1. Thus a complete linear polynomial is of the form c + c x + c y 0 1 2 and requires an element with three degrees of freedom to uniquely define c ,c ,c . 0 1 2 A complete cubic polynomial is of the form 2 2 3 2 2 3 c + c x + c y + c x + c xy + c y + c x + c x y + c xy + c y , 0 1 2 3 4 5 6 7 8 9 and this requires an element with ten degrees of freedom to uniquely define c ,...,c . 0 9 Although the form given by eqn (3.6) is of the type used in the pure weighted residual approach of Chapter 2, it is in fact better to use the nodal degrees of freedom as parameters rather than the generalized coordinates. Suppose element ehas p nodes with one degree of freedom at each node; see Fig. 3.5. Using eqn (3.6), the approximation in element eisgiven by e e e u ˜ (x,y)= P (x,y)Δ . Table 3.1 Complete polynomials up to order 5 1 xy 2 2 x xy y 2 2 2 3 x xyxy y 4 3 2 2 3 4 x xyx y xy y 5 4 3 2 2 3 4 5 x xyx y x y xy y78 The Finite Element Method U 2 2 U i i 1 U 1 p U p Fig. 3.5 Element ewith p nodes; U is the nodal variable associated with node i. i Therefore e U =˜ u (x ,y ) i i i  2 2 T = 1 x y x x y y ··· c ...c ,i=1,...,p. i i i i 0 p i i T e Thus the element vector U =U ...U may be written in the form 1 p e e (3.7) U = CΔ , where ⎡ ⎤ 2 1 x y x ··· 1 1 1 ⎢ ⎥ 2 ⎢ 1 x y x ··· ⎥ 2 2 2 ⎢ ⎥ C = . ⎢ ⎥ . . . . ⎢ ⎥ . . . . . . . . ⎣ ⎦ 2 1 x y x ··· p p p Then e −1 e (3.8) Δ = C U . Thus, for this element, eqn (3.6) yields e e −1 e u ˜ (x,y)= P C U , i.e. e e e (3.9) u ˜ (x,y)= N U , e e −1 where N = P C is called the shape function matrix for element e. It relates e e the unknown function u ˜ (x,y) to the nodal variables given by U . Outside e element e, N ≡ 0, so that eqn (3.4) is satisfied. There is a serious drawback to this approach, however, in that the matrix C has to be inverted, which may be computationally costly. Also, the matrix itself is often ill-conditioned and indeed may even be singular; see Exercise 3.1.The finite element method for elliptic problems 79 It is better to obtain eqn (3.9) directly by interpolation throughout the element. This is done by choosing a suitable set of interpolation polynomials e N (x,y), with the property that if (x ,y ) are the coordinates of node j,then j j i e (3.10) N (x ,y )= δ,i,j =1,...,p. j j ij i The shape function matrix is then given by  e e e e (3.11) N (x,y)= N (x,y) N (x,y) ··· N (x,y) . 1 2 p The shape functions must be such that conditions 1, 2 and 3 are satisfied. e e Condition 3 leads to relationships between N and N ; see Exercise 3.2. i j It is often helpful in setting up the shape functions, as well as in the following analysis, to use a set of local coordinates, say (ξ,η), for each element rather than the global coordinates (x,y). This will be illustrated in the examples in Sections 3.5–3.8. Indeed, when we deal with isoparametric elements in Chapter 4, this procedure is essential. At this stage, the discretization and the manner of element approximation have been decided. The next step in the procedure is to use the Galerkin method to set up the linear equations for the nodal variables U,i=1,...,n, and perhaps i the derivatives (∂u/∂x) , (∂u/∂y) , etc. This particular step involves two distinct i i sets of operations; setting up the equation for the individual elements, followed by the assembly of the overall system of equations. It is best demonstrated by way of examples, and these will be illustrated in the following sections. The solution of the linear equations yields the overall vector of unknowns T U=U U ... U 1 2 n as the approximate solution of eqn (3.1) at the nodal points. The approximate solution at non-nodal points is given by interpolation in each element as e e e u ˜(x,y)= u ˜ (x,y)= N (x,y)U . e e Finally, once the nodal values have been obtained, there may be field variables to detemine; for example, if u(x,y) is the velocity potential for a fluid flow, then the velocity vector at any point (x,y)isgiven by q(x,y)=−gradu. In particular, in element e,  −∂/∂x e e q = u ˜ −∂/∂y  e −∂N /∂x e = U , e −∂N /∂y80 The Finite Element Method i.e. e e e (3.12) q (x,y)= S U , where e e e T (3.13) S =−∂N /∂x − ∂N /∂y is the element stress matrix. For problems involving material anisotropy described by a tensor κ,itmay be shown (see Exercise 3.8) that T e e e (3.14) S = κ−∂N /∂x − ∂N /∂y . So far, we have dealt with the approximation expressed by eqn (3.5) in an element-by-element manner, since the approach is to develop the region as a set of elements. It is, however, more helpful to consider the approximation (3.5) in e a node-by-node manner. The shape functions N (x,y) have the property that i e N (x,y)=0 if i∈ e, so that we can define a nodal function w (x,y)which is i i local to node i,given by e (3.15) w (x,y)= N (x,y). i j ej So, instead of using an element-by-element approximation as in eqn (3.5), we use a node-by-node approximation n (3.16) u ˜(x,y)= w (x,y)U , j j j=1 where n is the number of nodes in the finite element approximation. The set of nodal functions is thus a linearly independent set of basis functions for the approximation. 3.5 Illustrative problem involving one independent variable The reason for the success of the finite element method is that it may be applied to problems of great complexity. However, to use such problems for illustrative purposes tends to obscure the underlying ideas. In this section, a simple problem involving an ordinary differential equation is solved. This is not to suggest that the finite element method is the best method for solving such problems; in fact, there are other numerical methods available which are superior. However, this particular problem involves only a small amount of algebra and the ‘mechanics’ of the method come through very well.The finite element method for elliptic problems 81 2 3 1 U U U 2 3 1 x 1 x =1 1 x =– x =0 2 2 Fig. 3.6 Finite element idealization for the two-point boundary-value problem (3.17), showing the elements 1 and 2 and the global node numbering 1, 2, 3. Example 3.1 Consider the two-point boundary-value problem  −u =2, 0x 1, (3.17)  u(0) = u (1) = 0. Consider a two-element discretization of 0, 1 with three nodes, as shown in Fig. 3.6, and suppose that there is one degree of freedom at each node. In each element there are just two nodal variables, and hence the interpola- tion polynomials for each element must be linear. Consider an element ewith midpoint x and length h, suppose that the nodes associated with element e m have local labels A, B,andthat ξ is a local coordinate as shown in Fig. 3.7 and given by 2 (3.18) ξ = (x− x ). m h The shape function matrix is e e e N (ξ)=N (ξ) N (ξ) , A B where e e N (−1) = N (1) = 1 A B and e e N (1) = N (−1) = 0; A B see Fig. 3.8. AB u u A B 2 ξ = -(x ¾ x ) m h x m h ξ=1 ξ=–1 Fig. 3.7 Element e, showing the local node labels A,B and the local coordinate ξ.82 The Finite Element Method 1 e e N (ξ) N (ξ) A B ξ 1 -1 e e Fig. 3.8 Linear shape functions N (ξ)and N (ξ). A B The shape functions are easily recognized as the Lagrange interpolation polynomials e 1 e 1 (3.19) N (ξ)= (1− ξ),N (ξ)= (1 + ξ). A B 2 2 Since the interpolation is linear, this is often called a linear element. Then, in element e, e e e (3.20) U = N U with T e (3.21) U =U U . A B Therefore the overall finite element approximation as given by eqn (3.5) is 2 e (3.22) u ˜ = u ˜ . e=1 In this case we have an essential Dirichlet condition at node 1 and a natural Neumann condition at node 3. Consequently, we must take U = u(0) and we 1 have just two unknown nodal variables, U and U . 2 3 The corresponding nodal functions are 1 1 2 2 w (x)= N (x),w (x)= N (x)+ N (x),w (x)= N (x), 1 2 3 A B A B where the global coordinate x is related to the local coordinate ξ by eqn (3.18), and (3.23) u ˜(x)= w (x)U + w (x)U + w (x)U . 1 1 2 2 3 3 The Galerkin formulation is obtained by using a weighted residual approach, taking the basis functions w and w as the weighting functions; we need only 2 3 those two functions associated with the unknowns U and U : 2 3The finite element method for elliptic problems 83  1  (−u ˜ − 2)w dx=0,i=2,3. i 0 We integrate by parts, to reduce the order of derivative required:   1 1 1    −u ˜ w + u ˜ w dx− 2w dx=0,i=2,3. i i 0 i 0 0 We now write down the two equations as follows:    1 1 1 du ˜ du ˜ dw 2 i=2: − w + dx− 2w dx=0. 2 2 dx dx dx 0 0 0 1 2 Now, w (0) = 0 since N (0) = 0 and N (x)≡ 0 in element 1, and w (1) = 0 2 2 B A 1 2 since N (x)≡ 0 in element 2 and N (1) = 0. Hence the first term is zero and B A we have    1/2 1 1 du ˜ dw du ˜ dw 2 2 dx + dx− 2w dx=0, 2 dx dx dx dx 0 1/2 0 i.e.     1/2 1/2 dw dw dw 1 2 2 U + U dx− 2w dx 1 2 2 dx dx dx 0 0     1 1 dw dw dw 2 3 2 + U + U − 2w dx=0; 2 3 2 dx dx dx 1 1/2 2    1 1 1 du ˜ du ˜ dw 3 i=3: − w + dx− 2w dx=0. 3 3 dx dx dx 0 0 0  Now, w (x)≡ 0 in element 1 and the Neumann condition u (0) = 0 is a natural 3 condition. Hence the first term is zero and we have     1 1 dw dw dw 2 3 3 U + U dx− 2w dx=0. 2 3 3 dx dx dx 1/2 1/2 It follows then that our two equations are of the form 1 1 2 2 1 2 k U +(k + k )U + k U = f + f , 1 2 3 21 22 11 12 2 1 (3.24) 2 2 2 k U + k U = f , 2 3 21 22 2 where the element stiffness matrices and force vectors are given respectively by  e e dN dN j e i k = dx ij dx dx e (3.25)  e and f = 2N dx. i i e84 The Finite Element Method Now, U = u(0), a known value, so it follows that eqns (3.24) are sufficent to find 1 U and U , and our problem is solved. However, since the idea will be helpful in 2 3 a more general setting, we develop the equation associated with the weighting function w (x). This equation has the form 1 1 1 1 (3.26) k U + k U = f , 1 2 11 12 1 using the notation of eqn (3.25). Consequently, the overall system of equations may be written as (3.27) KU = F, where ⎡ ⎤ ⎡ ⎤ 1 1 1 k k 0 f 11 12 1 ⎢ ⎥ ⎢ ⎥ 1 1 2 2 1 2 K = k k + k k and F = f + f . ⎣ ⎦ ⎣ ⎦ 21 22 11 12 2 1 2 2 2 0 k k f 21 22 2 Before proceeding further, some remarks regarding the overall stiffness matrix may be made here, since they are generally applicable. 1. It is clearly symmetric, as indeed are the element stiffness matrices. 2. K = K = 0, showing that there is no coupling between nodes 1 and 3, i.e. 13 31 the only coupling occurs between nodes associated with the same element. It is not difficult to see, then, that for a large number of elements K becomes sparse and banded. 3. K is singular. This, perhaps, is not so obvious. In structural terms, this simply says that K allows rigid-body motions, i.e. the structure is not fixed. This then suggests how the ‘singularity’ in K may be interpreted from the point of view of the solution of the boundary-value problem (3.17). Fixing a structure requires prescribing certain displacements, usually equal to zero, and this is equivalent to enforcing a Dirichlet boundary condition. By adding eqn (3.26) to the set of equations (3.24) we are, in essence, not enforcing the essential boundary condition. Before we see how to remove the singularity in K,the element stiffness and forces will be evaluated. Using the local coordinate ξ in the integrations in eqn (3.25) yields  1 2 dN 2 dN h i j k = dξ, i,j =1,2. ij h dξ h dξ 2 −1 Whence, using eqn (3.19), 1 1 k = k =,k = k =− . 11 22 12 21 h hThe finite element method for elliptic problems 85 Also,  1 h f = 2N dξ, i=1,2, i i 2 −1 whence f = f = h. 1 2 A convenient way to write the element matrices is to label the rows and columns according to the nodal variable with which they are associated, as shown below: 12 23   1 1 1 −1 1 1 −1 2 1 2 k = , k = , h −11 2 h −11 3 12 23   T T 1 2 f = h 11 , f = h 11 . The overall stiffness and force matrices are thus 12 3 ⎡ ⎤ 1 −10 1 1 ⎣ ⎦ K = −11+1 −1 2 h 0 −11 3 and 123  T F = h 11+1 1 . Notice that in K the contribution to the 2,2 position, and in F the contri- bution to the 2,1 position, is made up from two terms, one from each element. This is typical of the way that the overall matrices are assembled in general. If a node is associated with more than one element, then the contribution to the relevant parts of the overall matrix is merely a matter of addition of the corresponding terms. This will be discussed in Section 3.6. Notice also in this case that the element matrices have been defined in terms of h, the element length. Consequently, it is possible to solve this problem using elements of different lengths without the necessity of altering the analysis; see Exercise 3.3. 1 The equations (3.27) now become, with h = , 2 ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 −10 U 1 1 1 1 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ (3.28) −12 −1 U = 2 . 2 12 2 0 −11 U 1 386 The Finite Element Method The next step in the procedure is to solve eqn (3.28) for the unknown nodal variables. As was seen earlier, the overall stiffness matrix is singular, so that a solution of the equations is not possible as they stand. However, the boundary conditions have still to be imposed, and we return to the original set of equa- tions (3.24), which we obtain by removing row 1 from K and setting U = u(0) 1 where appropriate to obtain 4U − 2U =1, 2 3 (3.29) 1 −2U +2U = , 2 3 2 3 which gives U = , U =1. 2 3 4 To express the solution u ˜(x), it is usually more convenient to return to the element-by-element form, eqns (3.5) and (3.20):  T 1 1 3 u ˜ (x)= 1− ξ 1+ ξ 0 2 4   3x 2 1 = , since ξ = x− in element 1, 1 4 2 2 and  T 2 1 3 u ˜ (x)= 1− ξ 1+ ξ 1 2 4   x+1 2 3 = , since ξ = x− in element 2. 1 4 2 2 Then eqn (3.22) yields the finite element solution as  1 3x/2, 0≤ x≤ , 2 u ˜(x)= 1 (x+1)/2, ≤ x≤ 1. 2 This is compared with a four-element solution and the exact solution in Fig. 3.9. Now consider a four-element solution of eqn (3.17), with the discretization 1 as shown in Fig. 3.10. In this case, with all elements of length , the element 4 stiffness matrices are 12 23   1 −1 1 1 −1 2 1 2 k =4 , k =4 , −11 2 −11 3 34 45   1 −1 3 1 −1 4 3 4 k =4 , k =4 , −11 4 −11 4 and the element force vectors are     1 1 1 2 1 3 1 4 1 1 1 1 1 2 3 4 f = , f = , f = , f = , 4 4 4 4 1 2 1 3 1 4 1 5The finite element method for elliptic problems 87 u(x) 1 Exact solution 0.9 2elements 0.8 4elements 0.7 0.6 0.5 0.4 0.3 0.2 0.1 x 0 0 0.2 0.4 0.6 0.8 1 Fig. 3.9 Comparison of the two- and four-element solutions of eqn (3.17) with the exact 2 solution u (x)=2x− x . 0 U U 1 2 U 3 4 U 5 U 1 3 2 4 5 x 1 2 3 4 1 1 3 x =– x =– x =– x = 0 x = 1 4 2 4 Fig. 3.10 Four-element idealization for the two-point boundary-value problem (3.17). whence the overall system is ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ 1 −1000 U 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 2 −100 U 2 2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 4 2 −10 U = 2 . 3 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ 4 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 2 −1 U 2 4 symmetric 1 U 1 5 Enforcing the homogeneous Dirichlet boundary condition U = 0 then yields 1 the set of equations 1 8U −4U = , 2 3 2 1 −4U +8U −4U = , 2 3 4 2 1 −4U +8U −4U = , 3 4 5 2 1 −4U +4U = . 4 5 488 The Finite Element Method The solution is 7 3 15 U = ,U = ,U = ,U =1. 2 3 4 5 16 4 16 Then, as in the two-element case, the approximate solution throughout the interval 0,1 is found from eqns (3.20) and (3.22). It is easily seen to be given by ⎧ 7 1 ⎪ x, 0≤ x≤ , ⎪ 4 4 ⎪ ⎪ ⎪ ⎪ 1 1 1 ⎨ (10x+1), ≤ x≤ , 8 4 2 U(x)= 1 1 3 ⎪ ⎪ (6x+3), ≤ x≤ , ⎪ 8 2 4 ⎪ ⎪ ⎪ ⎩ 1 3 (x+3), ≤ x≤ 1. 4 4 This solution is compared with the two-element solution and the exact solution in Fig. 3.9. There is also a comparison with three other finite element solutions in Table 3.2. It should be noted here that the right-hand side in eqn (3.17) is a constant, so that the element force vector need only be found for element e. In general, it would be necessary to obtain the element force vectors separately for each element. The exact solution of eqn (3.17) is 2 u (x)=2x− x . 0 From Fig. 3.9 it is easy to see that the refined mesh, with four elements, gives a better approximation than does the original coarse two-element mesh. Indeed, this is the basis of the method; the more elements, the better the solution. In practice, a suitable number of elements is chosen to give satisfactory accuracy. Notice in this case that the convergence to the exact solution is monotonic; this is due to the fact that the refined mesh contains the original mesh as a subset. This is discussed in Chapter 7. Notice also that the mesh with graded elements Table 3.2 A comparison of the four finite element solutions to the problem in Example 3.1 and Exercise 3.3 x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 2 elements 0.15 0.3 0.45 0.6 0.75 0.8 0.85 0.9 0.95 1 4 elements 0.175 0.35 0.5 0.625 0.75 0.825 0.9 0.95 0.975 1 5 elements 0.18 0.36 0.5 0.64 0.74 0.84 0.91 0.96 0.99 1 10 elements 0.19 0.36 0.51 0.64 0.75 0.84 0.91 0.96 0.99 1 4graded elements 0.188 0.356 0.5 0.625 0.75 0.8 0.85 0.9 0.95 1 Exact solution 0.19 0.36 0.51 0.64 0.75 0.84 0.91 0.96 0.99 1The finite element method for elliptic problems 89 gives better results than the equivalent ungraded mesh in the region in which u 0 is changing most rapidly. There is one other point that may be made by reference to this example, which relates to the case of a non-homogeneous boundary condition. For a non- homogeneous Dirichlet condition, for example u(0) = 1, the overall set of equations is treated slightly differently. Consider the two- element approximation equations (3.28) as set up earlier. For the homogeneous Dirichlet condition U = 0, the row and column corresponding to U in the 1 1 stiffness matrix were deleted, leaving eqn (3.29); in the non-homogeneous case this is not done. The first equation in eqn (3.28) is replaced by U = 1, leaving 1 a reduced set of two equations. For a non-homogeneous Neumann boundary condition, for example  u (1) = 1, 1  the term −(du/dx ˜ )w is now replaced by −u (1)w (1), since the Neumann 3 3 0 condition is a natural boundary condition and w (x)≡ 0 in element 1. This 3 T givesrisetoaterm0 − 1 on the left-hand side of the reduced set of equations as follows: ⎡ ⎤    1 0 −12 −1 2 1 ⎣ ⎦ +2 U = , 2 2 −1 0 −11 1 U 3 which yields      4 −2 U 1 −1 0 2 = − 2 − , 1 −22 U 0 −1 3 2 9 whence U = , U = 3; see the similar procedure adopted in the classical 2 3 4 Rayleigh–Ritz method in Section 2.6. The resulting finite element solution is then  5 1 1+ x, 0≤ x≤ , 2 2 u ˜(x)= 3 3 1 + x, ≤ x≤ 1. 2 2 2  3  We note that u ˜ (1) = , a relatively poor approximation to u (1) = 1, but we 2 have used only two elements. This simple problem was deliberately chosen and worked through in detail to illustrate how the method works, step-by-step. Since the procedure for more general problems is identical with this one, the basic steps will be listed here and the generalization to problems involving partial differential equations will be presented in the next section.90 The Finite Element Method 1. Subdivide the region of interest into a finite number of subregions, the finite elements. In Example 3.1, the elements were simple subintervals of the interval 0,1. In two- and three-dimensional problems, there is a variety of elements to choose from, for example triangles and rectangles in two dimensions, and tetrahedra and rectangular bricks in three dimensions. 2. Choose nodal variables and shape functions so that the function behaviour throughout each element may be obtained. In Example 3.1, the nodal variables were simply the function values U , and the shape functions were chosen to i e give a linear variation for u ˜ . It is not necessary that only function values be determined at the nodes; derivatives may also be taken as nodal variables and suitable shape functions chosen; see Section 4.4. 3. Obtain the element stiffness and force matrices, and the stress matrices if field variables are required. In Chapter 5 we shall show how variational methods may be used as an alternative to the Galerkin method in certain cases. 4. Assemble the overall system of equations from the individual element matrices. 5. Enforce the essential boundary conditions. In Example 3.1, there were two boundary conditions: an essential homogeneous Dirichlet condition, which was enforced at this stage, and a natural homogeneous Neumann condition, which was handled automatically. 6. Solve the overall system of equations. In Example 3.1, four elements only were used and the resulting system of equations was easily solved by hand; for more equations, efficient computational methods must be used. 7. Compute further results. In many practical examples, field variables must be e found from the function u ˜ . These are obtained using the overall stress matrix. Example 3.2  x −u = e ,  u(0) = 1,u (1) = 2. We illustrate the use of a spreadsheet to solve this problem with five equal elements. Firstly, with the usual notation, we have element stiffness matrices  1 1−1 e k = . h −11 The element force vectors are given by (see Exercise 3.5)  (1 + 2/h)sinh(h/2)− cosh(h/2) H 1 x x m m e = e . H 2 (1− 2/h)sinh(h/2) + cosh(h/2)