Numerical Solution of Partial Differential Equations

Numerical Solution of Partial Differential Equations by Gordon C. Everstine 21 January 2010 Copyright c 2001–2010 by Gordon C. Everstine. All rights reserved. This book was typeset with LATEX 2ε (MiKTeX). Preface These lecture notes are intended to supplement a one-semester graduate-level engineering course at The George Washington University in numerical methods for the solution of partial differential equations. Both finite difference and finite element methods are included. The main prerequisite is a standard undergraduate calculus sequence including ordinary differential equations. In general, the mix of topics and level of presentation are aimed at upper-level undergraduates and first-year graduate students in mechanical, aerospace, and civil engineering. Gordon Everstine Gaithersburg, Maryland January 2010 iii Contents 1 Numerical Solution of Ordinary Differential 1.1 Euler’s Method . . . . . . . . . . . . . . . . 1.2 Truncation Error for Euler’s Method . . . . 1.3 Runge-Kutta Methods . . . . . . . . . . . . 1.4 Systems of Equations . . . . . . . . . . . . . 1.5 Finite Differences . . . . . . . . . . . . . . . 1.6 Boundary Value Problems . . . . . . . . . . 1.6.1 Example . . . . . . . . . . . . . . . . 1.6.2 Solving Tridiagonal Systems . . . . . 1.7 Shooting Methods . . . . . . . . . . . . . . . . . . . . . . . . 1 2 3 4 5 6 8 9 10 10 2 Partial Differential Equations 2.1 Classical Equations of Mathematical Physics . . . . . . . . . . . . . . . . . . 2.2 Classification of Partial Differential Equations . . . . . . . . . . . . . . . . . 2.3 Transformation to Nondimensional Form . . . . . . . . . . . . . . . . . . . . 11 11 14 15 3 Finite Difference Solution of Partial Differential Equations 3.1 Parabolic Equations . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Explicit Finite Difference Method . . . . . . . . . . . . 3.1.2 Crank-Nicolson Implicit Method . . . . . . . . . . . . . 3.1.3 Derivative Boundary Conditions . . . . . . . . . . . . . 3.2 Hyperbolic Equations . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 The d’Alembert Solution of the Wave Equation . . . . 3.2.2 Finite Differences . . . . . . . . . . . . . . . . . . . . . 3.2.3 Starting Procedure for Explicit Algorithm . . . . . . . 3.2.4 Nonreflecting Boundaries . . . . . . . . . . . . . . . . . 3.3 Elliptic Equations . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Derivative Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 16 16 18 20 21 21 25 26 27 31 33 4 Direct Finite Element Analysis 4.1 Linear Mass-Spring Systems . . . . . . . . . 4.2 Matrix Assembly . . . . . . . . . . . . . . . 4.3 Constraints . . . . . . . . . . . . . . . . . . 4.4 Example and Summary . . . . . . . . . . . . 4.5 Pin-Jointed Rod Element . . . . . . . . . . . 4.6 Pin-Jointed Frame Example . . . . . . . . . 4.7 Boundary Conditions by Matrix Partitioning 4.8 Alternative Approach to Constraints . . . . 4.9 Beams in Flexure . . . . . . . . . . . . . . . 4.10 Direct Approach to Continuum Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 34 36 36 37 38 41 42 43 44 45 v Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Change of Basis 5.1 Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Examples of Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Isotropic Tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 53 54 57 6 Calculus of Variations 6.1 Example 1: The Shortest Distance Between Two Points 6.2 Example 2: The Brachistochrone . . . . . . . . . . . . 6.3 Constraint Conditions . . . . . . . . . . . . . . . . . . 6.4 Example 3: A Constrained Minimization Problem . . . 6.5 Functions of Several Independent Variables . . . . . . . 6.6 Example 4: Poisson’s Equation . . . . . . . . . . . . . 6.7 Functions of Several Dependent Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 60 61 63 63 64 66 66 7 Variational Approach to the Finite Element Method 7.1 Index Notation and Summation Convention . . . . . . 7.2 Deriving Variational Principles . . . . . . . . . . . . . . 7.3 Shape Functions . . . . . . . . . . . . . . . . . . . . . . 7.4 Variational Approach . . . . . . . . . . . . . . . . . . . 7.5 Matrices for Linear Triangle . . . . . . . . . . . . . . . 7.6 Interpretation of Functional . . . . . . . . . . . . . . . 7.7 Stiffness in Elasticity in Terms of Shape Functions . . . 7.8 Element Compatibility . . . . . . . . . . . . . . . . . . 7.9 Method of Weighted Residuals (Galerkin’s Method) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 67 68 70 73 76 79 80 81 83 8 Potential Fluid Flow With Finite Elements 8.1 Finite Element Model . . . . . . . . . . . . . . . . . . . . 8.2 Application of Symmetry . . . . . . . . . . . . . . . . . . 8.3 Free Surface Flows . . . . . . . . . . . . . . . . . . . . . 8.4 Use of Complex Numbers and Phasors in Wave Problems 8.5 2-D Wave Maker . . . . . . . . . . . . . . . . . . . . . . 8.6 Linear Triangle Matrices for 2-D Wave Maker Problem . 8.7 Mechanical Analogy for the Free Surface Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 86 87 89 90 91 94 95 Bibliography 97 Index 99 List of Figures 1 2 3 4 5 1-DOF Mass-Spring System. . . . . . . . . . . . . Simply-Supported Beam With Distributed Load. . Finite Difference Approximations to Derivatives. . The Shooting Method. . . . . . . . . . . . . . . . Mesh for 1-D Heat Equation. . . . . . . . . . . . vi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 6 11 17 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Heat Equation Stencil for Explicit Finite Difference Algorithm. . Heat Equation Stencil for r = 1/10. . . . . . . . . . . . . . . . . Heat Equation Stencils for r = 1/2 and r = 1. . . . . . . . . . . Explicit Finite Difference Solution With r = 0.48. . . . . . . . . Explicit Finite Difference Solution With r = 0.52. . . . . . . . . Mesh for Crank-Nicolson. . . . . . . . . . . . . . . . . . . . . . . Stencil for Crank-Nicolson Algorithm. . . . . . . . . . . . . . . . Treatment of Derivative Boundary Conditions. . . . . . . . . . . Propagation of Initial Displacement. . . . . . . . . . . . . . . . Initial Velocity Function. . . . . . . . . . . . . . . . . . . . . . . Propagation of Initial Velocity. . . . . . . . . . . . . . . . . . . . Domains of Influence and Dependence. . . . . . . . . . . . . . . Mesh for Explicit Solution of Wave Equation. . . . . . . . . . . Stencil for Explicit Solution of Wave Equation. . . . . . . . . . . Domains of Dependence for r > 1. . . . . . . . . . . . . . . . . . Finite Difference Mesh at Nonreflecting Boundary. . . . . . . . . Finite Length Simulation of an Infinite Bar. . . . . . . . . . . . Laplace’s Equation on Rectangular Domain. . . . . . . . . . . . Finite Difference Grid on Rectangular Domain. . . . . . . . . . The Neighborhood of Point (i, j). . . . . . . . . . . . . . . . . . 20-Point Finite Difference Mesh. . . . . . . . . . . . . . . . . . . Laplace’s Equation With Dirichlet and Neumann B.C. . . . . . Treatment of Neumann Boundary Conditions. . . . . . . . . . . 2-DOF Mass-Spring System. . . . . . . . . . . . . . . . . . . . . A Single Spring Element. . . . . . . . . . . . . . . . . . . . . . . 3-DOF Spring System. . . . . . . . . . . . . . . . . . . . . . . . Spring System With Constraint. . . . . . . . . . . . . . . . . . . 4-DOF Spring System. . . . . . . . . . . . . . . . . . . . . . . . Pin-Jointed Rod Element. . . . . . . . . . . . . . . . . . . . . . Truss Structure Modeled With Pin-Jointed Rods. . . . . . . . . The Degrees of Freedom for a Pin-Jointed Rod Element in 2-D. Computing 2-D Stiffness of Pin-Jointed Rod. . . . . . . . . . . . Pin-Jointed Frame Example. . . . . . . . . . . . . . . . . . . . . Example With Reactions and Loads at Same DOF. . . . . . . . Large Spring Approach to Constraints. . . . . . . . . . . . . . . DOF for Beam in Flexure (2-D). . . . . . . . . . . . . . . . . . . The Beam Problem Associated With Column 1. . . . . . . . . . The Beam Problem Associated With Column 2. . . . . . . . . . DOF for 2-D Beam Element. . . . . . . . . . . . . . . . . . . . . Plate With Constraints and Loads. . . . . . . . . . . . . . . . . DOF for the Linear Triangular Membrane Element. . . . . . . . Element Coordinate Systems in the Finite Element Method. . . Basis Vectors in Polar Coordinate System. . . . . . . . . . . . . Change of Basis. . . . . . . . . . . . . . . . . . . . . . . . . . . Element Coordinate System for Pin-Jointed Rod. . . . . . . . . vii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 17 17 18 19 20 20 21 23 24 24 25 26 26 27 28 30 31 31 32 32 33 34 34 35 36 37 37 39 39 40 40 41 42 43 44 44 45 45 46 46 50 50 51 56 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 Minimum, Maximum, and Neutral Stationary Values. . . . . . Curve of Minimum Length Between Two Points. . . . . . . . . The Brachistochrone Problem. . . . . . . . . . . . . . . . . . . Several Brachistochrone Solutions. . . . . . . . . . . . . . . . . A Constrained Minimization Problem. . . . . . . . . . . . . . Two-Dimensional Finite Element Mesh. . . . . . . . . . . . . . Triangular Finite Element. . . . . . . . . . . . . . . . . . . . . Axial Member (Pin-Jointed Truss Element). . . . . . . . . . . Neumannn Boundary Condition at Internal Boundary. . . . . Two Adjacent Finite Elements. . . . . . . . . . . . . . . . . . Triangular Mesh at Boundary. . . . . . . . . . . . . . . . . . . Discontinuous Function. . . . . . . . . . . . . . . . . . . . . . Compatibility at an Element Boundary. . . . . . . . . . . . . . A Vector Analogy for Galerkin’s Method. . . . . . . . . . . . . Potential Flow Around Solid Body. . . . . . . . . . . . . . . . Streamlines Around Circular Cylinder. . . . . . . . . . . . . . Symmetry With Respect to y = 0. . . . . . . . . . . . . . . . . Antisymmetry With Respect to x = 0. . . . . . . . . . . . . . Boundary Value Problem for Flow Around Circular Cylinder. The Free Surface Problem. . . . . . . . . . . . . . . . . . . . . The Complex Amplitude. . . . . . . . . . . . . . . . . . . . . . Phasor Addition. . . . . . . . . . . . . . . . . . . . . . . . . . 2-D Wave Maker: Time Domain. . . . . . . . . . . . . . . . . Graphical Solution of ω 2 /α = g tanh(αd). . . . . . . . . . . . . 2-D Wave Maker: Frequency Domain. . . . . . . . . . . . . . . Single DOF Mass-Spring-Dashpot System. . . . . . . . . . . . viii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 60 61 63 64 70 70 72 74 75 78 81 82 83 86 87 87 88 88 89 90 91 91 93 93 95 1 Numerical Solution of Ordinary Differential Equations An ordinary differential equation (ODE) is an equation that involves an unknown function (the dependent variable) and some of its derivatives with respect to a single independent variable. An nth-order equation has the highest order derivative of order n:
(1.1) f x, y, y ′ , y ′′ , . . . , y (n) = 0 for a ≤ x ≤ b, where y = y(x), and y (n) denotes the nth derivative with respect to x. An nth-order ODE requires the specification of n conditions to assure uniqueness of the solution. If all conditions are imposed at x = a, the conditions are called initial conditions (I.C.), and the problem is an initial value problem (IVP). If the conditions are imposed at both x = a and x = b, the conditions are called boundary conditions (B.C.), and the problem is a boundary value problem (BVP). For example, consider the initial value problem  mü + ku = f (t) (1.2) u(0) = 5, u̇(0) = 0, where u = u(t), and dots denote differentiation with respect to the time t. This equation describes a one-degree-of-freedom mass-spring system which is released from rest and subjected to a time-dependent force, as illustrated in Fig. 1. Initial value problems generally arise in time-dependent situations. An example of a boundary value problem is shown in Fig. 2, for which the differential equation is  x FL  EIu′′ (x) = M (x) = x − Fx (1.3) 2 2  u(0) = u(L) = 0, where the independent variable x is the distance from the left end, u is the transverse displacement, and M (x) is the internal bending moment at x. Boundary value problems generally arise in static (time-independent) situations. As we will see, IVPs and BVPs must be treated differently numerically. A system of n first-order ODEs has the form  ′ y1 (x) = f1 (x, y1 , y2 , . . . , yn )     y ′ (x) = f2 (x, y1 , y2 , . . . , yn ) 2 (1.4) ..  .    ′ yn (x) = fn (x, y1 , y2 , . . . , yn ) ✲ k ✁❆❆✁✁❆❆✁✁ m ❡ ❡ u(t) ✲ f (t) Figure 1: 1-DOF Mass-Spring System. 1 F (N/m) ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❄ ❅ ❅ ✛ L ✲ Figure 2: Simply-Supported Beam With Distributed Load. for a ≤ x ≤ b. A single nth-order ODE is equivalent to a system of n first-order ODEs. This equivalence can be seen by defining a new set of unknowns y1 , y2 , . . . , yn such that y1 = y, y2 = y ′ , y3 = y ′′ , . . . , yn = y (n−1) . For example, consider the third-order IVP y ′′′ = xy ′ + ex y + x2 + 1, x ≥ 0 y(0) = 1, y ′ (0) = 0, y ′′ (0) = −1. To obtain an equivalent first-order system, define y1 = y, y2 = y ′ , y3 = y ′′ to obtain  ′  y1 = y2 y ′ = y3  2′ y3 = xy2 + ex y1 + x2 + 1 (1.5) (1.6) with initial conditions y1 (0) = 1, y2 (0) = 0, y3 (0) = −1. 1.1 Euler’s Method This method is the simplest of the numerical methods for solving initial value problems. Consider the IVP  ′ y (x) = f (x, y), x ≥ a (1.7) y(a) = η. To effect a numerical solution, we discretize the x-axis: a = x 0 < x1 < x 2 < · · · , where, for uniform spacing, xi − xi−1 = h, (1.8) and h is considered small. With this discretization, we can approximate the derivative y ′ (x) with the forward finite difference y ′ (x) ≈ y(x + h) − y(x) . h (1.9) If we let yk represent the numerical approximation to y(xk ), then y ′ (xk ) ≈ yk+1 − yk . h 2 (1.10) Thus, a numerical (difference) approximation to the ODE, Eq. 1.7, is or  yk+1 − yk = f (xk , yk ), k = 0, 1, 2, . . . h (1.11) yk+1 = yk + hf (xk , yk ), k = 0, 1, 2, . . . y0 = η. (1.12) This recursive algorithm is called Euler’s method. 1.2 Truncation Error for Euler’s Method There are two types of error that arise in numerical methods: truncation error (which arises primarily from a discretization process) and rounding error (which arises from the finiteness of number representations in the computer). Refining a mesh to reduce the truncation error often causes the rounding error to increase. To estimate the truncation error for Euler’s method, we first recall Taylor’s theorem with remainder, which states that a function f (x) can be expanded in a series about the point x = c: f (x) = f (c)+f ′ (c)(x−c)+ f (n) (c) f (n+1) (ξ) f ′′ (c) (x−c)2 +· · ·+ (x−c)n + (x−c)n+1 , (1.13) 2! n! (n + 1)! where ξ is between x and c. The last term in Eq. 1.13 is referred to as the remainder term. Note also that Eq. 1.13 is an equality, not an approximation. In Eq. 1.13, let x = xk+1 and c = xk , in which case 1 y(xk+1 ) = y(xk ) + hy ′ (xk ) + h2 y ′′ (ξk ), 2 (1.14) where xk ≤ ξk ≤ xk+1 . Since y satisfies the ODE, Eq. 1.7, y ′ (xk ) = f (xk , y(xk )), (1.15) where y(xk ) is the actual solution at xk . Hence, 1 y(xk+1 ) = y(xk ) + hf (xk , y(xk )) + h2 y ′′ (ξk ). 2 (1.16) Like Eq. 1.13, this equation is an equality, not an approximation. By comparing this last equation to Euler’s approximation, Eq. 1.12, it is clear that Euler’s method is obtained by omitting the remainder term 21 h2 y ′′ (ξk ) in the Taylor expansion of y(xk+1 ) at the point xk . The omitted term accounts for the truncation error in Euler’s method at each step. This error is a local error, since the error occurs at each step regardless of the error in the previous step. The accumulation of local errors is referred to as the global error, which is the more important error but much more difficult to compute. Most algorithms for solving ODEs are derived by expanding the solution function in a Taylor series and then omitting certain terms. 3 1.3 Runge-Kutta Methods Euler’s method is a first-order method, since it was obtained by omitting terms in the Taylor series expansion containing powers of h greater than one. To derive a second-order method, we again use Taylor’s theorem with remainder to obtain 1 1 y(xk+1 ) = y(xk ) + hy ′ (xk ) + h2 y ′′ (xk ) + h3 y ′′′ (ξk ) 2 6 (1.17) for some ξk such that xk ≤ ξk ≤ xk+1 . Since, from the ODE (Eq. 1.7), y ′ (xk ) = f (xk , y(xk )), (1.18) we can approximate y ′′ (x) = df (x, y(x)) f (x + h, y(x + h)) − f (x, y(x)) = + O(h) dx h (1.19) where we use the “big O” notation O(h) to represent terms of order h as h → 0. [For example, 2h3 = O(h3 ), 3h2 + 5h4 = O(h2 ), h2 O(h) = O(h3 ), and −287h4 e−h = O(h4 ).] From these last two equations, Eq. 1.17 can then be written as h y(xk+1 ) = y(xk ) + hf (xk , y(xk )) + [f (xk+1 , y(xk+1 )) − f (xk , y(xk ))] + O(h3 ), 2 (1.20) which leads (after combining terms) to the difference equation 1 yk+1 = yk + h[f (xk , yk ) + f (xk+1 , yk+1 )]. 2 (1.21) This formula is a second-order approximation to the original differential equation y ′ (x) = f (x, y) (Eq. 1.7), but it is an inconvenient approximation, since yk+1 appears on both sides of the formula. (Such a formula is called an implicit method, since yk+1 is defined implicitly. An explicit method would have yk+1 appear only on the left-hand side.) To obtain instead an explicit formula, we use the approximation yk+1 = yk + hf (xk , yk ) (1.22) to obtain 1 yk+1 = yk + h[f (xk , yk ) + f (xk+1 , yk + hf (xk , yk ))]. (1.23) 2 This formula is the Runge-Kutta formula of second order. Other higher-order formulas can be derived similarly. For example, a fourth-order formula turns out to be popular in applications. We illustrate the implementation of the second-order Runge-Kutta formula, Eq. 1.23, with the following algorithm. We first make the following three definitions: ak = f (xk , yk ), (1.24) bk = yk + hak , (1.25) 4 ck = f (xk+1 , bk ), (1.26) in which case 1 yk+1 = yk + h(ak + ck ). (1.27) 2 The calculations can then be performed conveniently with the following spreadsheet: k 0 1 2 .. . 1.4 xk yk ak bk ck yk+1 Systems of Equations The methods just derived can be extended directly to systems of equations. Consider the initial value problem involving two equations:  ′  y (x) = f (x, y(x), z(x)) (1.28) z ′ (x) = g(x, y(x), z(x))  y(a) = η, z(a) = ξ. We recall from Eq. 1.12 that, for one equation, Euler’s method uses the recursive formula yk+1 = yk + hf (xk , yk ). This formula is directly extendible to two equations as   yk+1 = yk + hf (xk , yk , zk ) z = zk + hg(xk , yk , zk )  k+1 y0 = η, z0 = ξ. (1.29) (1.30) We recall from Eq. 1.23 that, for one equation, the second-order Runge-Kutta method uses the recursive formula 1 yk+1 = yk + h[f (xk , yk ) + f (xk+1 , yk + hf (xk , yk ))]. 2 For two equations, this formula becomes  1   yk+1 = yk + 2 h[f (xk , yk , zk ) + f (xk+1 , yk + hf (xk , yk , zk ), zk + hg(xk , yk , zk ))] zk+1 = zk + 21 h[g(xk , yk , zk ) + g(xk+1 , yk + hf (xk , yk , zk ), zk + hg(xk , yk , zk ))]   y0 = η, z0 = ξ. 5 (1.31) (1.32) ✻ y(x) ✑ ✑ s ✑ central ✑ difference ❍❍ ❍❍ ❨ ❥✑ ✏ s✑ ✑ s ✑ ✏ ✑ ✏❍ ❨❍ ✏ forward difference backward difference ✲ x−h x x+h x Figure 3: Finite Difference Approximations to Derivatives. 1.5 Finite Differences Before addressing boundary value problems, we want to develop further the notion of finite difference approximation of derivatives. Consider a function y(x) for which we want to compute the derivative y ′ (x) at some point x. If we discretize the x-axis with uniform spacing h, we could approximate the derivative using the forward difference formula y ′ (x) ≈ y(x + h) − y(x) , h (1.33) which is the slope of the line to the right of x (Fig. 3). We could also approximate the derivative using the backward difference formula y ′ (x) ≈ y(x) − y(x − h) , h (1.34) which is the slope of the line to the left of x. Since, in general, there is no basis for choosing one of these approximations over the other, an intuitively more appealing approximation results from the average of these formulas:   1 y(x + h) − y(x) y(x) − y(x − h) ′ (1.35) y (x) ≈ + 2 h h or y(x + h) − y(x − h) . (1.36) 2h This formula, which is more accurate than either the forward or backward difference formulas, is the central finite difference approximation to the derivative. y ′ (x) ≈ 6 Similar approximations can be derived for second derivatives. Using forward differences, y ′ (x + h) − y ′ (x) y (x) ≈ h y(x + 2h) − y(x + h) y(x + h) − y(x) − ≈ h2 h2 y(x + 2h) − 2y(x + h) + y(x) = . h2 ′′ (1.37) This formula, which involves three points forward of x, is the forward difference approximation to the second derivative. Similarly, using backward differences, y ′ (x) − y ′ (x − h) h y(x) − y(x − h) y(x − h) − y(x − 2h) ≈ − h2 h2 y(x) − 2y(x − h) + y(x − 2h) = . h2 y ′′ (x) ≈ (1.38) This formula, which involves three points backward of x, is the backward difference approximation to the second derivative. The central finite difference approximation to the second derivative uses instead the three points which bracket x: y ′′ (x) ≈ y(x + h) − 2y(x) + y(x − h) . h2 (1.39) This last result can also be obtained by using forward differences for the second derivative followed by backward differences for the first derivatives, or vice versa. The central difference formula for second derivatives can alternatively be derived using Taylor series expansions: h3 h2 ′′ y (x) + y ′′′ (x) + O(h4 ). 2 6 (1.40) h2 ′′ h3 ′′′ y(x − h) = y(x) − hy (x) + y (x) − y (x) + O(h4 ). 2 6 (1.41) y(x + h) = y(x) + hy ′ (x) + Similarly, by replacing h by −h, ′ The addition of these two equations yields y(x + h) + y(x − h) = 2y(x) + h2 y ′′ (x) + O(h4 ) (1.42) or y(x + h) − 2y(x) + y(x − h) + O(h2 ), 2 h which, because of the error term, shows that the formula is second-order accurate. y ′′ (x) = 7 (1.43) 1.6 Boundary Value Problems The techniques for initial value problems (IVPs) are, in general, not directly applicable to boundary value problems (BVPs). Consider the BVP  ′′ y (x) = f (x, y, y ′ ), a ≤ x ≤ b (1.44) y(a) = η1 , y(b) = η2 . This equation could be nonlinear, depending on f . The methods used for IVPs started at one end (x = a) and computed the solution step by step for increasing x. For a BVP, not enough information is given at either endpoint to allow a step-by-step solution. Consider first a special case of Eq. 1.44 for which the right-hand side depends only on x and y:  ′′ y (x) = f (x, y), a ≤ x ≤ b (1.45) y(a) = η1 , y(b) = η2 . Subdivide the interval (a, b) into n equal subintervals: h= b−a , n (1.46) in which case xk = a + kh, k = 0, 1, 2, . . . , n. (1.47) Let yk denote the numerical approximation to the exact solution at xk . That is, yk ≈ y(xk ). (1.48) Then, if we use a central difference approximation to the second derivative in Eq. 1.45, the ODE can be approximated by y(xk−1 ) − 2y(xk ) + y(xk+1 ) ≈ f (xk , yk ), h2 which suggests the difference equation yk−1 − 2yk + yk+1 = h2 f (xk , yk ), k = 1, 2, 3, . . . , n − 1. (1.49) (1.50) Since this system of equations has n − 1 equations in n + 1 unknowns, the two boundary conditions are needed to obtain a nonsingular system: y 0 = η 1 , y n = η2 . (1.51) The resulting system is thus −2y1 + y2 y1 − 2y2 + y3 y2 − 2y3 + y3 − y4 2y4 + y5 = = = = .. . −η1 + h2 f (x1 , y1 ) h2 f (x2 , y2 ) h2 f (x3 , y3 ) h2 f (x4 , y4 ) (1.52) yn−2 − 2yn−1 = −η2 + h2 f (xn−1 , yn−1 ), which is a tridiagonal system of n − 1 equations in n − 1 unknowns. This system is linear or nonlinear, depending on f . 8 1.6.1 Example Consider  y ′′ = −y(x), 0 ≤ x ≤ π/2 y(0) = 1, y(π/2) = 0. (1.53) In Eq. 1.45, f (x, y) = −y, η1 = 1, η2 = 0. Thus, the right-hand side of the ith equation in Eq. 1.52 has −h2 yi , which can be moved to the left-hand side to yield the system −(2 − h2 )y1 + y2 y1 − (2 − h2 )y2 + y3 2 y2 − (2 − h )y3 + y4 = −1 = 0 = 0 .. . yn−2 − (2 − h2 )yn−1 = (1.54) 0. We first solve this tridiagonal system of simultaneous equations with n = 8 (i.e., h = π/16), and compare with the exact solution y(x) = cos x: k 0 1 2 3 4 5 6 7 8 xk 0 0.1963495 0.3926991 0.5890486 0.7853982 0.9817477 1.1780972 1.3744468 1.5707963 yk 1 0.9812186 0.9246082 0.8323512 0.7080045 0.5563620 0.3832699 0.1954016 0 Exact y(xk ) 1 0.9807853 0.9238795 0.8314696 0.7071068 0.5555702 0.3826834 0.1950903 0 Absolute Error 0 0.0004334 0.0007287 0.0008816 0.0008977 0.0007917 0.0005865 0.0003113 0 % Error 0 0.0441845 0.0788715 0.1060315 0.1269565 0.1425085 0.1532604 0.1595767 0 Absolute Error 0 0.0000172 0.0000290 0.0000351 0.0000357 0.0000315 0.0000233 0.0000124 0 % Error 0 0.0017571 0.0031363 0.0042161 0.0050479 0.0056660 0.0060933 0.0063443 0 We then solve this system with n = 40 (h = π/80): k 0 5 10 15 20 25 30 35 40 xk 0 0.1963495 0.3926991 0.5890486 0.7853982 0.9817477 1.1780972 1.3744468 1.5707963 yk 1 0.9808025 0.9239085 0.8315047 0.7071425 0.5556017 0.3827068 0.1951027 0 Exact y(xk ) 1 0.9807853 0.9238795 0.8314696 0.7071068 0.5555702 0.3826834 0.1950903 0 Notice that a mesh refinement by a factor of 5 has reduced the error by a factor of about 25. This behavior is typical of a numerical method which is second-order accurate. 9 1.6.2 Solving Tridiagonal Systems Tridiagonal systems are particularly easy (and fast) to solve using Gaussian elimination. It is convenient to solve such systems using the following notation:  d 1 x1 + u1 x2 = b1    l 2 x1 + d 2 x2 + u2 x 3 = b2    l 3 x2 + d 3 x 3 + u3 x 4 = b3 (1.55)  ..   .    ln xn−1 + dn xn = bn , where di , ui , and li are, respectively, the diagonal, upper, and lower matrix entries in Row i. All coefficients can now be stored in three one-dimensional arrays, D(·), U(·), and L(·), instead of a full two-dimensional array A(I,J). The solution algorithm (reduction to upper triangular form by Gaussian elimination followed by back-solving) can now be summarized as follows: 1. For k = 1, 2, . . . , n − 1: (a) m = −lk+1 /dk [m = multiplier needed to annihilate term below] (b) dk+1 = dk+1 + muk (c) bk+1 = bk+1 + mbk 2. xn = bn /dn [k = pivot row] [new diagonal entry in next row] [new rhs in next row] [start of back-solve] 3. For k = n − 1, n − 2, . . . , 1: [back-solve loop] (a) xk = (bk − uk xk+1 )/dk Tridiagonal systems arise in a variety of applications, including the Crank-Nicolson finite difference method for solving parabolic partial differential equations. 1.7 Shooting Methods Shooting methods provide a way to convert a boundary value problem to a trial-and-error initial value problem. It is useful to have additional ways to solve BVPs, particularly if the equations are nonlinear. Consider the following two-point BVP:  ′′ ′  y = f (x, y, y ), a ≤ x ≤ b (1.56) y(a) = A  y(b) = B. To solve this problem using the shooting method, we compute solutions of the IVP  ′′ ′  y = f (x, y, y ), x ≥ a (1.57) y(a) = A  ′ y (a) = M 10 y ✻ M2 r M1 A a ✻ B ❄ ✲ b Figure 4: The Shooting Method. x for various values of M (the slope at the left end of the domain) until two solutions, one with y(b) < B and the other with y(b) > B, have been found (Fig. 4). The initial slope M can then be interpolated until a solution is found (i.e., y(b) = B). 2 Partial Differential Equations A partial differential equation (PDE) is an equation that involves an unknown function (the dependent variable) and some of its partial derivatives with respect to two or more independent variables. An nth-order equation has the highest order derivative of order n. 2.1 Classical Equations of Mathematical Physics 1. Laplace’s equation (the potential equation) ∇2 φ = 0 (2.1) In Cartesian coordinates, the vector operator del is defined as ∇ = ex ∂ ∂ ∂ + ey + ez . ∂x ∂y ∂z (2.2) ∇2 is referred to as the Laplacian operator and given by ∂2 ∂2 ∂2 ∇ =∇·∇= + + . ∂x2 ∂y 2 ∂z 2 2 (2.3) Thus, Laplace’s equation in Cartesian coordinates is ∂ 2φ ∂ 2φ ∂ 2φ + + 2 = 0. ∂x2 ∂y 2 ∂z In cylindrical coordinates, the Laplacian is   ∂φ 1 ∂ 2φ ∂ 2φ 1 ∂ 2 r + 2 2 + 2. ∇ φ= r ∂r ∂r r ∂θ ∂z (2.4) (2.5) Laplace’s equation arises in incompressible fluid flow (in which case φ is the velocity potential), gravitational potential problems, electrostatics, magnetostatics, steady-state 11 heat conduction with no sources (in which case φ is the temperature), and torsion of bars in elasticity (in which case φ(x, y) is the warping function). Functions which satisfy Laplace’s equation are referred to as harmonic functions. 2. Poisson’s equation ∇2 φ + g = 0 (2.6) This equation arises in steady-state heat conduction with distributed sources (φ = temperature) and torsion of bars in elasticity (in which case φ(x, y) is the stress function). 3. wave equation 1 φ̈ c2 In this equation, dots denote time derivatives, e.g., ∇2 φ = φ̈ = ∂ 2φ , ∂t2 (2.7) (2.8) and c is the speed of propagation. The wave equation arises in several physical situations: (a) transverse vibrations of a string For this one-dimensional problem, φ = φ(x, t) is the transverse displacement of the string, and s T c= , (2.9) ρA where T is the string tension, ρ is the density of the string material, and A is the cross-sectional area of the string. The denominator ρA is mass per unit length. (b) longitudinal vibrations of a bar For this one-dimensional problem, φ = φ(x, t) represents the longitudinal displacement, and s E , (2.10) c= ρ where E and ρ are, respectively, the modulus of elasticity and density of the bar material. (c) transverse vibrations of a membrane For this two-dimensional problem, φ = φ(x, y, t) is the transverse displacement of the membrane (e.g., drum head), and r T c= , (2.11) m where T is the tension per unit length, m is the mass per unit area (i.e., m = ρt), and t is the membrane thickness. 12 (d) acoustics For this three-dimensional problem, φ = φ(x, y, z, t) is the fluid pressure or velocity potential, and c is the speed of sound, where s B c= , (2.12) ρ where B = ρc2 is the fluid bulk modulus, and ρ is the density. 4. Helmholtz equation (reduced wave equation) ∇2 φ + k 2 φ = 0 (2.13) The Helmholtz equation is the time-harmonic form of the wave equation, in which interest is restricted to functions which vary sinusoidally in time. To obtain the Helmholtz equation, we substitute φ(x, y, z, t) = φ0 (x, y, z) cos ωt (2.14) into the wave equation, Eq. 2.7, to obtain ∇2 φ0 cos ωt = − ω2 φ0 cos ωt. c2 (2.15) If we define the wave number k = ω/c, this equation becomes ∇2 φ0 + k 2 φ0 = 0. (2.16) With the understanding that the unknown depends only on the spacial variables, the subscript is unnecessary, and we obtain the Helmholtz equation, Eq. 2.13. This equation arises in steady-state (time-harmonic) situations involving the wave equation, e.g., steady-state acoustics. 5. heat equation ∇ · (k∇φ) + Q = ρcφ̇ (2.17) In this equation, φ represents the temperature T , k is the thermal conductivity, Q is the internal heat generation per unit volume per unit time, ρ is the material density, and c is the material specific heat (the heat required per unit mass to raise the temperature by one degree). The thermal conductivity k is defined by Fourier’s law of heat conduction: dT , (2.18) dx where q̂x is the rate of heat conduction (energy per unit time) with typical units J/s or BTU/hr, and A is the area through which the heat flows. Alternatively, Fourier’s law is written dT qx = −k , (2.19) dx where qx is energy per unit time per unit area (with typical units J/(s·m2 ). There are several special cases of the heat equation of interest: q̂x = −kA 13 (a) homogeneous material (k = constant): k∇2 φ + Q = ρcφ̇ (2.20) (b) homogeneous material, steady-state (time-independent): ∇2 φ = − Q k (Poisson’s equation) (2.21) (c) homogeneous material, steady-state, no sources (Q = 0): ∇2 φ = 0 2.2 (Laplace’s equation) (2.22) Classification of Partial Differential Equations Of the classical PDEs summarized in the preceding section, some involve time, and some don’t, so presumably their solutions would exhibit fundamental differences. Of those that involve time (wave and heat equations), the order of the time derivative is different, so the fundamental character of their solutions may also differ. Both these speculations turn out to be true. Consider the general, second-order, linear partial differential equation in two variables Auxx + Buxy + Cuyy + Dux + Euy + F u = G, (2.23) where the coefficients are functions of the independent variables x and y (i.e., A = A(x, y), B = B(x, y), etc.), and we have used subscripts to denote partial derivatives, e.g., uxx = ∂ 2u . ∂x2 (2.24) The quantity B 2 − 4AC is referred to as the discriminant of the equation. The behavior of the solution of Eq. 2.23 depends on the sign of the discriminant according to the following table: B 2 − 4AC <0 =0 >0 Equation Type Elliptic Parabolic Hyperbolic Typical Physics Described Steady-state phenomena Heat flow and diffusion processes Vibrating systems and wave motion The names elliptic, parabolic, and hyperbolic arise from the analogy with the conic sections in analytic geometry. Given these definitions, we can classify the common equations of mathematical physics already encountered as follows: Eq. Number Name Laplace Eq. 2.1 Eq. 2.6 Poisson Wave Eq. 2.7 Eq. 2.13 Helmholtz Heat Eq. 2.17 Eq. in Two Variables uxx + uyy = 0 uxx + uyy = −g uxx − uyy /c2 = 0 uxx + uyy + k 2 u = 0 kuxx − ρcuy = −Q 14 A, B, C A = C = 1, B = 0 A = C = 1, B = 0 A = 1, C = −1/c2 , B = 0 A = C = 1, B = 0 A = k, B = C = 0 Type Elliptic Elliptic Hyperbolic Elliptic Parabolic In the wave and heat equations in the above table, y represents the time variable. The behavior of the solutions of equations of different types differs. Elliptic equations characterize static (time-independent) situations, and the other two types of equations characterize timedependent situations. 2.3 Transformation to Nondimensional Form It is often convenient, when solving equations, to transform the equation to a nondimensional form. Consider, for example, the one-dimensional wave equation ∂ 2u 1 ∂ 2u = , ∂x2 c2 ∂t2 (2.25) where u is displacement (of dimension length), t is time, and c is the speed of propagation (of dimension length/time). Let L represent some characteristic length associated with the problem. We define the nondimensional variables x̄ = x , L ū = u , L t̄ = ct , L (2.26) in which case the derivatives in Eq. 2.25 become ∂(Lū) ∂ ū ∂u = = , ∂x ∂(Lx̄) ∂ x̄       ∂ ∂ ū ∂ ∂ ū dx̄ ∂ ∂u 1 ∂ 2 ū ∂ 2u = = = , = ∂x2 ∂x ∂x ∂x ∂ x̄ ∂ x̄ ∂ x̄ dx L ∂ x̄2 ∂u ∂(Lū) ∂ ū = =c , ∂t ∂(Lt̄/c) ∂ t̄       ∂ ∂ ū ∂ ∂ ū dt̄ ∂ 2u ∂ ∂u ∂ 2 ū c c2 ∂ 2 ū = c = c = . = c = ∂t2 ∂t ∂t ∂t ∂ t̄ ∂ t̄ ∂ t̄ dt ∂ t̄2 L L ∂ t̄2 Thus, from Eq. 2.25, 1 c2 ∂ 2 ū 1 ∂ 2 ū = · L ∂ x̄2 c2 L ∂ t̄2 (2.27) (2.28) (2.29) (2.30) (2.31) or ∂ 2 ū ∂ 2 ū = . (2.32) ∂ x̄2 ∂ t̄2 This is the nondimensional wave equation. This last equation can also be obtained more easily by direct substitution of Eq. 2.26 into Eq. 2.25 and factoring out the constants: or 1 ∂ 2 (Lū) ∂ 2 (Lū) = ∂(Lx̄)2 c2 ∂(Lt̄/c)2 (2.33) 1 c2 L ∂ 2 ū L ∂ 2 ū = . L2 ∂ x̄2 c2 L2 ∂ t̄2 (2.34) 15 3 3.1 Finite Difference Solution of Partial Differential Equations Parabolic Equations Consider the boundary-initial value problem (BIVP)  1    uxx = ut , u = u(x, t), 0 < x < 1, t > 0 c u(0, t) = u(1, t) = 0 (boundary conditions)    u(x, 0) = f (x) (initial condition), (3.1) where c is a constant. This problem represents transient heat conduction in a rod with the ends held at zero temperature and an initial temperature profile f (x). To solve this problem numerically, we discretize x and t such that xi = ih, i = 0, 1, 2, . . . tj = jk, j = 0, 1, 2, . . . . 3.1.1 (3.2) Explicit Finite Difference Method Let ui,j be the numerical approximation to u(xi , tj ). We approximate ut with the forward finite difference ui,j+1 − ui,j (3.3) ut ≈ k and uxx with the central finite difference uxx ≈ ui+1,j − 2ui,j + ui−1,j . h2 (3.4) The finite difference approximation to the PDE is then ui,j+1 − ui,j ui+1,j − 2ui,j + ui−1,j = . h2 ck Define the parameter r as r= ck c ∆t = , 2 h (∆x)2 (3.5) (3.6) in which case Eq. 3.5 becomes ui,j+1 = rui−1,j + (1 − 2r)ui,j + rui+1,j . (3.7) The domain of the problem and the mesh are illustrated in Fig. 5. Eq. 3.7 is a recursive relationship giving u in a given row (time) in terms of three consecutive values of u in the row below (one time step earlier). This equation is referred to as an explicit formula since one unknown value can be found directly in terms of several other known values. The recursive relationship can also be sketched with the stencil shown in Fig. 6. For example, for r = 1/10, we have the stencil shown in Fig. 7. That is, for r = 1/10, the solution (temperature) at 16 t ✻ i, j +1 s u=0 ✻ s i−1, j k ❄ ✛ s s i, j u=0 i+1, j h ✲ ✲ 0 u = f (x) 1 x Figure 5: Mesh for 1-D Heat Equation. ✓ ✒ r ✓ ✒ ✏✓ 1 ✑✒ ✏ ✑ 1 − 2r ✏✓ ✑✒ r ✏ ✑ Figure 6: Heat Equation Stencil for Explicit Finite Difference Algorithm. ✛✘ 1 ✚✙ ✬✩ ✬✩ ✬✩ 1 10 8 10 1 10 ✫✪ ✫✪ ✫✪ Figure 7: Heat Equation Stencil for r = 1/10. 1 r= 2 ✛✘ ✛✘ 1 1 ✚✙ r=1 ✚✙ ✬✩ ✬✩ ✬✩ ✬✩ ✬✩ ✬✩ 1 2 0 1 2 1 −1 1 ✫✪ ✫✪ ✫✪ ✫✪ ✫✪ ✫✪ Figure 8: Heat Equation Stencils for r = 1/2 and r = 1. 17 u ✻ 0.5 t = 0.048 0.4 0.3 0.2 0.1 0.0 ❝ ❝ ❝ ❝ ❝ ❝ 0.2 ❝ ❝ ❝ ❝ ❝ ❝ t = 0.096 ❝ ❝ ❝ t = 0.192 ❝ ❝ ❝ 0.4 0.6 ❝ ❝ ❝ r = 0.48 ❝ ❝ ❝ 0.8 ❝ ❝ Exact F.D. ❝ ❝ ❝ ✲ x 1.0 Figure 9: Explicit Finite Difference Solution With r = 0.48. the new point depends on the three points at the previous time step with a 1-8-1 weighting. Notice that, if r = 1/2, the solution at the new point is independent of the closest point, as illustrated in Fig. 8. For r > 1/2 (e.g., r = 1), the new point depends negatively on the closest point (Fig. 8), which is counter-intuitive. It can be shown that, for a stable solution, 0 < r ≤ 1/2. An unstable solution is one for which small errors grow rather than decay as the solution evolves. The instability which occurs for r > 1/2 can be illustrated with the following example. Consider the boundary-initial value problem (in nondimensional form) uxx = ut , 0 < x < 1, u = u(x, t) u(0, t) = u(1, t) = 0  2x, 0 ≤ x ≤ 1/2 u(x, 0) = f (x) = 2(1 − x), 1/2 ≤ x ≤ 1. (3.8) The physical problem is to compute the temperature history u(x, t) for a bar with a prescribed initial temperature distribution f (x), no internal heat sources, and zero temperature prescribed at both ends. We solve this problem using the explicit finite difference algorithm with h = ∆x = 0.1 and k = ∆t = rh2 = r(∆x)2 for two different values of r: r = 0.48 and r = 0.52. The two numerical solutions (Figs. 9 and 10) are compared with the analytic solution ∞ X nπ 8 2 sin (sin nπx) e−(nπ) t , (3.9) u(x, t) = 2 (nπ) 2 n=1 which can be obtained by the technique of separation of variables. The instability for r > 1/2 can be clearly seen in Fig. 10. Thus, a disadvantage of this explicit method is that a small time step ∆t must be used to maintain stability. This disadvantage will be removed with the Crank-Nicolson algorithm. 3.1.2 Crank-Nicolson Implicit Method The Crank-Nicolson method is a stable algorithm which allows a larger time step than could be used in the explicit method. In fact, Crank-Nicolson’s stability does not depend on the parameter r. 18 u ✻ 0.5 t = 0.052 ❝ ❝ 0.4 0.3 0.2 0.1 0.0 ❝ ❝ ❝ ❝ ❝ ❝ 0.2 ❝ ❝ ❝ ❝ t = 0.104 ❝ ❝ ❝t 0.4 ❝ ❝ ❝ ❝ = 0.208❝ ❝ 0.6 r = 0.52 ❝ ❝ ❝ 0.8 ❝ ❝ ❝ ❝ ❝ Exact F.D. ✲ 1.0 x Figure 10: Explicit Finite Difference Solution With r = 0.52. The basis for the Crank-Nicolson algorithm is writing the finite difference equation at a mid-level in time: i, j + 21 . The finite difference x derivative at j + 12 is computed as the average of the two central difference time derivatives at j and j + 1. Consider again the PDE of Eq. 3.1:  1    uxx = ut , u = u(x, t) c (3.10) u(0, t) = u(1, t) = 0 (boundary conditions)    u(x, 0) = f (x) (initial condition). The PDE is approximated numerically by   1 ui+1,j+1 − 2ui,j+1 + ui−1,j+1 ui+1,j − 2ui,j + ui−1,j ui,j+1 − ui,j + , = 2 2 2 h h ck (3.11) where the right-hand side is a central difference approximation to the time derivative at the middle point j + 21 . We again define the parameter r as r= ck c ∆t = , 2 h (∆x)2 (3.12) and rearrange Eq. 3.11 with all j + 1 terms on the left-hand side: − rui−1,j+1 + 2(1 + r)ui,j+1 − rui+1,j+1 = rui−1,j + 2(1 − r)ui,j + rui+1,j . (3.13) This formula is called the Crank-Nicolson algorithm. Fig. 11 shows the points involved in the Crank-Nicolson scheme. If we start at the bottom row (j = 0) and move up, the right-hand side values of Eq. 3.13 are known, and the left-hand side values of that equation are unknown. To get the process started, let j = 0, and write the C-N equation for each i = 1, 2, . . . , N to obtain N simultaneous equations in N unknowns, where N is the number of interior mesh points on the row. (The boundary points, with known values, are excluded.) This system of equations is a tridiagonal system, since each equation has three consecutive nonzeros centered around the diagonal. To advance in time, 19 t ✻ s s u=0 ✻ i−1, j k ❄ ✛ s i−1,j +1 ❞i, j +1 s i, j s i+1,j +1 s u=0 i+1, j h ✲ ✲ 0 u = f (x) 1 x Figure 11: Mesh for Crank-Nicolson. ✓ ✏✓ ✒ ✑✒ ✒ −r ✓ r ✑✒ 2(1 + r) ✏✓ 2(1 − r) ✏✓ ✏ ✑✒ ✑ ✑✒ −r ✏✓ r ✑ ✏ Figure 12: Stencil for Crank-Nicolson Algorithm. we then increment j to j = 1, and solve a new system of equations. An approach which requires the solution of simultaneous equations is called an implicit algorithm. A sketch of the C-N stencil is shown in Fig. 12. Note that the coefficient matrix of the C-N system of equations does not change from step to step. Thus, one could compute and save the LU factors of the coefficient matrix, and merely do the forward-backward substitution (FBS) at each new time step, thus speeding up the calculation. This speedup would be particularly significant in higher dimensional problems, where the coefficient matrix is no longer tridiagonal. It can be shown that the C-N algorithm is stable for any r, although better accuracy results from a smaller r. A smaller r corresponds to a smaller time step size (for a fixed spacial mesh). C-N also gives better accuracy than the explicit approach for the same r. 3.1.3 Derivative Boundary Conditions Consider the boundary-initial value problem  1    uxx = ut , u = u(x, t) c u(0, t) = 0, ux (1, t) = g(t) (boundary conditions)    u(x, 0) = f (x) (initial condition). (3.14) The only difference between this problem and the one considered earlier in Eq. 3.1 is the right-hand side boundary condition, which now involves a derivative (a Neumann boundary condition). Assume a mesh labeling as shown in Fig. 13. We introduce extra “phantom” points to the right of the boundary (outside the domain). Consider boundary Point 25, for example, 20 15 25 14 24 13 23 ♣ ♣ ♣ ❝35 ♣ ♣ ♣ ❝34 ♣ ♣ ♣ ❝33 Figure 13: Treatment of Derivative Boundary Conditions. and assume we use an explicit finite difference algorithm. Since u25 is not known, we must write the finite difference equation for u25 : u25 = ru14 + (1 − 2r)u24 + ru34 . (3.15) On the other hand, a central finite difference approximation to the x derivative at Point 24 is u34 − u14 = g24 . (3.16) 2h The phantom variable u34 can then be eliminated from the last two equations to yield a new equation for the boundary point u25 : u25 = 2ru14 + (1 − 2r)u24 + 2rhg24 . 3.2 3.2.1 (3.17) Hyperbolic Equations The d’Alembert Solution of the Wave Equation Before addressing the finite difference solution of hyperbolic equations, we review some background material on such equations. The time-dependent transverse response of an infinitely long string satisfies the onedimensional wave equation with nonzero initial displacement and velocity specified:  2 1 ∂ 2u ∂ u   = , −∞ < x < ∞, t > 0,    ∂x2 c2 ∂t2 ∂u(x, 0) (3.18) = g(x) u(x, 0) = f (x),   ∂t    lim u(x, t) = 0, x→±∞ where x is distance along the string, t is time, u(x, t) is the transverse displacement, f (x) is the initial displacement, g(x) is the initial velocity, and the constant c is given by s T , (3.19) c= ρA where T is the tension (force) in the string, ρ is the density of the string material, and A is the cross-sectional area of the string. Note that c has the dimension of velocity. This 21 equation assumes that all motion is vertical and that the displacement u and its slope ∂u/∂x are both small. It can be shown by direct substitution into Eq. 3.18 that the solution of this system is Z 1 1 x+ct u(x, t) = [f (x − ct) + f (x + ct)] + g(τ ) dτ. (3.20) 2 2c x−ct The differentiation of the integral in Eq. 3.20 is effected with the aid of Leibnitz’s rule: d dx Z B(x) h(x, t) dt = A(x) Z B A ∂h(x, t) dB dA dt + h(x, B) − h(x, A) . ∂x dx dx (3.21) Eq. 3.20 is known as the d’Alembert solution of the one-dimensional wave equation. For the special case g(x) = 0 (zero initial velocity), the d’Alembert solution simplifies to 1 u(x, t) = [f (x − ct) + f (x + ct)], 2 (3.22) which may be interpreted as two waves, each equal to f (x)/2, which travel at speed c to the right and left, respectively. For example, the argument x − ct remains constant if, as t increases, x also increases at speed c. Thus, the wave f (x − ct) moves to the right (increasing x) with speed c without change of shape. Similarly, the wave f (x + ct) moves to the left (decreasing x) with speed c without change of shape. The two waves [each equal to half the initial shape f (x)] travel in opposite directions from each other at speed c. If f (x) is nonzero only for a small domain, then, after both waves have passed the region of initial disturbance, the string returns to its rest position. For example, let f (x), the initial displacement, be given by f (x) =  −|x| + b, |x| ≤ b, 0, |x| ≥ b, (3.23) which is a triangular pulse of width 2b and height b (Fig. 14). For t > 0, half this pulse travels in opposite directions from the origin. For t > b/c, where c is the wave speed, the two half-pulses have completely separated, and the neighborhood of the origin has returned to rest. For the special case f (x) = 0 (zero initial displacement), the d’Alembert solution simplifies to Z 1 1 x+ct g(τ ) dτ = [G(x + ct) − G(x − ct)] , (3.24) u(x, t) = 2c x−ct 2 where 1 G(x) = c Z x g(τ ) dτ. (3.25) −∞ Thus, similar to the initial displacement special case, this solution, Eq. 3.24, may be interpreted as the combination (difference, in this case) of two identical functions G(x)/2, one moving left and one moving right, each with speed c. 22 ✻ ❅ (a) t = 0 ✟ −3 −2 (b) t = ✟✟ −1 ✟✟❍❍ −2 (d) t = −1 ✲ x/b ✲ x/b ✲ x/b ❍ ✲ x/b 2 3 0 1 2 3 ✟✟❍❍ ❍❍ ❍❍✟✟ 0 1 2 3 ✻ 2b c ✟ ✟✟ 1 ✻ b c −3 0 f (x)/2 ✟❍❍ ✟✟ ❍ ❍ ✟ ✟ ❍❍ ✟ ❍ ✟ −2 (c) t = −3 −1 ✟❍❅✠ ❅ ✙✟ ✟ ❍ ❅ ❍ ✻ b 2c −3 ✟ f (x) ✟❍ ❍❍ ❍ −2 −1 ✟✟ 0 ✟❍ ❍❍ ✟ 1 2 3 Figure 14: Propagation of Initial Displacement. For example, let the initial velocity g(x) be given by  M c, |x| < b, g(x) = 0, |x| > b, (3.26) which is a rectangular pulse of width 2b and height M c, where the constant M is the dimensionless Mach number, and c is the wave speed (Fig. 15). The travelling wave G(x) is given by Eq. 3.25 as  x ≤ −b,  0, G(x) = (3.27) M (x + b), −b ≤ x ≤ b,  2M b, x ≥ b. That is, half this wave travels in opposite directions at speed c. Even though g(x) is nonzero only near the origin, the travelling wave G(x) is constant and nonzero for x > b. Thus, as time advances, the center section of the string reaches a state of rest, but not in its original position (Fig. 16). From the preceding discussion, it is clear that disturbances travel with speed c. For an observer at some fixed location, initial displacements occurring elsewhere pass by after a finite time has elapsed, and then the string returns to rest in its original position. Nonzero initial velocity disturbances also travel at speed c, but, once having reached some location, will continue to influence the solution from then on. Thus, the domain of influence of the data at x = x0 , say, on the solution consists of all points closer than ct (in either direction) to x0 , the location of the disturbance (Fig. 17). 23 ✻ g(x) (a) −3 −2 −1 0 1 2 3 2 3 ✲ x/b ✲ x/b ✻ G(x) (b) ✘✘✘ ✘✘✘ ✘ ✘ −3 −2 −1 0 1 Figure 15: Initial Velocity Function. (a) t = 0 (displacement=0) −4 (b) t = −3 b 2c −4 (c) t = −4 (d) t = −4 −2 ✘ ✘ ❳ ❳ −1 G(x)/2 moves left ✻ ✘ ✘ ❳ ✻ ✘ ✘ ✘ ❳❳ ✘ ❳❳ ✘✘ ❳ ❳ ❳ 1❳ −1 −3 −2 −3 ✻❳ ✘ ✘❳ ✘ ❳❳❳ ✘ ✘ ❳❳❳ ✘✘✘ ❳ ❳ ❳ ❳ 2 −2 −1 b c 2b c ✘✘✘ ✘✘✘ −3 −2 ✘✘ ✻ ❳❳❳ ❳ ❳ −1 1 2 3 3 x/b ✲ x/b ✲ x/b ✲ x/b 4 4 ❳❳❳ ❳❳ ❳ ❳ 3 Figure 16: Propagation of Initial Velocity. 24 ✲ 2 3 4 ❳ 1 −G(x)/2 moves right 4 ✻ t ✻ t ❅ ❅ ❅ −1/c❅ ❅ ❅r 1/c slope=1/c ✲ (x0 , 0) r x r ❅ ❅ (x, t) (x − ct, 0) ❅−1/c ❅ ❅r ✲ (x + ct, 0) x (b) Domain of Dependence of Solution on Data. (a) Domain of Influence of Data on Solution. Figure 17: Domains of Influence and Dependence. Conversely, the domain of dependence of the solution on the initial data consists of all points within a distance ct of the solution point. That is, the solution at (x, t) depends on the initial data for all locations in the range (x − ct, x + ct), which are the limits of integration in the d’Alembert solution, Eq. 3.20. 3.2.2 Finite Differences From the preceding discussion of the d’Alembert solution, we see that hyperbolic equations involve wave motion. If the initial data are discontinuous (as, for example, in shocks), the most accurate and the most convenient approach for solving the equations is probably the method of characteristics. On the other hand, problems without discontinuities can probably be solved most conveniently using finite difference and finite element techniques. Here we consider finite differences. Consider the boundary-initial value problem (BIVP)  1    uxx = 2 utt , u = u(x, t), 0 < x < a, t > 0 c (3.28) u(0, t) = u(a, t) = 0 (boundary conditions)    u(x, 0) = f (x), u (x, 0) = g(x) (initial conditions). t This problem represents the transient (time-dependent) vibrations of a string fixed at the two ends with both initial displacement f (x) and initial velocity g(x) specified. A central finite difference approximation to the PDE, Eq. 3.28, yields ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1 = . 2 h c2 k 2 We define the parameter r= ck c ∆t = , h ∆x (3.29) (3.30) and solve for ui,j+1 : ui,j+1 = r2 ui−1,j + 2(1 − r2 )ui,j + r2 ui+1,j − ui,j−1 . (3.31) Fig. 18 shows the mesh points involved in this recursive scheme. If we know the solution for 25 t ✻ i, j +1 s u=0 ✻ s i−1, j k ❄ ✛ s i, j s u=0 i+1, j s i, j −1 h ✲ 0 ✲ u = f (x), ut = g(x) a x Figure 18: Mesh for Explicit Solution of Wave Equation. ✓ ✏ ✒ ✑ ✓ ✏ ✓ ✏✓ ✒ r2 ✑✒ 1 2(1 − r2 ) ✒ −1 ✑ ✏✓ ✏ ✑✒ r2 ✑ Figure 19: Stencil for Explicit Solution of Wave Equation. all time values up to the jth time step, we can compute the solution ui,j+1 at Step j + 1 in terms of known quantities. Thus, this algorithm is an explicit algorithm. The corresponding stencil is shown in Fig. 19. It can be shown that this finite difference algorithm is stable if r ≤ 1 and unstable if r > 1. This stability condition is know as the Courant, Friedrichs, and Lewy (CFL) condition or simply the Courant condition. It can be further shown that a theoretically correct solution is obtained when r = 1, and that the accuracy of the solution decreases as the value of r decreases farther and farther below the value of 1. Thus, the time step size ∆t should be chosen, if possible, so that r = 1. If that ∆t is inconvenient for a particular calculation, ∆t should be selected as large as possible without exceeding the stability limit of r = 1. An intuitive rationale behind the stability requirement r ≤ 1 can also be made using Fig. 20. If r > 1, the numerical domain of dependence (NDD) would be smaller than the actual domain of dependence (ADD) for the PDE, since NDD spreads by one mesh point at each level for earlier time. However, if NDD < ADD, the numerical solution would be independent of data outside NDD but inside ADD. That is, the numerical solution would ignore necessary information. Thus, to insure a stable solution, r ≤ 1. 3.2.3 Starting Procedure for Explicit Algorithm We note that the explicit finite difference scheme just described for the wave equation requires the numerical solution at two consecutive time steps to step forward to the next time. Thus, at t = 0, a special procedure is needed to advance from j = 0 to j = 1. 26 t ✻ ❝ Slope = −1/c ❅❝ ❅❝ ✠ ★ ❅❝ ★ ❅❝ ★ ❅❝ ❝ ★ ❅ ❝ ★ ★ ❅ ❝ ★ ❅ ❝ ★ ★ ❅ ❝❝ ✛ Numerical Domain of Dependence ✲ Slope = 1/c ✻ k ❄ ✛ ✲ h s ❝ ★❅ ✛ ❩★ ⑦ ★ ★ Actual Domain of Dependence ✲ x ✲ Figure 20: Domains of Dependence for r > 1. The explicit finite difference algorithm is given in Eq. 3.31: ui,j+1 = r2 ui−1,j + 2(1 − r2 )ui,j + r2 ui+1,j − ui,j−1 . (3.32) To compute the solution at the end of the first time step, let j = 0: ui,1 = r2 ui−1,0 + 2(1 − r2 )ui,0 + r2 ui+1,0 − ui,−1 , (3.33) where the right-hand side is known (from the initial condition) except for ui,−1 . However, we can write a central difference approximation to the first time derivative at t = 0: ui,1 − ui,−1 = gi 2k (3.34) ui,−1 = ui,1 − 2kgi , (3.35) or where gi is the initial velocity g(x) evaluated at the ith point, i.e., gi = g(xi ). If we substitute this last result into Eq. 3.33 (to eliminate ui,−1 ), we obtain ui,1 = r2 ui−1,0 + 2(1 − r2 )ui,0 + r2 ui+1,0 − ui,1 + 2kgi (3.36) or 1 1 (3.37) ui,1 = r2 ui−1,0 + (1 − r2 )ui,0 + r2 ui+1,0 + kgi . 2 2 This is the difference equation used for the first row. Thus, to implement the explicit finite difference algorithm, we use Eq. 3.37 for the first time step and Eq. 3.31 for all subsequent time steps. 3.2.4 Nonreflecting Boundaries In some applications, it is of interest to model domains that are large enough to be considered infinite in extent. In a finite difference representation of the domain, an infinite boundary has 27 s s i, j s s ❝ ❝ s ❝ ❝ Figure 21: Finite Difference Mesh at Nonreflecting Boundary. to be truncated at some sufficiently large distance. At such a boundary, a suitable boundary condition must be imposed to ensure that outgoing waves are not reflected. Consider a vibrating string which extends to infinity for large x. We truncate the computational domain at some finite x. Let the initial velocity be zero. The d’Alembert solution, Eq. 3.20, of the one-dimensional wave equation c2 uxx = utt can be written in the form u(x, t) = F1 (x − ct) + F2 (x + ct), (3.38) where F1 represents a wave advancing at speed c toward the boundary, and F2 represents the returning wave, which should not exist if the boundary is nonreflecting. With F2 = 0, we differentiate u with respect to x and t to obtain ∂u ∂u = F1′ , = −cF1′ , ∂x ∂t where the prime denotes the derivative with respect to the argument. Thus, (3.39) ∂u 1 ∂u =− . (3.40) ∂x c ∂t This is the one-dimensional nonreflecting boundary condition. Note that the x direction is normal to the boundary. The boundary condition, Eq. 3.40, must be imposed to inhibit reflections from the truncated boundary. This condition is exact in 1-D (i.e., plane waves) and approximate in higher dimensions, where the nonreflecting condition is written ∂u ∂u + = 0, (3.41) ∂n ∂t where n is the outward unit normal to the boundary. The nonreflecting boundary condition, Eq. 3.40, can be approximated in the finite difference method with central differences expressed in terms of the phantom point outside the boundary. For example, at the typical point (i, j) on the nonreflecting boundary in Fig. 21, the general recursive formula is given by Eq. 3.31: c ui,j+1 = r2 ui−1,j + 2(1 − r2 )ui,j + r2 ui+1,j − ui,j−1 . (3.42) The central difference approximation to the nonreflecting boundary condition, Eq. 3.40, is c ui,j+1 − ui,j−1 ui+1,j − ui−1,j =− 2h 2k 28 (3.43) or r (ui+1,j − ui−1,j ) = − (ui,j+1 − ui,j−1 ) . (3.44) The substitution of Eq. 3.44 into Eq. 3.42 (to eliminate the phantom point) yields (1 + r)ui,j+1 = 2r2 ui−1,j + 2(1 − r2 )ui,j − (1 − r)ui,j−1 . (3.45) For the first time step (j = 0), the last term in this relation is evaluated using the central difference approximation to the initial velocity, Eq. 3.35. Note also that, for r = 1, Eq. 3.45 takes a particularly simple (and perhaps unexpected) form: ui,j+1 = ui−1,j . (3.46) To illustrate the perfect wave absorption that occurs in a one-dimensional finite difference model, consider an infinitely-long vibrating string with a nonzero initial displacement and zero initial velocity. The initial displacement is a triangular-shaped pulse in the middle of the string, similar to Fig. 14. According to the d’Alembert solution, half the pulse should propagate at speed c to the left and right and be absorbed into the boundaries. We solve the problem with the explicit central finite difference approach with r = ck/h = 1. With r = 1, the finite difference formulas 3.31 and 3.37 simplify to ui,j+1 = ui−1,j + ui+1,j − ui,j−1 (3.47) ui,1 = (ui−1,0 + ui+1,0 )/2. (3.48) and On the right and left nonreflecting boundaries, Eq. 3.46 implies un,j+1 = un−1,j , u0,j+1 = u1,j , (3.49) where the mesh points in the x direction are labeled 0 to n. The finite difference calculation for this problem results in the following spreadsheet: t x=0 x=1 x=2 x=3 x=4 x=5 x=6 x=7 x=8 x=9 x=10 --------------------------------------------------------------------0 0.000 0.000 0.000 1.000 2.000 3.000 2.000 1.000 0.000 0.000 0.000 1 0.000 0.000 0.500 1.000 2.000 2.000 2.000 1.000 0.500 0.000 0.000 2 0.000 0.500 1.000 1.500 1.000 1.000 1.000 1.500 1.000 0.500 0.000 3 0.500 1.000 1.500 1.000 0.500 0.000 0.500 1.000 1.500 1.000 0.500 4 1.000 1.500 1.000 0.500 0.000 0.000 0.000 0.500 1.000 1.500 1.000 5 1.500 1.000 0.500 0.000 0.000 0.000 0.000 0.000 0.500 1.000 1.500 6 1.000 0.500 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.500 1.000 7 0.500 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.500 8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 29 ρcA ❅ ❅ ❅ Figure 22: Finite Length Simulation of an Infinite Bar. Notice that the triangular wave is absorbed without any reflection from the two boundaries. For steady-state wave motion, the solution u(x, t) is time-harmonic, i.e., u = u0 eiωt , (3.50) and the nonreflecting boundary condition, Eq. 3.40, becomes ∂u0 iω = − u0 ∂x c (3.51) or, in general, iω ∂u0 = − u0 = −iku0 , (3.52) ∂n c where n is the outward unit normal at the nonreflecting boundary, and k = ω/c is the wave number. This condition is exact in 1-D (i.e., plane waves) and approximate in higher dimensions. The nonreflecting boundary condition can be interpreted physically as a damper (dashpot). Consider, for example, a bar undergoing longitudinal vibration and terminated on the right end with the nonreflecting boundary condition, Eq. 3.40 (Fig. 22). The internal longitudinal force F in the bar is given by F = Aσ = AEε = AE ∂u , ∂x (3.53) where A is the cross-sectional area of the bar, σ is the stress, E is the Young’s modulus of the bar material, and u is displacement. Thus, from Eq. 3.40, the nonreflecting boundary condition is equivalent to applying an end force given by F =− AE v, c (3.54) where v = ∂u/∂t is the velocity. Since, from Eq. 2.10, E = ρc2 , (3.55) F = −(ρcA)v, (3.56) Eq. 3.54 becomes which is a force proportional to velocity. A mechanical device which applies a force proportional to velocity is a dashpot. The minus sign in this equation means that the force opposes the direction of motion, as required to be physically realizable. The dashpot constant is ρcA. Thus, the application of this dashpot to the end of a finite length bar simulates exactly a bar of infinite length (Fig. 22). Since, in acoustics, the ratio of pressure (force/area) to velocity is referred to as impedance, we see that the characteristic impedance of an acoustic medium is ρc. 30 b T = g1 (y) ✻ y T = f2 (x) ∇2 T = 0 T = f1 (x) T = g2 (y) a ✲ x Figure 23: Laplace’s Equation on Rectangular Domain. ✻ j+1 j j−1 s✠ Point (i, j) ✲ i−1 i i+1 Figure 24: Finite Difference Grid on Rectangular Domain. 3.3 Elliptic Equations Consider Laplace’s equation on the two-dimensional rectangular domain shown in Fig. 23:  2 ∇ T (x, y) = 0 (0 < x < a, 0 < y < b), (3.57) T (0, y) = g1 (y), T (a, y) = g2 (y), T (x, 0) = f1 (x), T (x, b) = f2 (x). This problem corresponds physically to two-dimensional steady-state heat conduction over a rectangular plate for which temperature is specified on the boundary. We attempt an approximate solution by introducing a uniform rectangular grid over the domain, and let the point (i, j) denote the point having the ith value of x and the jth value of y (Fig. 24). Then, using central finite difference approximations to the second derivatives (Fig. 25), Ti−1,j − 2Ti,j + Ti+1,j ∂ 2T ≈ , (3.58) 2 ∂x h2 Ti,j−1 − 2Ti,j + Ti,j+1 ∂ 2T ≈ . (3.59) ∂y 2 h2 The finite difference approximation to Laplace’s equation thus becomes Ti−1,j − 2Ti,j + Ti+1,j Ti,j−1 − 2Ti,j + Ti,j+1 + =0 h2 h2 (3.60) 4Ti,j − (Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1 ) = 0. (3.61) or 31 (i − 1, j + 1) (i, j + 1) ✻ h s (i − 1, j) (i, j) ✛ (i − 1, j − 1) (i + 1, j + 1) (i + 1, j) ❄ h ✲ (i, j − 1) (i + 1, j − 1) Figure 25: The Neighborhood of Point (i, j). 4 8 12 16 20 3 7 11 15 19 2 6 10 14 18 1 5 9 13 17 Figure 26: 20-Point Finite Difference Mesh. That is, for Laplace’s equation with the same uniform mesh in each direction, the solution at a typical point (i, j) is the average of the four neighboring points. For example, consider the mesh shown in Fig. 26. Although there are 20 mesh points, 14 are on the boundary, where the temperature is known. Thus, the resulting numerical problem has only six degrees of freedom (unknown variables). The application of Eq. 3.61 to each of the six interior points yields  4T6 − T7 − T10 = T2 + T5      −T6 + 4T7 − T11 = T3 + T8    −T6 + 4T10 − T11 − T14 = T9 (3.62)  −T7 − T10 + 4T11 − T15 = T12     −T10 + 4T14 − T15 = T13 + T18    −T11 − T14 + 4T15 = T16 + T19 , where all known quantities have been placed on the right-hand side. This linear system of six equations in six unknowns can be solved with standard equation solvers. Because the central difference operator is a 5-point operator, systems of equations of this type would have at most five nonzero terms in each equation, regardless of how large the mesh is. Thus, for large meshes, the system of equations is sparsely populated, so that sparse matrix solution techniques would be applicable. 32 b ✻ y T = f2 (x) ∂T = g(y) ∂x ∇2 T = 0 T = g1 (y) T = f1 (x) a ✲ x Figure 27: Laplace’s Equation With Dirichlet and Neumann B.C. Since the numbers assigned to the mesh points in Fig. 26 are merely labels for identification, the pattern of nonzero coefficients appearing on the left-hand side of Eq. 3.62 depends on the choice of mesh ordering. Some equation solvers based on Gaussian elimination operate more efficiently on systems of equations for which the nonzeros in the coefficient matrix are clustered near the main diagonal. Such a matrix system is called banded. Systems of this type can also be solved using an iterative procedure known as relaxation, which uses the following general algorithm: 1. Initialize the boundary points to their prescribed values, and initialize the interior points to zero or some other convenient value (e.g., the average of the boundary values). 2. Loop systematically through the interior mesh points, setting each interior point to the average of its four neighbors. 3. Continue this process until the solution converges to the desired accuracy. 3.3.1 Derivative Boundary Conditions The approach of the preceding section must be modified for Neumann boundary conditions, in which the normal derivative is specified. For example, consider again the problem of the last section but with a Neumann, rather than Dirichlet, boundary condition on the right side (Fig. 27):   ∇2 T (x, y) = 0 (0 < x < a, 0 < y < b), (3.63) ∂T (a, y)  T (0, y) = g1 (y), = g(y), T (x, 0) = f1 (x), T (x, b) = f2 (x). ∂x We extend the mesh to include additional points to the right of the boundary at x = a (Fig. 28). At a typical point on the boundary, a central difference approximation yields g18 = ∂T18 T22 − T14 ≈ ∂x 2h (3.64) or T22 = T14 + 2hg18 . (3.65) On the other hand, the equilibrium equation for Point 18 is 4T18 − (T14 + T22 + T17 + T19 ) = 0, 33 (3.66) ... 16 20 ... 15 19 ... 14 18 ... 13 17 ♣ ♣ ♣ ❜24 ♣ ♣ ♣ ❜23 ♣ ♣ ♣ ❜22 ♣ ♣ ♣ ❜21 Figure 28: Treatment of Neumann Boundary Conditions. ✲ k1 ✁✁❆❆✁✁❆❆✁✁ m2 ❡ ✲ u2 (t) k2 ❡ ✁✁❆❆✁✁❆❆✁✁ m3 ❡ u3 (t) ❡ ✲ f3 (t) Figure 29: 2-DOF Mass-Spring System. which, when combined with Eq. 3.65, yields 4T18 − 2T14 − T17 − T19 = 2hg18 , (3.67) which imposes the Neumann boundary condition on Point 18. 4 Direct Finite Element Analysis The finite element method is a numerical procedure for solving partial differential equations. The procedure is used in a variety of applications, including structural mechanics and dynamics, acoustics, heat transfer, fluid flow, electric and magnetic fields, and electromagnetics. Although the main theoretical bases for the finite element method are variational principles and the weighted residual method, it is useful to consider discrete systems first to gain some physical insight into some of the procedures. 4.1 Linear Mass-Spring Systems Consider the two-degree-of-freedom (DOF) system shown in Fig. 29. We let u2 and u3 denote the displacements from the equilibrium of the two masses m2 and m3 . The stiffnesses of the two springs are k1 and k2 . The dynamic equilibrium equations could be obtained from Newton’s second law (F=ma): m2 ü2 + k1 u2 − k2 (u3 − u2 ) = 0 m3 ü3 + k2 (u3 − u2 ) = f3 (t) 34 (4.1) ❝ ❝ ✁✁❆❆✁✁❆❆✁✁ 1 2 k Figure 30: A Single Spring Element. or m2 ü2 + (k1 + k2 )u2 − k2 u3 = 0 m3 ü3 − k2 u2 + k2 u3 = f3 (t). (4.2) Mü + Ku = F(t), (4.3)  (4.4) This system could be rewritten in matrix notation as where u= u2 u3  is the displacement vector (the vector of unknown displacements),   k1 + k2 −k2 K= −k2 k2 (4.5) is the system stiffness matrix, M=   m2 0 0 m3 (4.6) is the system mass matrix, and F=  0 f3  (4.7) is the force vector. This approach would be very tedious and error-prone for more complex systems involving many springs and masses. To develop instead a matrix approach, we first isolate one element, as shown in Fig. 30. The stiffness matrix Kel for this element satisfies Kel u = F, or  k11 k12 k21 k22  u1 u2  (4.8) =  f1 f2  . (4.9) By expanding this equation, we obtain k11 u1 + k12 u2 = f1 k21 u1 + k22 u2 = f2  . (4.10) From this equation, we observe that k11 can be defined as the force on DOF 1 corresponding to enforcing a unit displacement on DOF 1 and zero displacement on DOF 2: k11 = f1 u1 =1,u2 =0 35 = k. (4.11) ❝ 1 ✁✁❆❆✁✁❆❆✁✁ k1 ❝ ❝ ✁✁❆❆✁✁❆❆✁✁ 2 3 k2 Figure 31: 3-DOF Spring System. Similarly, k12 = f1 u1 =0,u2 =1 = −k, k21 = f2 u1 =1,u2 =0 = −k, k22 = f2 u1 =0,u2 =1 = k. (4.12) Thus, el K =  k −k −k k  =k  1 −1 −1 1  . (4.13) In general, for a larger system with many more DOF, K11 u1 + K12 u2 + K13 u3 + · · · + K1n un = F1 K21 u1 + K22 u2 + K23 u3 + · · · + K2n un = F2 K31 u1 + K32 u2 + K33 u3 + · · · + K3n un = F3 .. . Kn1 u1 + Kn2 u2 + Kn3 u3 + · · · + Knn un = Fn        , (4.14)       in which case we can interpret an individual element Kij in the stiffness matrix as the force at DOF i if uj = 1 and all other displacement components ui = 0: Kij = Fi 4.2 uj =1, others=0 . (4.15) Matrix Assembly We now return to the stiffness part of the original problem shown in Fig. 29. In the absence of the masses and constraints, this system is shown in Fig. 31. Since there are three points in this system, each with one DOF, the system is a 3-DOF system. The system stiffness matrix can be assembled for this system by adding the 3 × 3 stiffness matrices for each element:      k1 −k1 0 0 0 0 k1 −k1 0 K =  −k1 k2 −k2  =  −k1 k1 + k2 −k2  . k1 0  +  0 0 −k2 k2 0 −k2 k2 0 0 0  (4.16) The justification for this assembly procedure is that forces are additive. For example, k22 is the force at DOF 2 when u2 = 1 and u1 = u3 = 0; both elements which connect to DOF 2 contribute to k22 . This matrix corresponds to the unconstrained system. 4.3 Constraints The system in Fig. 29 has a constraint on DOF 1, as shown in Fig. 32. If we expand the 36 fixed r✠ ✁❆ ✁❆ ✁ ❆✁ ❆✁ 1 k1 ❝ 2 ❝ ✁✁❆❆✁✁❆❆✁✁ 3 k2 Figure 32: Spring System With Constraint. k1 1 ❝ ✁✁❆❆✁✁❆❆✁✁ 2❝ k2 ✁✁❆❆✁✁❆❆✁✁ 3 ✁✁❆❆✁✁❆❆✁✁ ❝ k3 ✁✁❆❆✁✁❆❆✁✁ ❝4 k4 Figure 33: 4-DOF Spring System. (unconstrained) matrix system Ku = F into  K11 u1 + K12 u2 + K13 u3 = F1  , K21 u1 + K22 u2 + K23 u3 = F2  K31 u1 + K32 u2 + K33 u3 = F3 (4.17) we see that Row i corresponds to the equilibrium equation for DOF i. Thus, if ui = 0, we do not need Row i of the matrix (although that equation can be saved to recover later the constraint force). Also, Column i of the matrix multiplies ui , so that, if ui = 0, we do not need Column i of the matrix. That is, if ui = 0 for some system, we can enforce that constraint by deleting Row i and Column i from the unconstrained matrix. Hence, with u1 = 0 in Fig. 32, we delete Row 1 and Column 1 in Eq. 4.16 to obtain the reduced matrix   k1 + k2 −k2 , (4.18) K= −k2 k2 which is the same matrix obtained previously in Eq. 4.5. Notice that, from Eqs. 4.16 and 4.18, det K3×3 = k1 [(k1 + k2 )k2 − k22 ] + k1 (−k1 k2 ) = 0, (4.19) det K2×2 = (k1 + k2 )k2 − k22 = k1 k2 6= 0. (4.20) whereas That is, the unconstrained matrix K3×3 is singular, but the constrained matrix K2×2 is nonsingular. Without constraints, K is singular (and the solution of the mechanical problem is not unique) because of the presence of rigid body modes. 4.4 Example and Summary Consider the 4-DOF spring system shown in Fig. 33. The unconstrained stiffness matrix for 37 this system is    K=  k1 + k4 −k1 −k4 0 −k1 k1 + k2 −k2 0 −k4 −k2 k2 + k3 + k4 −k3 0 0 −k3 k3 We summarize several properties of stiffness matrices:    .  (4.21) 1. K is symmetric. This property is a special case of the Betti reciprocal theorem in mechanics. 2. An off-diagonal term is zero unless the two points are common to the same element. Thus, K is sparse in general and usually banded. 3. K is singular without enough constraints to eliminate rigid body motion. For spring systems, that have only one DOF at each point, the sum of any matrix column or row is zero. This property is a consequence of equilibrium, since the matrix entries in Column i consist of all the grid point forces when ui = 1 and other DOF are fixed. The forces must sum to zero, since the object is in static equilibrium. We summarize the solution procedure for spring systems: 1. Generate the element stiffness matrices. 2. Assemble the system K and F. 3. Apply constraints. 4. Solve Ku = F for u. 5. Compute reactions, spring forces, and stresses. 4.5 Pin-Jointed Rod Element Consider the pin-jointed rod element (an axial member) shown in Fig. 34. From mechanics of materials, we recall that the change in displacement u for a rod of length L subjected to an axial force F is FL , (4.22) u= AE where A is the rod cross-sectional area, and E is the Young’s modulus for the rod material. Thus, the axial stiffness is F AE k= = . (4.23) u L The rod is therefore equivalent to a scalar spring with k = AE/L, and AE K = L el  1 −1 −1 1 38  . (4.24) ❛ ✒ F AE, L ✠ ❛ F Figure 34: Pin-Jointed Rod Element. r ❅ r ❅ ❅ ❅ r ❅ ❅ ❅ ❅ ❅ ❅ ❅ r ❅ ❅ ❅ ❅r ❅ ❅ ❅r ❅ ❅ ❅ ❅ ❅ ❅r ❄ Figure 35: Truss Structure Modeled With Pin-Jointed Rods. 39 4 ✻ ❛ ✲3 θ 2 ✻ ❛ ✲1 Figure 36: The Degrees of Freedom for a Pin-Jointed Rod Element in 2-D. ✚ ✚ ✚ ✚ ✚ r (x2 , y2 ) ✚ 1 cos θ ✚ ❆ θ (x1 , y1 ) ❆❯ ✚✚ ❘ r✚ ❅ F y ✻✚ ❃ 1 ✚✲ Fx r Figure 37: Computing 2-D Stiffness of Pin-Jointed Rod. However, matrix assembly for a truss structure (a structure made of pin-jointed rods) differs from that for a collection of springs, since the rod elements are not all colinear (e.g., Fig. 35). Thus the elements must be transformed to a common coordinate system before the element matrices can be assembled into the system stiffness matrix. A typical rod element in 2-D is shown in Fig. 36. To use this element in 2-D trusses requires expanding the 2 × 2 matrix Kel into a 4 × 4 matrix. The four DOF for the element are also shown in Fig. 36. To compute the entries in the first column of the 4 × 4 Kel , we set u1 = 1, u2 = u3 = u4 = 0, and compute the four grid point forces F1 , F2 , F3 , F4 , as shown in Fig. 37. For example, at Point 1, Fx = F cos θ = (k · 1 cos θ) cos θ = k cos2 θ = k11 = −k31 (4.25) Fy = F sin θ = (k · 1 cos θ) sin θ = k cos θ sin θ = k21 = −k41 , (4.26) where k = AE/L, and the forces at the opposite end of the rod (x2 , y2 ) are equal in magnitude and opposite in sign. These four values complete the first column of the matrix. Similarly we can find the rest of the matrix to obtain   c2 cs −c2 −cs  AE  s2 −cs −s2   cs el K = (4.27) ,  c2 cs  L  −c2 −cs −cs −s2 cs s2 40 10,000 lb ❆ ❆ s=c 3 4 ❆ ✻ 60◦ ❆❆❯ ❛ ✲3 ❅ 2❅❅❅❅ = − La ❍ s = −c = La ❅ ❅✟ ❥ ❍ ✙✟ ❅ ❅ ❅ ❅ ❅ ❅ 45◦ 6 2 ❅ ❅ ✻ ❅ ❅ ✻ ❅ ❅ ❛ ✲5 1❅ ❛ ✲1 ✛ 2a ✲ Figure 38: Pin-Jointed Frame Example. where y2 − y1 x2 − x1 , s = sin θ = . (4.28) L L Note that the angle θ is not of interest; only c and s are needed. Later, in the discussion of matrix transformations, we will derive a more convenient way to obtain this matrix. c = cos θ = 4.6 Pin-Jointed Frame Example We illustrate the matrix assembly and solution procedures for truss structures with the simple frame shown in Fig. 38. For this frame, we assume E = 30 × 106 psi, L=10 in, and A=1 in2 . Before the application of constraints, the stiffness matrix is   1 −1 −1 1 0 0  −1 1 1 −1 0 0      1 1 + 1 −1 + 1 −1 −1   −1 (4.29) K = k0  ,  1 −1 −1 + 1 1 + 1 −1 −1     0 0 −1 −1 1 1  0 0 −1 −1 1 1 where AE  a 2 k0 = = 1.5 × 106 lb/in, L L 1 a =√ . L 2 After constraints, the system of equations is      u3 5000 2 0 √ , = k0 −5000 3 0 2 u4 (4.30) (4.31) whose solution is u3 = 1.67 × 10−3 in, u4 = −2.89 × 10−3 in. 41 (4.32) ❄ ❄ ❄ ❄ ❄ ❄ ❄ r r r r r r ❅ ❄ ❄ ❄ ❄ r r r ❅ Figure 39: Example With Reactions and Loads at Same DOF. The reactions can be recovered as R1 = k0 (−u3 + u4 ) = −6840 lb, R5 = k0 (−u3 − u4 ) = 1830 lb, 4.7 R2 = k0 (u3 − u4 ) = 6840 lb R6 = k0 (−u3 − u4 ) = 1830 lb. (4.33) Boundary Conditions by Matrix Partitioning We recall that, to enforce ui = 0, we could delete Row i and Column i from the stiffness matrix. By using matrix partitioning, we can treat nonzero constraints and recover the forces of constraint. Consider the finite element matrix system Ku = F, where some DOF are specified, and some are free. We partition the unknown displacement vector into   uf , (4.34) u= us where uf and us denote, respectively, the free and constrained DOF. A partitioning of the matrix system then yields      Ff uf Kf f Kf s , (4.35) = Fs us Ksf Kss where us is prescribed. From the first partition, Kf f uf + Kf s us = Ff (4.36) Kf f uf = Ff − Kf s us , (4.37) or which can be solved for the unknown uf . The second partition then yields the forces at the constrained DOF: Fs = Ksf uf + Kss us , (4.38) where uf is now known, and us is prescribed. The grid point forces Fs at the constrained DOF consist of both the reactions of constraint and the applied loads, if any, at those DOF. For example, consider the beam shown in Fig. 39. When the applied load is distributed to the grid points, the loads at the two end grid points would include both reactions and a portion of the applied load. Thus, Fs , which includes all loads at the constrained DOF, can be written as Fs = Ps + R, 42 (4.39) i r ✲ ui = U prescribed F i = k0 U ✲ k0 ✁✁❆❆✁✁❆❆✁✁ ❅ ❅ ❅ (b) (a) Figure 40: Large Spring Approach to Constraints. where Ps is the applied load, and R is the vector of reactions. The reactions can then be recovered as R = Ksf uf + Kss us − Ps . (4.40) The disadvantage of using the partitioning approach to constraints is that the writing of computer programs is made more complicated, since one needs the general capability to partition a matrix into smaller matrices. However, for general purpose software, such a capability is essential. 4.8 Alternative Approach to Constraints Consider a structure (Fig. 40) for which we want to prescribe a value for the ith DOF: ui = U . An alternative approach to enforce this constraint is to attach a large spring k0 between DOF i and ground (a fixed point) and to apply a force Fi to DOF i equal to k0 U . This spring should be many orders of magnitude larger than other typical values in the stiffness matrix. A large number placed on the matrix diagonal will have no adverse effects on the conditioning of the matrix. For example, if we prescribe DOF 3 in a 4-DOF system, the modified system becomes         k11 k12 k13 k14  u F 1  1             k k k k u F  21 22 23 24  2 2 = . (4.41)    u3  F 3 + k0 U   k31 k32 k33 + k0 k34            k41 k42 k43 k44  u4   F4 For large k0 , the third equation is k0 u3 ≈ k0 U or u3 ≈ U, (4.42) as required. The large spring approach can be used for any system of equations for which one wants to enforce a constraint on a particular unknown. The main advantage of this approach is that computer coding is easier, since matrix sizes do not have to change. We can summarize the algorithm for applying the large-spring approach for applying constraints to the linear system Ku = F as follows: 1. Choose the large spring k0 to be, say, 10,000 times the largest diagonal entry in the unconstrained K matrix. 43 u u ✻1 u2 r ✲ EI, L u4 x 1 ✻3 r 2 Figure 41: DOF for Beam in Flexure (2-D). F1 ✻ M1 ✻ w1 = 1 r F2 ✻ M2 ❄ r Figure 42: The Beam Problem Associated With Column 1. 2. For each DOF i which is to be constrained (zero or not), add k0 to the diagonal entry Kii , and add k0 U to the corresponding right-hand side term Fi , where U is the desired constraint value for DOF i. 4.9 Beams in Flexure Like springs and pin-jointed rods, beam elements also have stiffness matrices which can be computed exactly. In two dimensions, we designate the four DOF associated with flexure as shown in Fig. 41. The stiffness matrix would therefore be of the form       F1  w1  k11 k12 k13 k14               M1 θ k k k k  21 22 23 24  1 = , (4.43)   F2  w2   k31 k32 k33 k34           M   k41 k42 k43 k44  θ2  2 where wi and θi are, respectively, the transverse displacement and rotation at the ith point. Rotations are expressed in radians. To compute the first column of K, set w1 = 1, θ1 = w2 = θ2 = 0, (4.44) and compute the four forces. Thus, for an Euler beam, we solve the beam differential equation EIy ′′ = M (x), (4.45) subject to the boundary conditions, Eq. 4.44, as shown in Fig. 42. For Column 2, we solve the beam equation subject to the boundary conditions θ1 = 1, w1 = w2 = θ2 = 0, (4.46) as shown in Fig. 43. If we then combine the resulting flexural stiffnesses with the axial 44 F1 ✻ M1 F2 ✻ M2 r r Figure 43: The Beam Problem Associated With Column 2. u3 u2 ✻ r ✲ u1 u5 ✻ u6 A, E, I, L r 1 ✲ u4 2 Figure 44: DOF for 2-D Beam Element. stiffnesses previously derived for the axial member, matrix for the Euler beam:      AE AE     F 0 0 − 1      L  L          6EI 12EI     0 F2    0    L3 L2           4EI 6EI     F3  0  0 2 L L =  AE   AE   −    F 0 0   4    L L           6EI 12EI     − 2 0 − 3 F5     0   L L           2EI 6EI    F6  0 0 2 L L we obtain the two-dimensional stiffness 0 0 12EI L3 6EI − 2 L 6EI L2 2EI L − 0 12EI L3 6EI − 2 L 0 6EI L2 4EI L −                                            u1         u2         u3  , (4.47)   u4         u5         u6  where the six DOF are shown in Fig. 44. For this element, note that there is no coupling between axial and transverse behavior. Transverse shear, which was ignored, can be added. The three-dimensional counterpart to this matrix would have six DOF at each grid point: three translations and three rotations (ux , uy , uz , Rx , Ry , Rz ). Thus, in 3-D, K is a 12 × 12 matrix. In addition, for bending in two different planes, there would have to be two moments of inertia I1 and I2 , in addition to a torsional constant J and (possibly) a product of inertia I12 . For transverse shear, “area factors” would also be needed to reflect the fact that two different beams with the same cross-sectional area, but different cross-sectional shapes, would have different shear stiffnesses. 4.10 Direct Approach to Continuum Problems Stiffness matrices for springs, rods, and beams can be derived exactly. For most problems of interest in engineering, exact stiffness matrices cannot be derived. Consider the 2-D problem of computing the displacements and stresses in a thin plate with applied forces and 45 Constraint ✕ ✁ ✁ ✁ ✁✕ ✁ ✕ ✁ ✍✂ ✍✂ Load ✁ ✁ ✂ ✂ ✂ ✂ Figure 45: Plate With Constraints and Loads. v3 ✻ (x3 , y3 ) v1 ✄ ✄ ✄ ✄ ✄ ✄ ✄❅ ✄3 ❅ ✲ u3 ❅ ❅ ❅ ❅ v2 ❅ ✻ ❅ ✲ u2 2✘ ✄ ✻ ✘❅ ✘ ✘ ✄ (x2 , y2 ) ✘✘ ✄1 ✘✘✘ ✘ ✘ ✲ u1 ✄ (x1 , y1 ) ✘ Figure 46: DOF for the Linear Triangular Membrane Element. constraints (Fig. 45). This figure also shows the domain modeled with several triangular elements. A typical element is shown in Fig. 46. With three grid points and two DOF at each point, this is a 6-DOF element. The displacement and force vectors for the element are     u1  F1x                   F v 1  1y             u2 F2x u= , F= . (4.48)   v2  F2y                u3  F3x              v3 F3y Since we cannot determine exactly a stiffness matrix that relates these two vectors, we approximate the displacement field over the element as  u(x, y) = a1 + a2 x + a3 y (4.49) v(x, y) = a4 + a5 x + a6 y, where u and v are the x and y components of displacement, respectively. Note that the number of undetermined coefficients equals the number of DOF (6) in the element. We 46 choose polynomials for convenience in the subsequent mathematics. The linear approximation in Eq. 4.49 is analogous to the piecewise linear approximation used in trapezoidal rule integration. At the vertices, the displacements in Eq. 4.49 must match the grid point values, e.g., u1 = a1 + a2 x1 + a3 y1 . We can write five more equations of this type to obtain    1 x1 y1 0 0 0 u1           v1    0 0 0 1 x1 y1       u2  1 x2 y2 0 0 0 =   0 0 0 1 x2 y2 v2           1 x3 y3 0 0 0 u3        v3 0 0 0 1 x3 y3 or u = γa. (4.50)                    a1 a2 a3 a4 a5 a6                  (4.51) (4.52) Since the x and y directions uncouple in Eq. 4.51, we can write this last equation in uncoupled form:      1 x1 y1  a1   u1  =  1 x2 y2  a2 (4.53) u  2    u3 a3 1 x3 y3      1 x 1 y 1  a4   v1  . (4.54) =  1 x 2 y 2  a5 v    2  a6 v3 1 x3 y3 We now observe that  1 x1 y1 det  1 x2 y2  = 2A 6= 0, 1 x3 y3  (4.55) since |A| is (from geometry) the area of the triangle. A is positive for counterclockwise ordering and negative for clockwise ordering. Thus, unless the triangle is degenerate, we conclude that γ is invertible, and      x2 y 3 − x3 y2 x3 y1 − x 1 y 3 x 1 y 2 − x2 y 1  u1   a1  1  . (4.56) = y2 − y3 y3 − y1 y 1 − y2  u2 a2    2A  x3 − x2 x1 − x3 x2 − x1 u3 a3 Similarly,      x2 y3 − x 3 y2 x 3 y 1 − x1 y 3 x1 y 2 − x2 y 1  v1   a4  1  = . a5 y2 − y3 y3 − y1 y 1 − y 2  v2  2A    x3 − x2 x1 − x3 x2 − x1 a6 v3 47 (4.57) Thus, we can write a = γ −1 u, (4.58) and treat the element’s unknowns as being either the physical displacements u or the nonphysical coefficients a. The strain components of interest in 2-D are     u,x  εxx    ε= = . (4.59) εyy v,y     γxy u,y + v,x From Eq. 4.49, we evaluate the strains for this element as   a1             a 2      0 1 0 0 0 0 a 2     a 3   = 0 0 0 0 0 1 ε= = Ba = Bγ −1 u = Cu, a6     a 4    0 0 1 0 1 0  a3 + a5    a  5       a6 (4.60) where C = Bγ −1 . (4.61) Eq. 4.60 is a formula to compute element strains given the grid point displacements. Note that, for this element, the strains are independent of position in the element. Thus, this element is referred to as the Constant Strain Triangle (CST). From generalized Hooke’s law, each stress component is a linear combination of all the strain components: σ = Dε = DBγ −1 u = DCu, (4.62) where, for example, for an isotropic material in plane stress,  1 ν 0 E  , D= ν 1 0 1 − ν2 0 0 (1 − ν)/2  (4.63) where E is Young’s modulus, and ν is Poisson’s ratio. Eq. 4.62 is a formula to compute element stresses given the grid point displacements. For this element, the stresses are constant (independent of position) in the element. We now derive the element stiffness matrix using the Principle of Virtual Work. Consider an element in static equilibrium with a set of applied loads F and displacements u. According to the Principle of Virtual Work, the work done by the applied loads acting through the displacements is equal to the increment in stored strain energy during an arbitrary virtual (imaginary) displacement δu: Z δεT σ dV. δuT F = V 48 (4.64) We then substitute Eqs. 4.60 and 4.62 into this equation to obtain  Z Z T T T T C DC dV u, δu F = (C δu) (DCu) dV = δu (4.65) V V where δu and u can be removed from the integral, since they are grid point quantities independent of position. Since δu is arbitrary, it follows that  Z T C DC dV u. (4.66) F= V Therefore, the stiffness matrix is given by Z Z T C DC dV = K= V Bγ −1 V
T
D Bγ −1 dV, (4.67) where the integral is over the volume of the element. Note that, since the material property matrix D is symmetric, K is symmetric. The substitution of the constant B, γ, and D matrices into this equation yields the stiffness matrix for the Constant Strain Triangle: Et K= × 4|A|(1 − ν 2 )   2 y23 +λx232 µx32 y23 y23 y31 +λx13 x32 νx13 y23 +λx32 y31 y12 y23 +λx21 x32 νx21 y23 +λx32 y12     2 2  x32 +λy23 νx32 y31 +λx13 y23 x32 x13 +λy31 y23 νx32 y12 +λx21 y23 x21 x32 +λy12 y23      2 2  y31 +λx13 µx13 y31 y12 y31 +λx21 x13 νx21 y31 +λx13 y12   ,   2 2  x13 +λy31 νx13 y12 +λx21 y31 x13 x21 +λy12 y31      2 2   Sym y12 +λx21 µx21 y12   2 x221 +λy12 (4.68) where xij = xi − xj , yij = yi − yj , λ= 1−ν , 2 µ= 1+ν , 2 (4.69) and t is the element thickness. 5 Change of Basis On many occasions in engineering applications, the need arises to transform vectors and matrices from one coordinate system to another. For example, in the finite element method, it is frequently more convenient to derive element matrices in a local element coordinate system and then transform those matrices to a global system (Fig. 47). Transformations are also needed to transform from other orthogonal coordinate systems (e.g., cylindrical or spherical) to Cartesian coordinates (Fig. 48). 49 y y ✻ ▼❇ ✻ ȳ ❇ ❇ r✏✏✏ ✶ x̄ ✏ r✏✏ ✏ ✏✏ ▼❇ axial member ✲ r❳❳ ✶ x̄ ✏ ❳❳r✏✏ ✏ ✏ ❇ ✡ ✏✏ ❇✏ r✡ ✏ triangular membrane ȳ ✡ ✲ x x Figure 47: Element Coordinate Systems in the Finite Element Method. y ✻ eθ ❇▼ ❇ r ❇ ✶ er ✏✏ ❇✏ q ✏ θ ✲ x Figure 48: Basis Vectors in Polar Coordinate System. Let the vector v be given by v = v1 e1 + v2 e2 + v 3 e3 = 3 X vi ei , (5.1) i=1 where ei are the basis vectors, and vi are the components of v. With the summation convention, we can omit the summation sign and write v = vi ei , (5.2) where, if a subscript appears exactly twice, a summation is implied over the range. An orthonormal basis is defined as a basis whose basis vectors are mutually orthogonal unit vectors (i.e., vectors of unit length). If ei is an orthonormal basis, ei · ej = δij =  1, i = j, 0, i = 6 j, (5.3) where δij is the Kronecker delta. Since bases are not unique, we can express v in two different orthonormal bases: v= 3 X v i ei = i=1 3 X v̄i ēi , (5.4) i=1 where vi are the components of v in the unbarred coordinate system, and v̄i are the components in the barred system (Fig. 49). If we take the dot product of both sides of Eq. 5.4 with ej , we obtain 3 3 X X v i ei · ej = (5.5) v̄i ēi · ej , i=1 i=1 50 x̄2 x2 ✻ ❇▼ ❇ v ✒ ❇ ❇ ❇ ❇ ✏ ❇✏✏ ✏✏ ✏✏ ✏✏ ✶ ✏✏ θ x̄1 ✲ x1 Figure 49: Change of Basis. where ei · ej = δij , and we define the 3 × 3 matrix R as Rij = ēi · ej . (5.6) Thus, from Eq. 5.5, vj = 3 X Rij v̄i = i=1 Since the matrix product 3 X T Rji v̄i . (5.7) i=1 C = AB (5.8) can be written using subscript notation as Cij = 3 X Aik Bkj , (5.9) k=1 Eq. 5.7 is equivalent to the matrix product v = RT v̄. (5.10) Similarly, if we take the dot product of Eq. 5.4 with ēj , we obtain 3 X i=1 vi ei · ēj = 3 X i=1 v̄i ēi · ēj , (5.11) where ēi · ēj = δij , and ei · ēj = Rji . Thus, v̄j = 3 X Rji vi or v̄ = Rv or v = R−1 v̄. (5.12) i=1 A comparison of Eqs. 5.10 and 5.12 yields R−1 = RT or RRT = I or 3 X k=1 51 Rik Rjk = δij , (5.13) where I is the identity matrix (Iij = δij ):  1 0 0 I =  0 1 0 . 0 0 1  (5.14) This type of transformation is called an orthogonal coordinate transformation (OCT). A matrix R satisfying Eq. 5.13 is said to be an orthogonal matrix. That is, an orthogonal matrix is one whose inverse is equal to the transpose. R is sometimes called a rotation matrix. For example, for the coordinate rotation shown in Fig. 49, in 3-D,   cos θ sin θ 0 (5.15) R =  − sin θ cos θ 0  . 0 0 1 In 2-D, R= and   cos θ sin θ − sin θ cos θ  vx = v̄x cos θ − v̄y sin θ vy = v̄x sin θ + v̄y cos θ. (5.16) (5.17) We recall that the determinant of a matrix product is equal to the product of the determinants. Also, the determinant of the transpose of a matrix is equal to the determinant of the matrix itself. Thus, from Eq. 5.13, det(RRT ) = (det R)(det RT ) = (det R)2 = det I = 1, (5.18) and we conclude that, for an orthogonal matrix R, det R = ±1. (5.19) The plus sign occurs for rotations, and the minus sign occurs for combinations of rotations and reflections. For example, the orthogonal matrix   1 0 0 R= 0 1 (5.20) 0  0 0 −1 indicates a reflection in the z direction (i.e., the sign of the z component is changed). We note that the length of a vector is unchanged under an orthogonal coordinate transformation, since the square of the length is given by v̄i v̄i = Rij vj Rik vk = δjk vj vk = vj vj , (5.21) where the summation convention was used. That is, the square of the length of a vector is the same in both coordinate systems. 52 To summarize, under an orthogonal coordinate transformation, vectors transform according to the rule 3 X v̄ = Rv or v̄i = Rij vj , (5.22) j=1 where Rij = ēi · ej , (5.23) RRT = RT R = I. (5.24) and 5.1 Tensors A vector which transforms under an orthogonal coordinate transformation according to the rule v̄ = Rv is defined as a tensor of rank 1. A tensor of rank 0 is a scalar (a quantity which is unchanged by an orthogonal coordinate transformation). For example, temperature and pressure are scalars, since T̄ = T and p̄ = p. We now introduce tensors of rank 2. Consider a matrix M = (Mij ) which relates two vectors u and v by 3 X v = Mu or vi = Mij uj (5.25) j=1 (i.e., the result of multiplying a matrix and a vector is a vector). Also, in a rotated coordinate system, v̄ = M̄ū. (5.26) Since both u and v are vectors (tensors of rank 1), Eq. 5.25 implies RT v̄ = MRT ū or v̄ = RMRT ū. (5.27) By comparing Eqs. 5.26 and 5.27, we obtain M̄ = RMRT (5.28) or, in index notation, M̄ij = 3 X 3 X Rik Rjl Mkl , (5.29) k=1 l=1 which is the transformation rule for a tensor of rank 2. In general, a tensor of rank n, which has n indices, transforms under an orthogonal coordinate transformation according to the rule 3 X 3 3 X X Āij···k = ··· Rip Rjq · · · Rkr Apq···r . (5.30) p=1 q=1 r=1 53 5.2 Examples of Tensors 1. Stress and strain in elasticity The stress tensor σ is  σ11 σ12 σ13 σ =  σ21 σ22 σ23  , σ31 σ32 σ33  (5.31) where σ11 , σ22 , σ33 are the direct (normal) stresses, and σ12 , σ13 , and σ23 are the shear stresses. The corresponding strain tensor is  ε11 ε12 ε13 ε =  ε21 ε22 ε23  , ε31 ε32 ε33  (5.32) where, in terms of displacements, 1 1 εij = (ui,j + uj,i ) = 2 2  ∂ui ∂uj + ∂xj ∂xi  . (5.33) The shear strains in this tensor are equal to half the corresponding engineering shear strains. Both σ and ε transform like second rank tensors. 2. Generalized Hooke’s law According to Hooke’s law in elasticity, the extension in a bar subjected to an axial force is proportional to the force, or stress is proportional to strain. In 1-D, σ = Eε, (5.34) where E is Young’s modulus, an experimentally determined material property. In general three-dimensional elasticity, there are nine components of stress σij and nine components of strain εij . According to generalized Hooke’s law, each stress component can be written as a linear combination of the nine strain components: σij = cijkl εkl , (5.35) where the 81 material constants cijkl are experimentally determined, and the summation convention is being used. We now prove that cijkl is a tensor of rank 4. We can write Eq. 5.35 in terms of stress and strain in a second orthogonal coordinate system as Rki Rlj σ̄kl = cijkl Rmk Rnl ε̄mn . (5.36) If we multiply both sides of this equation by Rpj Roi , and sum repeated indices, we obtain Rpj Roi Rki Rlj σ̄kl = Roi Rpj Rmk Rnl cijkl ε̄mn , 54 (5.37) or, because R is an orthogonal matrix, δok δpl σ̄kl = σ̄op = Roi Rpj Rmk Rnl cijkl ε̄mn . (5.38) Since, in the second coordinate system, σ̄op = c̄opmn ε̄mn , (5.39) c̄opmn = Roi Rpj Rmk Rnl cijkl , (5.40) we conclude that which proves that cijkl is a tensor of rank 4. 3. Stiffness matrix in finite element analysis In the finite element method for linear analysis of structures, the forces F acting on an object in static equilibrium are a linear combination of the displacements u (or vice versa): Ku = F, (5.41) where K is referred to as the stiffness matrix (with dimensions of force/displacement). In this equation, u and F contain several subvectors, since u and F are the displacement and force vectors for all grid points, i.e.,     Fa  ua           Fb     ub  (5.42) , F= u= Fc  uc             ..   ..  . . for grid points a, b, c, . . ., where ūa = Ra ua , ūb = Rb ub , · · · . (5.43) Thus,    ūa       ūb     ū = = ū   c      ..   . Ra 0 0 0 Rb 0 0 0 Rc .. .. .. . . . where Γ is an orthogonal block-diagonal matrix  Ra 0  0 Rb  Γ=  0 0 .. .. . . and            ua ub uc .. .      = Γu, (5.44)     consisting of rotation matrices R:  0 ··· 0 ···   , Rc · · ·   .. . . . . ΓT Γ = I. 55 ··· ··· ··· ... (5.45) (5.46) q 4 ✻y θ ȳ ■ ❅ ❅q ✒ DOF x̄ 2 ✲ ■ ❅ x ❅q ✒ ■ ❅ ❅q ✒ 3 θ 1 Figure 50: Element Coordinate System for Pin-Jointed Rod. Similarly, F̄ = ΓF. (5.47) K̄ū = F̄, (5.48) K̄Γu = ΓF (5.49) Thus, if or
ΓT K̄Γ u = F. (5.50) K = ΓT K̄Γ. (5.51) That is, the stiffness matrix transforms like other tensors of rank 2: We illustrate the transformation of a finite element stiffness matrix by transforming the stiffness matrix for the pin-jointed rod element from a local element coordinate system to a global Cartesian system. Consider the rod shown in Fig. 50. For this element, the 4 × 4 2-D stiffness matrix in the element coordinate system is   1 0 −1 0  AE  0 0   0 0 (5.52) K̄ = ,  1 0  L  −1 0 0 0 0 0 where A is the cross-sectional area of the rod, E is Young’s modulus of the rod material, and L is the rod length. In the global coordinate system, K = ΓT K̄Γ, (5.53) where the 4 × 4 transformation matrix is Γ= and the rotation matrix is R=   R 0 0 R  , cos θ sin θ − sin θ cos θ 56 (5.54)  . (5.55) Thus, K= = AE L AE L   0 −1 0 c s 0  0 0 0   −s c 0  0 1 0  0 0 c 0 0 0 0 0 −s  c −s 0 0 s c 0 0 0 0 c −s 0 0 s c   c2 cs −c2 −cs  cs s2 −cs −s2  , −c2 −cs c2 cs  −cs −s2 cs s2             1 0 −1 0 0 0 s c      (5.56) where c = cos θ and s = sin θ. This result agrees with that found earlier in Eq. 4.27. 5.3 Isotropic Tensors An isotropic tensor is a tensor which is independent of coordinate system (i.e., invariant under an orthogonal coordinate transformation). The Kronecker delta δij is a second rank tensor and isotropic, since δ̄ij = δij , and Ī = RIRT = RRT = I. (5.57) That is, the identity matrix I is invariant under an orthogonal coordinate transformation. It can be shown in tensor analysis that δij is the only isotropic tensor of rank 2 and, moreover, δij is the characteristic tensor for all isotropic tensors: Rank 1 2 3 4 odd Isotropic Tensors none cδij none aδij δkl + bδik δjl + cδil δjk none That is, all isotropic tensors of rank 4 must be of the form shown above, which has three constants. For example, in generalized Hooke’s law (Eq. 5.35), the material property tensor cijkl has 34 = 81 constants (assuming no symmetry). For an isotropic material, cijkl must be an isotropic tensor of rank 4, thus implying at most three material constants (on the basis of tensor analysis alone). The actual number of independent material constants for an isotropic material turns out to be two rather than three, a result which depends on the existence of a strain energy function, which implies the additional symmetry cijkl = cklij . 6 Calculus of Variations Recall from calculus that a function of one variable y = f (x) attains a stationary value (minimum, maximum, or neutral) at the point x = x0 if the derivative f ′ (x) = 0 at x = x0 . 57 ✻ r ✻ r ✲ ✻ r ✲ (a) Minimum (b) Maximum ✲ (c) Neutral Figure 51: Minimum, Maximum, and Neutral Stationary Values. Alternatively, the first variation of f (which is similar to a differential) must vanish for arbitrary changes δx in x:   df δf = δx = 0. (6.1) dx The second derivative determines what kind of stationary value one has: d2 f > 0 at x = x0 dx2 d2 f maximum: < 0 at x = x0 dx2 d2 f neutral: = 0 at x = x0 , dx2 minimum: (6.2) (6.3) (6.4) as shown in Fig. 51. For a function of two variables z = f (x, y), z is stationary at (x0 , y0 ) if ∂f ∂f = 0 and = 0 at (x0 , y0 ). ∂x ∂y Alternatively, for a stationary value,     ∂f ∂f δf = δx + δy = 0 at (x0 , y0 ). ∂x ∂y A function with n independent variables f (x1 , x2 , . . . , xn ) is stationary at a point P if n X ∂f δxi = 0 at P δf = ∂x i i=1 for arbitrary variations δxi . Hence, ∂f = 0 at P ∂xi (i = 1, 2, . . . , n). Variational calculus is concerned with finding stationary values, not of functions, but of functionals. In general, a functional is a function of a function. More precisely, a functional is an integral which has a specific numerical value for each function that is substituted into it. 58 Consider the functional I(φ) = Z b F (x, φ, φx ) dx, (6.5) a where x is the independent variable, φ(x) is the dependent variable, and φx = dφ . dx (6.6) The variation in I due to an arbitrary small change in φ is  Z b Z b ∂F ∂F δF dx = δI = δφ + δφx dx, ∂φ ∂φx a a (6.7) where, with an interchange of the order of d and δ,   dφ d δφx = δ = (δφ). dx dx (6.8) With this equation, the second term in Eq. 6.7 can be integrated by parts to obtain   Z b Z b Z b b ∂F ∂F d ∂F d ∂F dx. δφx dx = (δφ) dx = δφ − δφ ∂φx dx ∂φx a ∂φx a ∂φx dx a a Thus, δI = Z b a d ∂F − ∂φ dx  ∂F ∂φx  δφ dx +  ∂F δφ ∂φx  (6.9) b . (6.10) a Since δφ is arbitrary (within certain limits of admissibility), δI = 0 implies that both terms in Eq. 6.10 must vanish. Moreover, for arbitrary δφ, a zero integral in Eq. 6.10 implies a zero integrand:   ∂F ∂F d = 0, (6.11) − dx ∂φx ∂φ which is referred to as the Euler-Lagrange equation. The second term in Eq. 6.10 must also vanish for arbitrary δφ: b  ∂F δφ = 0. ∂φx a (6.12) If φ is specified at the boundaries x = a and x = b, δφ(a) = δφ(b) = 0, (6.13) and Eq. 6.12 is automatically satisfied. This type of boundary condition is called a geometric or essential boundary condition. If φ is not specified at the boundaries x = a and x = b, Eq. 6.12 implies ∂F (b) ∂F (a) = = 0. (6.14) ∂φx ∂φx This type of boundary condition is called a natural boundary condition. 59 s (1, 1) y ✻ r ds dy r dx (0, 0) ✲ x Figure 52: Curve of Minimum Length Between Two Points. 6.1 Example 1: The Shortest Distance Between Two Points We wish to find the equation of the curve y(x) along which the distance from the origin (0, 0) to (1, 1) in the xy-plane is least. Consider the curve in Fig. 52. The differential arc length along the curve is given by p p (6.15) ds = dx2 + dy 2 = 1 + (y ′ )2 dx. We seek the curve y(x) which minimizes the total arc length Z 1p I(y) = 1 + (y ′ )2 dx, (6.16) 0 where y(0) = 0 and y(1) = 1. For this problem, 1/2  , F (x, y, y ′ ) = 1 + (y ′ )2 (6.17) which depends only on y ′ explicitly. Thus, the Euler-Lagrange equation, Eq. 6.11, is   d ∂F = 0, (6.18) dx ∂y ′ where  y′ 1 ∂F ′ ′ 2 −1/2 2y = 1 + (y ) = . (6.19) ∂y ′ 2 [1 + (y ′ )2 ]1/2 Hence, y′ = c, (6.20) [1 + (y ′ )2 ]1/2 where c is a constant of integration. If we solve this equation for y ′ , we obtain r c2 ′ = a, (6.21) y = 1 − c2 where a is another constant. Thus, y(x) = ax + b, (6.22) y(x) = x, (6.23) and, with the boundary conditions, which is the straight line joining (0, 0) and (1, 1), as expected. 60 O ✲ y y(x) sA x ❄ Figure 53: The Brachistochrone Problem. 6.2 Example 2: The Brachistochrone We wish to determine the path y(x) from the origin O to the point A along which a bead will slide under the force of gravity and without friction in the shortest time (Fig. 53). Assume that the bead starts from rest. The velocity v along the path y(x) is p p 1 + (y ′ )2 dx dx2 + dy 2 ds = = , (6.24) v= dt dt dt where s denotes distance along the path. Thus, p 1 + (y ′ )2 dx dt = . v (6.25) As the bead slides down the path, potential energy is converted to kinetic energy. At any location x, the kinetic energy equals the reduction in potential energy: 1 2 mv = mgx 2 (6.26) or v= Thus, from Eq. 6.25, dt = s p 2gx. 1 + (y ′ )2 dx. 2gx The total time for the bead to fall from O to A is then Z xA s 1 + (y ′ )2 t= dx, 2gx 0 (6.27) (6.28) (6.29) which is the integral to be minimized. From Eq. 6.11, the Euler-Lagrange equation is   ∂F d ∂F − = 0, (6.30) ′ dx ∂y ∂y where F (x, y, y ′ ) = s 61 1 + (y ′ )2 . 2gx (6.31) Thus, the Euler-Lagrange equation becomes "  −1/2 ′ # d 1 1 + (y ′ )2 2y =0 dx 2 2gx 2gx or y′ d dx {x [1 + (y ′ )2 ]}1/2 ! = 0. (6.32) (6.33) To solve this equation, we integrate Eq. 6.33 to obtain y′ {x [1 + (y ′ )2 ]}1/2 = c, (6.34) where c is a constant of integration. We then square both sides of this equation, and solve for y ′ :   (6.35) (y ′ )2 = c2 x 1 + (y ′ )2 or r c2 x dy y = = . dx 1 − c2 x This equation can be integrated using the change of variable ′ 1 1 sin2 (θ/2), dx = 2 sin(θ/2) cos(θ/2) dθ. 2 c c The integral implied by Eq. 6.36 then transforms to Z Z 1 sin(θ/2) 1 · 2 sin(θ/2) cos(θ/2) dθ = 2 sin2 (θ/2) dθ y= cos(θ/2) c c Z 1 1 = 2 (1 − cos θ) dθ = 2 (θ − sin θ) + y0 , 2c 2c x= (6.36) (6.37) (6.38) (6.39) where y0 is a constant of integration. Since the curve must pass through the origin, we must have y0 = 0. Also, from Eq. 6.37, 1 (1 − cos θ). (6.40) 2c2 If we then replace the constant a = 1/(2c2 ), we obtain the final result in parametric form  x = a(1 − cos θ) (6.41) y = a(θ − sin θ), x= which is a cycloid generated by the motion of a fixed point on the circumference of a circle of radius a which rolls on the positive side of the line x = 0. To solve Eq. 6.41 for a specific cycloid (defined by the two end points), one can first eliminate the circle radius a from Eq. 6.41 to solve (iteratively) for θA (the value of θ at Point A), after which either of the two equations in Eq. 6.41 can be used to determine the constant a. The resulting cycloids for several end points are shown in Fig. 54. This brachistochrone solution is valid for any end point which the bead can reach. Thus, the end point must not be higher than the starting point. The end point may be at the same elevation since, without friction, there is no loss of energy during the motion. 62 1 2 1 r r 2 r ❄ x 3 4 ✲ y r Figure 54: Several Brachistochrone Solutions. 6.3 Constraint Conditions Suppose we want to extremize the functional Z b F (x, φ, φx ) dx I(φ) = (6.42) a subject to the additional constraint condition Z b G(x, φ, φx ) dx = J, (6.43) a where J is a specified constant. We recall that the Euler-Lagrange equation was found by requiring that δI = 0. However, since J is a constant, δJ = 0. Thus, δ(I + λJ) = 0, (6.44) where λ is an unknown constant referred to as a Lagrange multiplier. Thus, to enforce the constraint in Eq. 6.43, we can replace F in the Euler-Lagrange equation with H = F + λG. 6.4 (6.45) Example 3: A Constrained Minimization Problem We wish to find the function y(x) which minimizes the integral Z 1 (y ′ )2 dx I(y) = (6.46) 0 subject to the end conditions y(0) = y(1) = 0 and the constraint Z 1 y dx = 1. 0 That is, the area under the curve y(x) is unity (Fig. 55). For this problem, 63 (6.47) y ✻ Area=1 r r (0, 0) (1, 0) ✲ x Figure 55: A Constrained Minimization Problem. H(x, y, y ′ ) = (y ′ )2 + λy. (6.48) The Euler-Lagrange equation, d dx yields or  ∂H ∂y ′  − ∂H = 0, ∂y (6.49) d (2y ′ ) − λ = 0 dx (6.50) y ′′ = λ/2. (6.51) After integrating this equation, we obtain y= λ 2 x + Ax + B. 4 (6.52) With the boundary conditions y(0) = y(1) = 0, we obtain B = 0 and A = −λ/4. Thus, y = −λx(1 − x)/4. The area constraint, Eq. 6.47, is used to evaluate the Lagrange multiplier λ:  1 Z 1 Z 1
λ λ x2 x3 λ 2 y dx = − x − x dx = − 1= − =− . 4 2 3 0 24 0 4 0 (6.53) (6.54) Thus, with λ = −24, we obtain the parabola y(x) = 6x(1 − x). 6.5 (6.55) Functions of Several Independent Variables Consider I(φ) = Z F (x, y, z, φ, φx , φy , φz ) dV, (6.56) V a functional with three independent variables (x, y, z). Note that, in 3-D, the integration is a volume integral. The variation in I due to an arbitrary small change in φ is  Z  ∂F ∂F ∂F ∂F (6.57) δφ + δφx + δφy + δφz dV, δI = ∂φ ∂φx ∂φy ∂φz V 64 where, with an interchange of the order of ∂ and δ,  Z  ∂F ∂ ∂F ∂ ∂F ∂ ∂F δφ + (δφ) + (δφ) + (δφ) dV. δI = ∂φ ∂φx ∂x ∂φy ∂y ∂φz ∂z V (6.58) To perform a three-dimensional integration by parts, we note that the second term in Eq. 6.58 can be expanded to yield     ∂F ∂ ∂ ∂ ∂F ∂F δφ. (6.59) (δφ) = δφ − ∂φx ∂x ∂x ∂φx ∂x ∂φx If we then expand the third and fourth terms similarly, Eq. 6.58 becomes, after regrouping,       Z  ∂F ∂ ∂ ∂ ∂F ∂F ∂F δI = δφ − δφ − δφ δφ − ∂φ ∂x ∂φx ∂y ∂φy ∂z ∂φz V       ∂ ∂ ∂ ∂F ∂F ∂F + δφ + δφ + δφ dV. (6.60) ∂x ∂φx ∂y ∂φy ∂z ∂φz The last three terms in this integral can then be converted to a surface integral using the divergence theorem, which states that, for a vector field f , I Z f · n dS, (6.61) ∇ · f dV = S V where S is the closed surface which encloses the volume V . Eq. 6.60 then becomes       Z  ∂F ∂ ∂ ∂ ∂F ∂F ∂F δI = − − δφ dV − ∂φ ∂x ∂φx ∂y ∂φy ∂z ∂φz V  I  ∂F ∂F ∂F nx + ny + nz δφ dS, + ∂φx ∂φy ∂φz S (6.62) where n = (nx , ny , nz ) is the unit outward normal on the surface. Since δI = 0, both integrals in this equation must vanish for arbitrary δφ. It therefore follows that the integrand in the volume integral must also vanish to yield the Euler-Lagrange equation for three independent variables:       ∂ ∂ ∂F ∂F ∂F ∂F ∂ + + − = 0. (6.63) ∂x ∂φx ∂y ∂φy ∂z ∂φz ∂φ If φ is specified on the boundary S, δφ = 0 on S, and the boundary integral in Eq. 6.62 automatically vanishes. If φ is not specified on S, the integrand in the boundary integral must vanish: ∂F ∂F ∂F nx + ny + nz = 0 on S. (6.64) ∂φx ∂φy ∂φz This is the natural boundary condition when φ is not specified on the boundary. 65 6.6 Example 4: Poisson’s Equation Consider the functional I(φ) = Z  V  1 2 2 2 (φ + φy + φz ) − Qφ dV, 2 x (6.65) in which case 1 F (x, y, z, φ, φx , φy , φz ) = (φ2x + φ2y + φ2z ) − Qφ. 2 The Euler-Lagrange equation, Eq. 6.63, implies (6.66) ∂ ∂ ∂ (φx ) + (φy ) + (φz ) − (−Q) = 0. ∂x ∂y ∂z (6.67) φxx + φyy + φzz = −Q (6.68) ∇2 φ = −Q, (6.69) That is, or which is Poisson’s equation. Thus, we have shown that solving Poisson’s equation is equivalent to minimizing the functional I in Eq. 6.65. In general, a key conclusion of this discussion of variational calculus is that solving a differential equation is equivalent to minimizing some functional involving an integral. 6.7 Functions of Several Dependent Variables Consider the functional I(φ1 , φ2 , φ3 ) = Z b a F (x, φ1 , φ2 , φ3 , φ′1 , φ′2 , φ′3 ) dx, (6.70) where x is the independent variable, and φ1 (x), φ2 (x), and φ3 (x) are the dependent variables. It can be shown that the generalization of the Euler-Lagrange equation, Eq. 6.11, for this case is   ∂F d ∂F − =0 (6.71) ′ dx ∂φ1 ∂φ1   ∂F d ∂F − =0 (6.72) ′ dx ∂φ2 ∂φ2   ∂F d ∂F − = 0. (6.73) ′ dx ∂φ3 ∂φ3 7 Variational Approach to the Finite Element Method In modern engineering analysis, one of the most important applications of the energy theorems, such as the minimum potential energy theorem in elasticity, is the finite element method, a numerical procedure for solving partial differential equations. For linear equations, finite element solutions are often based on variational methods. 66 7.1 Index Notation and Summation Convention Let a be the vector a = a1 e1 + a2 e2 + a3 e3 , (7.1) where ei is the unit vector in the ith Cartesian coordinate direction, and ai is the ith component of a (the projection of a on ei ). Consider a second vector b = b1 e1 + b2 e2 + b 3 e3 . (7.2) Then, the dot product (scalar product) is a · b = a1 b1 + a2 b2 + a3 b3 = 3 X ai bi . (7.3) i=1 Using the summation convection, we can omit the summation symbol and write a · b = ai bi , (7.4) where, if a subscript appears exactly twice, a summation is implied over the range. The range is 1,2,3 in 3-D and 1,2 in 2-D. Let f (x1 , x2 , x3 ) be a scalar function of the Cartesian coordinates x1 , x2 , x3 . We define the shorthand comma notation for derivatives: ∂f = f,i . ∂xi (7.5) That is, the subscript “, i” denotes the partial derivative with respect to the ith Cartesian coordinate direction. Using the comma notation and the summation convention, a variety of familiar quantities can be written in compact form. For example, the gradient can be written ∇f = ∂f ∂f ∂f e1 + e2 + e3 = f,i ei . ∂x1 ∂x2 ∂x3 (7.6) For a vector-valued function F(x1 , x2 , x3 ), the divergence is written ∇·F= ∂F1 ∂F2 ∂F3 + + = Fi,i , ∂x1 ∂x2 ∂x3 (7.7) and the Laplacian of the scalar function f (x1 , x2 , x3 ) is ∇2 f = ∇ · ∇f = ∂ 2f ∂ 2f ∂ 2f + + = f,ii . ∂x21 ∂x22 ∂x23 (7.8) The divergence theorem states that, for any vector field F(x1 , x2 , x3 ) defined in a volume V bounded by a closed surface S, Z I ∇ · F dV = F · n dS (7.9) V S 67 or, in index notation, Z Fi,i dV = V I Fi ni dS. The normal derivative can be written in index notation as ∂φ = ∇φ · n = φ,i ni . ∂n The dot product of two gradients can be written ∇φ · ∇φ = (φ,i ei ) · (φ,j ej ) = φ,i φ,j ei · ej = φ,i φ,j δij = φ,i φ,i , where δij is the Kronecker delta defined as  1, i = j, δij = 0, i = 6 j. 7.2 (7.10) S (7.11) (7.12) (7.13) Deriving Variational Principles For each partial differential equation of interest, one needs a functional which a solution makes stationary. Given the functional, it is generally easy to see what partial differential equation corresponds to it. The harder problem is to start with the equation and derive the variational principle (i.e., derive the functional which is to be minimized). To simplify this discussion, we will consider only scalar, rather than vector, problems. A scalar problem has one dependent variable and one partial differential equation, whereas a vector problem has several dependent variables coupled to each other through a system of partial differential equations. Consider Poisson’s equation subject to both Dirichlet and Neumann boundary conditions,  2   ∇ φ + f = 0 in V,  φ = φ0 on S1 , (7.14)  ∂φ   + g = 0 on S2 . ∂n This problem arises, for example, in (1) steady-state heat conduction, where the temperature and heat flux are specified on different parts of the boundary, (2) potential fluid flow, where the velocity potential and velocity are specified on different parts of the boundary, and (3) torsion in elasticity. Recall that, in index notation, ∇2 φ = φ,ii . We wish to find a functional, I(φ) say, whose first variation δI vanishes for φ satisfying the above boundary value problem. From Eq. 7.14a, Z (7.15) 0 = (φ,ii + f ) δφ dV Z ZV δ(f φ) dV (7.16) = [(φ,i δφ),i − φ,i δφ,i ] dV + V V Z I Z 1 φ,i ni δφ dS − = f φ dV (7.17) δ(φ,i φ,i ) dV + δ V 2 S V Z I Z 1 f φ dV, (7.18) = (∇φ · n) δφ dS − δ φ,i φ,i dV + δ V 2 S V 68 where ∇φ · n = ∂φ/∂n = −g on S2 , and, since φ is specified on S1 , δφ = 0 on S1 . Then, Z Z Z 1 φ,i φ,i dV + δ f φ dV (7.19) 0=− g δφ dS − δ V S2 V 2 Z Z Z 1 gφ dS − δ = −δ f φ dV (7.20) φ,i φ,i dV + δ V 2 S2 V   Z  Z 1 gφ dS = −δI(φ). (7.21) φ,i φ,i − f φ dV + = −δ 2 S2 V Thus, the functional for this boundary value problem is  Z  Z 1 I(φ) = gφ dS φ,i φ,i − f φ dV + 2 V S2 (7.22) or, in vector notation, I(φ) = Z  V  Z 1 gφ dS ∇φ · ∇φ − f φ dV + 2 S2 or, in expanded form, ) Z Z ( "  2  2  2 # 1 ∂φ ∂φ ∂φ − f φ dV + gφ dS. I(φ) = + + 2 ∂x ∂y ∂z S2 V (7.23) (7.24) If we were given the functional, Eq. 7.22, we could determine which partial differential equation corresponds to it by computing and setting to zero the first variation of the functional. From Eq. 7.22,  Z Z  1 1 g δφ dS (7.25) δφ,i φ,i + φ,i δφ,i − f δφ dV + δI = 2 2 S2 V Z Z g δφ dS (7.26) = (φ,i δφ,i − f δφ) dV + S2 V Z Z Z g δφ dS (7.27) f δφ dV + = [(φ,i δφ),i − φ,ii δφ] dV − V Z S2 Z IV g δφ dS, (7.28) φ,i ni δφ dS − (φ,ii + f ) δφ dV + = S S2 V where δφ = 0 on S1 , and φ,i ni = ∇φ · n = ∂φ/∂n. Hence,  Z  Z ∂φ δI = + g δφ dS − (φ,ii + f ) δφ dV. ∂n S2 V (7.29) Stationary I (δI = 0) for arbitrary admissible δφ thus yields the original partial differential equation and Neumann boundary condition, Eq. 7.14. 69 y✻ ✲ x Figure 56: Two-Dimensional Finite Element Mesh. y✻ 3 ✄❅ ✄ ❅ ✄ ❅ ✘ ✄✘✘✘ 2 1 ✲ x Figure 57: Triangular Finite Element. 7.3 Shape Functions Consider a two-dimensional field problem for which we seek the function φ(x, y) satisfying some (unspecified) partial differential equation in a domain. Let us subdivide the domain into triangular finite elements, as shown in Fig. 56. A typical element, shown in Fig. 57, has three degrees of freedom (DOF): φ1 , φ2 , and φ3 . The three DOF are the values of the fundamental unknown φ at the three grid points. (The number of DOF for an element represents the number of pieces of data needed to determine uniquely the solution for the element.) For numerical solution, we approximate φ(x, y) using a polynomial in two variables having the same number of terms as the number of DOF. Thus, we assume, for each element, φ(x, y) = a1 + a2 x + a3 y, (7.30) where a1 , a2 , and a3 are unknown coefficients. Eq. 7.30 is a complete linear polynomial in two variables. We can evaluate this equation at the three grid points to obtain      1 x1 y1  a1   φ1  . (7.31) =  1 x2 y2  a2 φ    2  a3 1 x3 y3 φ3 This matrix equation can be inverted to yield      x2 y3 − x3 y2 x3 y1 − x1 y3 x1 y2 − x2 y1  φ1   a1  1  , = y2 − y3 y3 − y1 y1 − y2  φ2 a2    2A  x3 − x2 x1 − x3 x2 − x1 φ3 a3 (7.32) = x2 y3 − x3 y2 − x1 y3 + x3 y1 + x1 y2 − x2 y1 . (7.33) where 2A = 1 x1 y1 1 x2 y2 1 x3 y3 70 Note that the area of the triangle is |A|. The determinant 2A is positive or negative depending on whether the grid point ordering in the triangle is counter-clockwise or clockwise, respectively. Since Eq. 7.31 is invertible, we can interpret the element DOF as either the grid point values φi , which have physical meaning, or the nonphysical polynomial coefficients ai . From Eq. 7.30, the scalar unknown φ can then be written in the matrix form   a1     φ(x, y) = 1 x y (7.34) a  2  a3    x2 y3 − x3 y2 x3 y1 − x1 y3 x1 y2 − x2 y1  φ 1   1  = (7.35) 1 x y  y2 − y3 y3 − y1 y1 − y 2  φ 2   2A x3 − x2 x1 − x3 x2 − x1 φ3   φ   1  = N1 (x, y) N2 (x, y) N3 (x, y) (7.36) φ  2  φ3 or φ(x, y) = N1 (x, y)φ1 + N2 (x, y)φ2 + N3 (x, y)φ3 = Ni φi , where   N1 (x, y) = [(x2 y3 − x3 y2 ) + (y2 − y3 )x + (x3 − x2 )y]/(2A) N (x, y) = [(x3 y1 − x1 y3 ) + (y3 − y1 )x + (x1 − x3 )y]/(2A)  2 N3 (x, y) = [(x1 y2 − x2 y1 ) + (y1 − y2 )x + (x2 − x1 )y]/(2A). (7.37) (7.38) Notice that all the subscripts in this equation form a cyclic permutation of the numbers 1, 2, 3. That is, knowing N1 , we can write down N2 and N3 by a cyclic permutation of the subscripts. Alternatively, if we let i, j, and k denote the three grid points of the triangle, the above equation can be written in the compact form Ni (x, y) = 1 [(xj yk − xk yj ) + (yj − yk )x + (xk − xj )y], 2A (7.39) where (i, j, k) can be selected to be any cyclic permutation of (1, 2, 3) such as (1, 2, 3), (2, 3, 1), or (3, 1, 2). In general, the functions Ni are called shape functions or interpolation functions and are used extensively in finite element theory. Eq. 7.37 implies that N1 can be thought of as the shape function associated with Point 1, since N1 is the form (or “shape”) that φ would take if φ1 = 1 and φ2 = φ3 = 0. In general, Ni is the shape function for Point i. From Eq. 7.37, if φ1 6= 0, and φ2 = φ3 = 0, φ(x, y) = N1 (x, y)φ1 . (7.40) φ1 = φ(x1 , y1 ) = N1 (x1 , y1 )φ1 (7.41) N1 (x1 , y1 ) = 1. (7.42) In particular, at Point 1, or 71 1 ✞ q ✝ 2☎ q L ✆ ✲ x Figure 58: Axial Member (Pin-Jointed Truss Element). Also, at Point 2, 0 = φ2 = φ(x2 , y2 ) = N1 (x2 , y2 )φ1 (7.43) N1 (x2 , y2 ) = 0. (7.44) N1 (x3 , y3 ) = 0. (7.45) or Similarly, That is, a shape function takes the value unity at its own grid point and zero at the other grid points in an element: Ni (xj , yj ) = δij . (7.46) Eq. 7.38 gives the shape functions for the linear triangle. An interesting observation is that, if the shape functions are evaluated at the triangle centroid, N1 = N2 = N3 = 1/3, since each grid point has equal influence on the centroid. To prove this result, substitute the coordinates of the centroid, x̄ = (x1 + x2 + x3 )/3, ȳ = (y1 + y2 + y3 )/3, (7.47) into N1 in Eq. 7.38, and use Eq. 7.33: 1 N1 (x̄, ȳ) = [(x2 y3 − x3 y2 ) + (y2 − y3 )x̄ + (x3 − x2 )ȳ]/(2A) = . 3 (7.48) The preceding discussion assumed an element having three grid points. In general, for an element having r points, φ(x, y) = N1 (x, y)φ1 + N2 (x, y)φ2 + · · · + Nr (x, y)φr = Ni φi . (7.49) This is a convenient way to represent the finite element approximation within an element, since, once the grid point values are known, φ can be evaluated anywhere in the element. We note that, since the shape functions for the 3-point triangle are linear in x and y, the gradients in the x or y directions are constant. Thus, from Eq. 7.37, ∂N1 ∂φ = φ1 + ∂x ∂x ∂φ ∂N1 = φ1 + ∂y ∂y ∂N2 φ2 + ∂x ∂N2 φ2 + ∂y ∂N3 φ3 = [(y2 − y3 )φ1 + (y3 − y1 )φ2 + (y1 − y2 )φ3 ]/(2A) (7.50) ∂x ∂N3 φ3 = [(x3 − x2 )φ1 + (x1 − x3 )φ2 + (x2 − x1 )φ3 ]/(2A). (7.51) ∂y A constant gradient within any element implies that many small elements would have to be used wherever φ is changing rapidly. Example. The displacement u(x) of an axial structural member (a pin-jointed truss member) (Fig. 58) is given by u(x) = N1 (x)u1 + N2 (x)u2 , 72 (7.52) where u1 and u2 are the axial displacements at the two end points. For a linear variation of displacement along the length, it follows that Ni must be a linear function which is unity at Point i and zero at the other end. Thus,   N1 (x) = 1 − x , L (7.53)  N2 (x) = x L or  x x u(x) = 1 − u1 + u2 . (7.54) L L 7.4 Variational Approach Consider Poisson’s equation in a two-dimensional domain subject to both Dirichlet and Neumann boundary conditions. We consider a slight generalization of the problem of Eq. 7.14:  2 ∂ φ ∂ 2φ     ∂x2 + ∂y 2 + f = 0 in A, (7.55) φ = φ0 on S1 ,    ∂φ  + g + hφ = 0 on S2 , ∂n where S1 and S2 are curves in 2-D, and f , g, and h may depend on position. The difference between this problem and that of Eq. 7.14 is that the h term has been added to the boundary condition on S2 , where the gradient ∂φ/∂n is specified. The h term could arise in a variety of physical situations, including heat transfer due to convection (where the heat flux is proportional to temperature) and free surface flow problems (where the free surface and radiation boundary conditions are both of this type). Assume that the domain has been subdivided into a mesh of triangular finite elements similar to that shown in Fig. 56. The functional which must be minimized for this boundary value problem is similar to that of Eq. 7.24 except that an additional term must be added to account for the h term in the boundary condition on S2 . It can be shown that the functional which must be minimized for this problem is )  Z ( "  2  2 # Z  1 2 ∂φ ∂φ 1 gφ + hφ dS, + − f φ dA + I(φ) = (7.56) 2 ∂x ∂y 2 A S2 where A is the domain. With a finite element discretization, we do not define a single smooth function φ over the entire domain, but instead define φ over individual elements. Thus, since I is an integral, it can be evaluated by summing over the elements: I = I1 + I2 + I3 + · · · = E X Ie , (7.57) e=1 where E is the number of elements. The variation of I is also computed by summing over the elements: E X δI = δIe , (7.58) e=1 73 29 ∂φ ✒ ∂n ✄❅ ✄ ❧ ❅ ✄ 8 ✘❅ ✘ 17 ✄✘✘ 15 29✘✘✘✘✄ 35 ❅ 9❧✄ ∂φ ✠ ❅❅✄✄ 17 ∂n Figure 59: Neumannn Boundary Condition at Internal Boundary. which must vanish for I to be a minimum. We thus can consider a typical element. For one element, in index notation,   Z  Z  1 1 2 Ie = gφ + hφ dS. φ,k φ,k − f φ dA + 2 2 A S2 (7.59) The last two terms in this equation are, from Eq. 7.55c, integrals of the form Z ∂φ φ dS. S2 ∂n As can be seen from Fig. 59, for two elements which share a common edge, the unknown φ is continuous along that edge, and the normal derivative ∂φ/∂n is of equal magnitude and opposite sign, so that the individual contributions cancel each other out. Thus, the last two terms in the functional make a nonzero contribution to the functional only if an element abuts S2 (i.e., if one edge of an element coincides with S2 ). The degrees of freedom which define the function φ over the entire domain are the nodal (grid point) values φi , since, given all the φi , φ is known everywhere using the element shape functions. Therefore, to minimize I, we differentiate with respect to each φi , and set the result to zero:   Z  Z  ∂Ie ∂φ 1 ∂ ∂φ ∂φ 1 ∂ g dA + dS (7.60) = (φ,k )φ,k + φ,k (φ,k ) − f + hφ ∂φi 2 ∂φi ∂φi ∂φi ∂φi S2 A 2 ∂φi   Z  Z  ∂ ∂φ ∂φ ∂φ = φ,k dA + dS. (7.61) g (φ,k ) − f + hφ ∂φi ∂φi ∂φi ∂φi A S2 where φ = Nj φj implies φ,k = (Nj φj ),k = Nj,k φj , ∂ ∂φj ∂ (φ,k ) = (Nj,k φj ) = Nj,k = Nj,k δij = Ni,k , ∂φi ∂φi ∂φi (7.62) (7.63) and ∂ ∂φj ∂φ = (Nj φj ) = Nj = Nj δij = Ni . ∂φi ∂φi ∂φi We now substitute these last three equations into Eq. 7.61 to obtain Z Z ∂Ie (gNi + hNj φj Ni ) dS = (Nj,k φj Ni,k − f Ni ) dA + ∂φi S2 A   Z  Z Z Z hNi Nj dS φj gNi dS + f Ni dA − Ni,k Nj,k dA φj − = A S2 A = Kij φj − Fi + Hij φj , 74 (7.64) (7.65) (7.66) S2 (7.67) 29✘✘✘✘✄ 35 ✄❅ 9❧✄ ✄ ❧ ❅ ✄ ✄ 8 ✘❅ ✘✄ 17 ✄✘✘ 15 Figure 60: Two Adjacent Finite Elements. where, for each element, Kije Fie = Z A Hije = Z Ni,k Nj,k dA (7.68) Z (7.69) A f Ni dA − = Z gNi dS S2 hNi Nj dS. (7.70) S2 The second term of the “load” F is present only if Point i is on S2 , and the matrix entries in H apply only if Points i, j are both on the boundary. The two terms in F correspond to the body force and Neumann boundary condition, respectively. The superscript e in the preceding equations indicates that we have computed the contribution for one element. We must combine the contributions for all elements. Consider two adjacent elements, as illustrated in Fig. 60. In that figure, the circled numbers are the element labels. From Eq. 7.67, for Element 8, for the special case h = 0, ∂I8 8 8 8 8 = K15,15 φ15 + K15,17 φ17 + K15,29 φ29 − F15 , ∂φ15 (7.71) ∂I8 8 8 8 8 = K17,15 φ15 + K17,17 φ17 + K17,29 φ29 − F17 , ∂φ17 (7.72) ∂I8 8 8 8 8 = K29,15 φ15 + K29,17 φ17 + K29,29 φ29 − F29 . ∂φ29 (7.73) Similarly, for Element 9, To evaluate ∂I9 9 9 9 9 = K17,17 φ17 + K17,29 φ29 + K17,35 φ35 − F17 , ∂φ17 (7.74) ∂I9 9 9 9 9 = K29,17 φ17 + K29,29 φ29 + K29,35 φ35 − F29 , ∂φ29 (7.75) ∂I9 9 9 9 9 = K35,17 φ17 + K35,29 φ29 + K35,35 φ35 − F35 . ∂φ35 (7.76) X ∂Ie e ∂φi 75 , we sum on e for fixed i. For example, for i = 17, ∂I8 ∂I9 8 8 9 8 9 + =K17,15 φ15 + (K17,17 + K17,17 )φ17 + (K17,29 + K17,29 )φ29 ∂φ17 ∂φ17 9 8 9 + K17,35 φ35 − F17 − F17 . (7.77) This process is then repeated for all grid points and all elements. To minimize the functional, the individual sums are set to zero, resulting in a set of linear algebraic equations of the form Kφ = F (7.78) or, for the more general case where h 6= 0, (K + H)φ = F, (7.79) where φ is the vector of unknown nodal values of φ, and F is the vector of forcing functions at the nodes. Because of the historical developments in structural mechanics, K is sometimes referred to as the stiffness matrix. For each element, the contributions to K, F, and H are Z Ni,k Nj,k dA, (7.80) Kij = A Fi = Z A f Ni dA − Hij = Z Z gNi dS, (7.81) S2 hNi Nj dS. (7.82) S2 The sum on k in the first integral can be expanded to yield  Z  ∂Ni ∂Nj ∂Ni ∂Nj dA. + Kij = ∂x ∂x ∂y ∂y A (7.83) Note that, in three dimensions, Eq. 7.80 still applies, except that the sum extends over the range 1–3, and the integration is over the element volume. Note also in Eq. 7.81 that, if g = 0, there is no contribution to the right-hand side vector F. Thus, to implement the Neumann boundary condition ∂φ/∂n = 0 at a boundary, the boundary is left free. The zero gradient boundary condition is the natural boundary condition for this formulation, since a zero gradient results by default if the boundary is left free. The Dirichlet boundary condition φ = φ0 must be explicitly imposed and is referred to as an essential boundary condition. 7.5 Matrices for Linear Triangle Consider the linear three-node triangle, for which, from Eq. 7.39, the shape functions are given by 1 Ni (x, y) = [(xj yk − xk yj ) + (yj − yk )x + (xk − xj )y], (7.84) 2A 76 where ijk can be any cyclic permutation of 123. To compute K from Eq. 7.83, we first compute the derivatives ∂Ni = (yj − yk )/(2A) = bi /(2A), (7.85) ∂x ∂Ni = (xk − xj )/(2A) = ci /(2A), (7.86) ∂y where bi and ci are defined as bi = yj − yk , ci = x k − x j . (7.87) By a cyclic permutation of the indices, we obtain ∂Nj = (yk − yi )/(2A) = bj /(2A), ∂x (7.88) ∂Nj = (xi − xk )/(2A) = cj /(2A). ∂y (7.89) For the linear triangle, these derivatives are all constant and hence can be removed from the integral to yield 1 Kij = (bi bj + ci cj )|A|, (7.90) 4A2 where |A| is the area of the triangular element. Thus, the i, j contribution for one element is 1 (bi bj + ci cj ), (7.91) Kij = 4|A| where i and j each have the range 123, since there are three grid points in the element. However, bi and ci are computed from Eq. 7.87 using the shorthand notation that ijk is a cyclic permutation of 123; that is, ijk = 123, ijk = 231, or ijk = 312. Thus, K11 K22 K33 K12 K13 K23 = (b21 + c21 )/(4|A|) = [(y2 − y3 )2 + (x3 − x2 )2 ]/(4|A|), = (b22 + c22 )/(4|A|) = [(y3 − y1 )2 + (x1 − x3 )2 ]/(4|A|), = (b23 + c23 )/(4|A|) = [(y1 − y2 )2 + (x2 − x1 )2 ]/(4|A|), = (b1 b2 + c1 c2 )/(4|A|) = [(y2 − y3 )(y3 − y1 ) + (x3 − x2 )(x1 − x3 )]/(4|A|), = (b1 b3 + c1 c3 )/(4|A|) = [(y2 − y3 )(y1 − y2 ) + (x3 − x2 )(x2 − x1 )]/(4|A|), = (b2 b3 + c2 c3 )/(4|A|) = [(y3 − y1 )(y1 − y2 ) + (x1 − x3 )(x2 − x1 )]/(4|A|). (7.92) (7.93) (7.94) (7.95) (7.96) (7.97) Note that K depends only on differences in grid point coordinates. Thus, two elements that are geometrically congruent and differ only by a translation will have the same element matrix. To compute the right-hand side contributions as given in Eq. 7.81, consider first the contribution of the source term f . From Eq. 7.81, we calculate F1 , which is typical, as Z Z f f N1 dA = F1 = [(x2 y3 − x3 y2 ) + (y2 − y3 )x + (x3 − x2 )y] dA. (7.98) A 2A A 77 ✛ ❅ s ❅ 1 ❅ ❅ ❅ L ❅ ✲ s 2❅❅ ❅ Figure 61: Triangular Mesh at Boundary. We now consider the special case f = fˆ (a constant), and note that Z Z Z y dA = ȳ |A|, x dA = x̄ |A|, dA = |A|, A (7.99) A A where (x̄, ȳ) is the centroid of the triangle given by Eq. 7.47. Thus, F1 = fˆ [(x2 y3 − x3 y2 ) + (y2 − y3 )x̄ + (x3 − x2 )ȳ] |A|. 2A (7.100) From Eq. 7.48, the expression in brackets is given by 2A/3. Hence, 1 F1 = fˆ|A|, 3 (7.101) and similarly 1 (7.102) F1 = F2 = F3 = fˆ|A|. 3 That is, for a uniform f , 1/3 of the total element “load” is applied to each grid point. For nonuniform f , the integral could be computed using natural (or area) coordinates for the triangle [7]. It is also of interest to calculate, for the linear triangle, the right-hand side contributions for the second term of Eq. 7.81 for the uniform special case g = ĝ, where ĝ is a constant. For the second term of Eq. 7.81 to contribute, Point i must be on an element edge which lies on the boundary S2 . Since the integration involving g is on the boundary, the only shape functions needed are those which describe the interpolation of φ on the boundary. Thus, since the triangular shape functions are linear, the boundary shape functions are N1 = 1 − x x , N2 = , L L (7.103) where the subscripts 1 and 2 refer to the two element edge grid points on the boundary S2 (Fig. 61), x is a local coordinate along the element edge, and L is the length of the element edge. Thus, for Point i on the boundary, Z gNi dS. (7.104) Fi = − S2 Since the boundary shape functions are given by Eq. 7.103, Z L x ĝL 1− F1 = −ĝ dx = − , L 2 0 78 (7.105) L x ĝL . L 2 0 Thus, for the two points defining a triangle edge on the boundary S2 , F2 = −ĝ Z dx = − (7.106) ĝL . (7.107) 2 That is, for a uniform g, 1/2 of the total element “load” is applied to each grid point. To calculate H using Eq. 7.82, we first note that Points i, j must both be on the boundary for this matrix to contribute. Consider some triangular elements adjacent to a boundary, as shown in Fig. 61. Since the integration in Eq. 7.82 is on the boundary, the only shape functions needed are those which describe the interpolation of φ on the boundary. Thus, using the boundary shape functions of Eq. 7.103 in Eq. 7.82, for constant h = ĥ, Z L x 2 ĥL 1− H11 = ĥ dx = (7.108) L 3 0 Z L  2 ĥL x (7.109) dx = H22 = ĥ L 3 0 Z L xx ĥL 1− H12 = H21 = ĥ dx = . (7.110) L L 6 0 Hence, for an edge of a linear 3-node triangle,   ĥL 2 1 . (7.111) H= 1 2 6 F1 = F2 = − 7.6 Interpretation of Functional Now that the variational problem has been solved (i.e., the functional I has been minimized), we can evaluate I. We recall from Eq. 7.59 (with h = 0) that, for the two-dimensional Poisson’s equation,  Z  Z 1 I(φ) = gφ dS, (7.112) φ,k φ,k − f φ dA + 2 A S2 where φ = Ni φ i . (7.113) Since φi is independent of position, it follows that φ,k = Ni,k φi and  Z  Z 1 I(φ) = gNi φi dS, Ni,k φi Nj,k φj − f Ni φi dA + 2 A S2   Z Z Z 1 gNi dS φi f Ni dA − Ni,k Nj,k dA φj − = φi 2 S2 A A 1 = φi Kij φj − Fi φi (index notation) 2 1 = φT Kφ − φT F (matrix notation) 2 1 = φ · Kφ − φ · F (vector notation), 2 79 (7.114) (7.115) (7.116) (7.117) (7.118) where the last result has been written in index, matrix, and vector notations. The first term for I is a quadratic form. In solid mechanics, I corresponds to the total potential energy. The first term is the strain energy, and the second term is the potential of the applied loads. Since strain energy is zero only for zero strain (corresponding to rigid body motion), it follows that the stiffness matrix K is positive semi-definite. For well-posed problems (which allow no rigid body motion), K is positive definite. (By definition, a matrix K is positive definite if φT Kφ > 0 for all φ 6= 0.) Three consequences of positive-definiteness are 1. All eigenvalues of K are positive. 2. The matrix is non-singular. 3. Gaussian elimination can be performed without pivoting. 7.7 Stiffness in Elasticity in Terms of Shape Functions In §4.10 (p. 45), the Principle of Virtual Work was used to obtain the element stiffness matrix for an elastic finite element as (Eq. 4.67) Z CT DC dV, (7.119) K= V where D is the symmetric matrix of material constants relating stress and strain, and C is the matrix converting grid point displacements to strain: ε = Cu. (7.120) For the constant strain triangle (CST) in 2-D, for example, the fundamental unknowns u and v (the x and y components of displacement) can both be expressed in terms of the three shape functions defined in Eq. 7.38: u = Ni u i , v = Ni v i , (7.121) where the summation convention is used. In general, in 2-D, the strains can be expressed as       Ni,x ui u,x     εxx   = = ε= Ni,y vi v,y ε     yy   Ni,y ui + Ni,x vi u,y + v,x γxy   u   1           v1     N1,x 0 N2,x 0 N3,x 0   u 2   = . (7.122) 0 N1,y 0 N2,y 0 N3,y  v2      N1,y N1,x N2,y N2,x N3,y N3,x    u3        v3 80 φ ✻ ♣♣♣ ✲ x a b Figure 62: Discontinuous Function. Thus, from Eq. 7.120, C in Eq. 7.119 is a matrix of shape function derivatives:  N1,x 0 N2,x 0 N3,x 0 C =  0 N1,y 0 N2,y 0 N3,y  . N1,y N1,x N2,y N2,x N3,y N3,x  (7.123) In general, C would have as many rows as there are independent components of strain (3 in 2-D and 6 in 3-D) and as many columns as there are DOF in the element. 7.8 Element Compatibility In the variational approach to the finite element method, an integral was minimized. It was also assumed that the integral evaluated over some domain was equal to the sum of the integrals over the elements. We wish to investigate briefly the conditions which must hold for this assumption to be valid. That is, what condition is necessary for the integral over a domain to be equal to the sum of the integrals over the elements? Consider the one-dimensional integral Z b φ(x) dx. (7.124) I= a For the integral I to be well-defined, simple jump discontinuities in φ are allowed, as illustrated in Fig. 62. Singularities in φ, on the other hand, will not be allowed, since some singularities cannot be integrated. Thus, we conclude that, for any functional of interest in finite element analysis, the integrand may be discontinuous, but we do not allow singularities in the integrand. Now consider the one-dimensional integral Z b dφ(x) I= dx. (7.125) dx a Since the integrand dφ(x)/dx may be discontinuous, φ(x) must be continuous, but with kinks (slope discontinuities). In the integral I= Z b a d2 φ(x) dx, dx2 (7.126) the integrand φ′′ may have simple discontinuities, in which case φ′ is continuous with kinks, and φ is smooth (i.e., the slope is continuous). 81 2✘✘✘✘✄ ✄❅ ✄ ✄ ❅rP ✄ ✄ ✘❅ ✘✄ 1 ✄✘✘ Figure 63: Compatibility at an Element Boundary. Thus, we conclude that the smoothness required of φ depends on the highest order derivative of φ appearing in the integrand. If φ′′ is in the integrand, φ′ must be continuous. If φ′ is in the integrand, φ must be continuous. Therefore, in general, we conclude that, at element interfaces, the field variable φ and any of its partial derivatives up to one order less than the highest order derivative appearing in I(φ) must be continuous. This requirement is referred to as the compatibility requirement or the conforming requirement. For example, consider the Poisson equation in 2-D, ∂ 2φ ∂ 2φ + + f = 0, ∂x2 ∂y 2 for which the functional to be minimized is ) Z ( "  2  2 # 1 ∂φ ∂φ I(φ) = − f φ dA. + 2 ∂x ∂y A (7.127) (7.128) Since the highest derivative in I is first order, we conclude that, for I to be well-defined, φ must be continuous in the finite element approximation. For the 3-node triangular element already formulated, the shape function is linear, which implies that, in Fig. 63, given φ1 and φ2 , φ varies linearly between φ1 and φ2 along the line 1–2. That is, φ is the same at the mid-point P for both elements; otherwise, I(φ) might be infinite, since there could be a gap or overlap in the model along the line 1–2. Note also that the Poisson equation is a second order equation, but φ need only be continuous in the finite element approximation. That is, the first derivatives φ,x and φ,y may have simple discontinuities, and the second derivatives φ,xx and φ,yy that appear in the partial differential equation may not even exist at the element interfaces. Thus, one of the strengths of the variational approach is that I(φ) involves derivatives of lower order than in the original PDE. In elasticity, the functional I(φ) has the physical interpretation of total potential energy, including strain energy. A nonconforming element would result in a discontinuous displacement at the element boundaries (i.e., a gap or overlap), which would correspond to infinite strain energy. However, note that having displacement continuous implies that the displacement gradients (which are proportional to the stresses) are discontinuous at the element boundaries. This property is one of the fundamental characteristics of an approximate numerical solution. If all quantities of interest (e.g., displacements and stresses) were continuous, the solution would be an exact solution rather than an approximate solution. Thus, the approximation inherent in a displacement-based finite element method is that the displacements are continuous, and the stresses (displacement gradients) are discontinuous at element boundaries. In fluid mechanics, a discontinuous φ (which is not allowed) corresponds to a singularity in the velocity field. 82 z ✻ v ✟ ✟✟ ✟ ❍ ❍ ❍❍ x ✠ ✯✻ ✟✟ R = v − u (error) ✲ y ❍ ❥ ❍ u Figure 64: A Vector Analogy for Galerkin’s Method. 7.9 Method of Weighted Residuals (Galerkin’s Method) Here we discuss an alternative to the use of a variational principle when the functional is either unknown or does not exist (e.g., nonlinear equations). We consider the Poisson equation ∇2 φ + f = 0. (7.129) ∇2 φ̃ + f = R 6= 0, (7.130) For an approximate solution φ̃, where R is referred to as the residual or error. The best approximate solution will be one which in some sense minimizes R at all points of the domain. We note that, if R = 0 in the domain, Z RW dV = 0, (7.131) V where W is any function of the spatial coordinates. W is referred to as a weighting function. With n DOF in the domain, n functions W can be chosen: Z RWi dV = 0, i = 1, 2, . . . , n. (7.132) V This approach is called the method of weighted residuals. Various choices of Wi are possible. When Wi = Ni (the shape functions), the process is called Galerkin’s method. The motivation for using shape functions Ni as the weighting functions is that we want the residual (the error) orthogonal to the shape functions. In the finite element approximation, we are trying to approximate an infinite DOF problem (the PDE) with a finite number of DOF (the finite element model). Consider an analogous problem in vector analysis, where we want to approximate a vector v in 3-D with another vector in 2-D. That is, we are attempting to approximate v with a lesser number of DOF, as shown in Fig. 64. The “best” 2-D approximation to the 3-D vector v is the projection u in the plane. The error in this approximation is R = v − u, (7.133) which is orthogonal to the xy-plane. That is, the error R is orthogonal to the basis vectors ex and ey (the vectors used to approximate v): R · ex = 0 and R · ey = 0. 83 (7.134) In the finite element problem, the approximating functions are the shape functions Ni . The counterpart to the dot product is the integral Z RNi dV = 0. (7.135) V That is, the residual R is orthogonal to its approximating functions, the shape functions. The integral in Eq. 7.135 must hold over the entire domain V or any portion of the domain, e.g., an element. Thus, for Poisson’s equation, for one element, Z (7.136) 0 = (∇2 φ + f )Ni dV ZV = (φ,kk + f )Ni dV (7.137) V Z Z f Ni dV. (7.138) = [(φ,k Ni ),k − φ,k Ni,k ] dV + V V The first term is converted to a surface integral using the divergence theorem, Eq. 7.10, to obtain I Z Z 0= φ,k Ni nk dS − φ,k Ni,k dV + f Ni dV, (7.139) S V V where φ,k nk = ∇φ · n = ∂φ ∂n (7.140) and φ,k = (Nj φj ),k = Nj,k φj . (7.141) Hence, for each i, Z Z ∂φ f Ni dV 0= Nj,k φj Ni,k dV + Ni dS − S ∂n V V  Z  Z I ∂φ f Ni dV + Ni,k Nj,k dV φj + Ni dS =− S ∂n V V = −Kij φj + Fi . I (7.142) (7.143) (7.144) Thus, in matrix notation, where Kij = Fi = Z Kφ = F, (7.145) Z (7.146) Ni,k Nj,k dV V f Ni dV + V I S ∂φ Ni dS. ∂n (7.147) From Eq. 7.55, ∂φ/∂n is specified on S2 and unknown a priori on S1 , where φ = φ0 is specified. On S1 , ∂φ/∂n is the “reaction” to the specified φ. At points where φ is specified, the Dirichlet boundary conditions are handled like displacement boundary conditions in structural problems. 84 Galerkin’s method thus results in algebraic equations identical to those derived from a variational principle. However, Galerkin’s method is more general, since sometimes a variational principle may not exist for a given problem. When a principle does exist, the two approaches yield the same results. When the variational principle does not exist or is unknown, Galerkin’s method can still be used to derive a finite element model. 8 Potential Fluid Flow With Finite Elements In potential flow, the fluid is assumed to be inviscid and incompressible. Since there are no shearing stresses in the fluid, the fluid slips tangentially along boundaries. This mathematical model of fluid behavior is useful for some situations. Define a velocity potential φ such that velocity v = ∇φ. That is, in 3-D, ∂φ ∂φ ∂φ vx = , vy = , vz = . (8.1) ∂x ∂y ∂z It can be shown that, within the domain occupied by the fluid, ∇2 φ = 0. (8.2) Various boundary conditions are of interest. At a fixed boundary, where the normal velocity vanishes, ∂φ vn = v · n = ∇φ · n = = 0, (8.3) ∂n where n is the unit outward normal on the boundary. On a boundary where the velocity is specified, ∂φ = v̂n . (8.4) ∂n We will see later in the discussion of symmetry that, at a plane of symmetry for the potential φ, ∂φ = 0, (8.5) ∂n where n is the unit normal to the plane. At a plane of antisymmetry, φ = 0. The boundary value problem for flow around a solid body is illustrated in Fig. 65. In this example, far away from the body, where the velocity is known, ∂φ = v∞ , (8.6) vx = ∂x which is specified. As is, the problem posed in Fig. 65 is not a well-posed problem, because only conditions on the derivative ∂φ/∂n are specified. Thus, for any solution φ, φ + c is also a solution for any constant c. Therefore, for uniqueness, we must specify φ somewhere in the domain. Thus, the potential flow boundary value problem is that the velocity potential φ satisfies   ∇2 φ = 0 in V   φ = φ̂ on S1 (8.7)  ∂φ   = v̂n on S2 , ∂n where v = ∇φ in V . 85 ∂φ ∂φ = =0 ∂n ∂y n✛ n ✻ ∇2 φ = 0 ✒ ∂φ ∂φ =− = −v∞ ∂n ∂x ✲ n n ❅ ■ ❅ ∂φ ∂φ =− =0 ∂n ∂y ∂φ =0 ∂n ∂φ ∂φ = = v∞ ∂n ∂x n ❄ Figure 65: Potential Flow Around Solid Body. 8.1 Finite Element Model A finite element model of the potential flow problem results in the equation Kφ = F, where the contributions to K and F for each element are  Z  Z ∂Ni ∂Nj ∂Ni ∂Nj Ni,k Nj,k dA = dA, Kij = + ∂x ∂x ∂y ∂y A A Z v̂n Ni dS, Fi = (8.8) (8.9) (8.10) S2 where v̂n is the specified velocity on S2 in the outward normal direction, and Ni is the shape function for the ith grid point in the element. Once the velocity potential φ is known, the pressure can be found using the steady-state Bernoulli equation 1 2 p v + gy + = c = constant, (8.11) 2 ρ where v is the velocity magnitude given by  2  2 ∂φ ∂φ 2 + , (8.12) v = ∂x ∂y gy is the (frequently ignored) body force potential, g is the acceleration due to gravity, y is the height above some reference plane, p is pressure, and ρ is the fluid density. The constant c is evaluated using a location where v is known (e.g., v∞ ). For example, at infinity, if we ignore gy and pick p = 0 (ambient), 1 2 c = v∞ , (8.13) 2 86 ✲ ✬✩ ✲ y ✻ ✲ ✲ ✲ ✲ ✲ x ✫✪ Figure 66: Streamlines Around Circular Cylinder. y✻ rP ❍ ✤✜ ❥ ❍ ✲x ✣✢ ✯ ✟ r ✟ ′ P Figure 67: Symmetry With Respect to y = 0. and 8.2 p 1 2 1 2 + v = v∞ . ρ 2 2 (8.14) Application of Symmetry Consider the 2-D potential flow around a circular cylinder, as shown in Fig. 66. The velocity field is symmetric with respect to y = 0. For example, in Fig. 67, the velocity vectors at P and its image P ′ are mirror images of each other. As P and P ′ get close to the axis y = 0, the velocities at P and P ′ must converge to each other, since P and P ′ are the same point in the plane y = 0. Thus, ∂φ = 0 for y = 0. (8.15) vy = ∂y The y-direction in this case is the normal to the symmetry plane y = 0. Thus, in general, we conclude that, for points in a symmetry plane with normal n, ∂φ = 0. (8.16) ∂n Note from Fig. 65 that the specified normal velocities at the two x extremes are of opposite signs. Thus, from Eq. 8.10, the right-hand side “loads” in Eq. 8.8 are equal in magnitude and opposite in sign for the left and right boundaries, and the velocity field is antisymmetric with respect to the plane x = 0, as shown in Fig. 68. That is, the velocity vectors at P and P ′ can be transformed into each other by a reflection and a negation of sign. Thus, (vx )P = (vx )P ′ . (8.17) If we change the direction of flow (i.e., make it right to left), then (φP )flow right = (φP ′ )flow left . 87 (8.18) y ✻ P✟ P ✯ ❍r❍ r ✟✤✜ ❥ ′ ✲x ✣✢ Figure 68: Antisymmetry With Respect to x = 0. ∂φ =0 ∂n Antisym: φ=0 ∇2 φ = 0 ✒ a ∂φ = v∞ ∂n ∂φ =0 ∂n r Sym: ∂φ =0 ∂n Figure 69: Boundary Value Problem for Flow Around Circular Cylinder. However, changing the direction of flow also means that F in Eq. 8.8 becomes −F, since the only nonzero contributions to F occur at the left and right boundaries. However, if the sign of F changes, the sign of the solution also changes, i.e., (φP ′ )flow right = −(φP ′ )flow left . (8.19) Combining the last two equations yields (φP )flow right = −(φP ′ )flow right . (8.20) If we now let P and P ′ converge to the plane x = 0 and become the same point, we obtain (φ)x=0 = −(φ)x=0 (8.21) or (φ)x=0 = 0. Thus, in general, we conclude that, for points in a symmetry plane with normal n for which the solution is antisymmetric, φ = 0. The key to recognizing antisymmetry is to have a symmetric geometry with an antisymmetric “loading” (RHS). For example, for 2-D potential flow around a circular cylinder, Fig. 66, the boundary value problem is shown in Fig. 69. This problem has two planes of geometric symmetry, y = 0 and x = 0, and can be solved using a one-quarter model. Since φ is specified on x = 0, this problem is well-posed. The boundary condition ∂φ/∂n = 0 is the natural boundary condition (the condition that results if the boundary is left free). 88 y ✻ ❄ ❃ ✚ ✚ ✚ ✻ ❃ ✚ ✚ ✚ ✻ η = deflection of ✚ ✚ free surface disturbed free undisturbed free surface surface g d = depth ❄ ❄ Figure 70: The Free Surface Problem. 8.3 Free Surface Flows Consider an inviscid, incompressible fluid in an irrotational flow field with a free surface, as shown in Fig. 70. The equations of motion and continuity reduce to ∇2 φ = 0, (8.22) where φ is the velocity potential, and the velocity is v = ∇φ. (8.23) The pressure p can be determined from the time-dependent Bernoulli equation − ∂φ 1 2 p = + v + gy, ρ ∂t 2 (8.24) where ρ is the fluid density, g is the acceleration due to gravity, and y is the vertical coordinate. If we let η denote the deflection of the free surface, the vertical velocity on the free surface is ∂φ ∂η = on y = 0. (8.25) ∂y ∂t If we assume small wave height (i.e., η is small compared to the depth d), the velocity v on the free surface is also small, and we can ignore the velocity term in Eq. 8.24. If we also take the pressure p = 0 on the free surface, Bernoulli’s equation implies ∂φ + gη = 0 on y = 0. (8.26) ∂t This equation can be viewed as an equation for the surface elevation η given φ. We can then eliminate η from the last two equations by differentiating Eq. 8.26: 0= ∂ 2φ ∂η ∂φ ∂ 2φ + g +g . = 2 2 ∂t ∂t ∂t ∂y (8.27) Hence, on the free surface y = 0, ∂φ 1 ∂ 2φ . =− ∂y g ∂t2 This equation is referred to as the linearized free surface boundary condition. 89 (8.28) Im ✻ ✶A ✏✏ |A| ✏✏✏ ✏ ✏ ✏✏ ✏✏ ✏✏ ✏ ✏ θ ✲ Re Figure 71: The Complex Amplitude. 8.4 Use of Complex Numbers and Phasors in Wave Problems The wave maker problem considered in the next section will involve a forcing function which is sinusoidal in time (i.e., time-harmonic). It is common in engineering analysis to represent time-harmonic signals using complex numbers, since amplitude and phase information can be included in a single complex number. Such an approach is used with A.C. circuits, steady-state acoustics, and mechanical vibrations. Consider a sine wave φ(t) of amplitude Â, circular frequency ω, and phase θ: φ(t) =  cos(ωt + θ), (8.29) where all quantities in this equation are real, and  can be taken as positive. Using complex notation,  i i h h (8.30) φ(t) = Re Âei(ωt+θ) = Re Âeiθ eiωt , √ where i = −1. If we define the complex amplitude A = Âeiθ , (8.31)
φ(t) = Re Aeiωt , (8.32) |A| = Âeiθ = Â eiθ = Â, (8.33) then where the magnitude of the complex amplitude is given by which is the actual amplitude, and arg(A) = θ, (8.34) which is the actual phase angle. The complex amplitude A is thus a complex number which embodies both the amplitude and the phase of the sinusoidal signal, as shown in Fig. 71. The directed vector in the complex plane is called a phasor by electrical engineers. It is common practice, when dealing with these sinsusoidal functions, to drop the “Re”, and agree that it is only the real part which is of interest. Thus, we write φ(t) = Aeiωt (8.35) with the understanding that it is the real part of this signal which is of interest. In this equation, A is the complex amplitude. Two sinusoids of the same frequency add just like vectors in geometry. For example, consider the sum Â1 cos(ωt + θ1 ) + Â2 cos(ωt + θ2 ). 90 Im ✻ ✒ A1 + A2 A2 ✂ ✂ ✂ ✂ ✂✍ ✂ ✂ ✏✏✏ ✂ ✏ ✏ ✶ A1 ✏✏✏ ✲ Re Figure 72: Phasor Addition. y ✻ 1 ∂φ = − φ̈ (linearized free surface condition) ∂n g ✲ ✻ x ∇2 φ = 0 ∂φ = vn = v0 cos ωt (specified) ∂n d ✲ ∞ n✛ ∂φ = 0 (rigid bottom) ∂n ❄ Figure 73: 2-D Wave Maker: Time Domain. In terms of complex arithmetic, A1 eiωt + A2 eiωt = (A1 + A2 )eiωt , (8.36) where A1 and A2 are complex amplitudes given by A1 = Â1 eiθ1 , A2 = Â2 eiθ2 . (8.37) This addition, referred to as phasor addition, is illustrated in Fig. 72. 8.5 2-D Wave Maker Consider a semi-infinite body of water with a wall oscillating in simple harmonic motion, as shown in Fig. 73. In the time domain, the problem is 91     ∇2 φ = 0       ∂φ   = v0 cos ωt on x = 0 ∂n  ∂φ   = 0 on y = −d   ∂n     1 ∂φ   = − φ̈ on y = 0, ∂n g (8.38) where dots denote differentiation with respect to time, and an additional boundary condition is needed for large x at the location where the model is terminated. For this problem, the excitation frequency ω is specified. The solution of this problem is a function φ(x, y, t). The forcing function is the oscillating wall. We first write Eq. 8.38b in the form where i = √
∂φ = v0 cos ωt = Re v0 eiωt , ∂n (8.39) −1, and v0 is real. We therefore look for solutions in the form φ(x, y, t) = φ0 (x, y)eiωt , where φ0 (x, y) is the complex amplitude. Eq. 8.38 then becomes     ∇2 φ0 = 0       ∂φ0   = v0 on x = 0 ∂n  ∂φ0   = 0 on y = −d   ∂n     ω2 ∂φ0   = φ0 on y = 0. ∂n g (8.40) (8.41) It can be shown that, for large x, ∂φ0 = −iαφ0 , ∂x (8.42) ω2 = α tanh(αd), g (8.43) where α is the positive solution of and ω is the fixed excitation frequency. The graphical solution of this equation is shown in Fig. 74. Thus, for a finite element solution to the 2-D wave maker problem, we truncate the domain “sufficiently far” from the wall, and impose a boundary condition, Eq. 8.42, referred to as a radiation boundary condition which accounts approximately for the fact that the 92 y ✻ ✲ α Figure 74: Graphical Solution of ω 2 /α = g tanh(αd). ω2 ∂φ0 = φ0 ∂n g ✻ ∂φ0 = v0 ∂n d ✛ ∂φ0 = −iαφ0 ∂n ∇2 φ0 = 0 ✲ W ❄ ∂φ0 =0 ∂n Figure 75: 2-D Wave Maker: Frequency Domain. fluid extends to infinity. If the radiation boundary is located at x = W , the boundary value problem in the frequency domain becomes     ∇2 φ 0 = 0       ∂φ0   = v0 on x = 0 (oscillating wall)   ∂n   ∂φ0 (8.44) = 0 on y = −d (rigid bottom)  ∂n     ω2 ∂φ0   = φ0 on y = 0 (linearized free surface condition)   ∂n g       ∂φ0 = −iαφ0 on x = W (radiation condition), ∂x as summarized in Fig. 75. Note the similarity between this radiation condition and the nonreflecting boundary condition for the Helmholtz equation, Eq. 3.52. 93 8.6 Linear Triangle Matrices for 2-D Wave Maker Problem The boundary value problem defined in Eq. 8.44 has two boundaries where ∂φ/∂n is specified and two boundaries on which ∂φ/∂n is proportional to the unknown φ. Thus, this boundary value problem is a special case of Eq. 7.55 (p. 73), so that the finite element formulas derived in §7.5 are all applicable. Note that the function g appearing in Eq. 7.55 is not the acceleration due to gravity appearing in the formulation of the free surface flow problem. The matrix system for the wave maker problem is therefore, from Eq. 7.79, (K + H)φ = F, (8.45) where, for each element, K, H, and F are given by Eqs. 7.80–7.82. Thus, from Eq. 7.91, Kij = 1 (bi bj + ci cj ), 4|A| (8.46) where i and j each have the range 123, and bi = y j − y k , c i = xk − xj . (8.47) For bi and ci , the symbols ijk refer to the three nodes 123 in a cyclic permutation. For example, if j = 1, then k = 2 and i = 3. Thus, K11 K22 K33 K12 K13 K23 = (b21 + c21 )/(4|A|), = (b22 + c22 )/(4|A|), = (b23 + c23 )/(4|A|), = K21 = (b1 b2 + c1 c2 )/(4|A|), = K31 = (b1 b3 + c1 c3 )/(4|A|), = K32 = (b2 b3 + c2 c3 )/(4|A|). (8.48) (8.49) (8.50) (8.51) (8.52) (8.53) H is calculated using Eq. 7.111, where ĥ = −ω 2 /g on the free surface, and ĥ = iα on the radiation boundary. Thus, for triangular elements adjacent to the free surface,   ω2L 2 1 HFS = − (8.54) 1 2 6g for each free surface edge. For elements adjacent to the radiation boundary,   iαL 2 1 HRB = 1 2 6 (8.55) for each radiation boundary edge. Note that, since H is purely imaginary on the radiation boundary, the coefficient matrix K+H in Eq. 8.45 is complex, and the solution φ is complex. The solution of free surface problems thus requires either the use of complex arithmetic or separating the matrix system into real and imaginary parts. The right-hand side F is calculated using Eq. 7.107, where ĝ = −v0 on the oscillating wall. Thus, for two points on an element edge on the oscillating wall, F1 = F2 = 94 v0 L . 2 (8.56) ✲ k ✁✁❆❆✁✁❆❆✁✁ c u(t) ✲ m ❡ f (t) ❡ Figure 76: Single DOF Mass-Spring-Dashpot System. The solution vector φ obtained from Eq. 8.45 is the complex amplitude of the velocity potential. The time-dependent velocity potential is given by
(8.57) Re φeiωt = Re [(φR + iφI )(cos ωt + i sin ωt)] = φR cos ωt − φI sin ωt, where φR and φI are the real and imaginary parts of the complex amplitude. It is this function which is displayed in computer animations of the time-dependent response of the velocity potential. 8.7 Mechanical Analogy for the Free Surface Problem Consider the single DOF spring-mass-dashpot system shown in Fig. 76. The application of Newton’s second law of motion (F=ma) to this system yields the differential equation of motion mü + cu̇ + ku = f (t), (8.58) where m is mass, c is the viscous dashpot constant, k is the spring stiffness, u is the displacement from the equilibrium, f is the applied force, and dots denote differentiation with respect to the time t. For a sinusoidal force, f (t) = f0 eiωt , (8.59) where ω is the excitation frequency, and f0 is the complex amplitude of the force. The displacement solution is also sinusoidal: u(t) = u0 eiωt , (8.60) where u0 is the complex amplitude of the displacement response. If we substitute the last two equations into the differential equation, we obtain − ω 2 mu0 eiωt + iωcu0 eiωt + ku0 eiωt = f0 eiωt (8.61) (−ω 2 m + iωc + k)u0 = f0 . (8.62) or We make two observations from this last equation: 1. The inertia force is proportional to ω 2 and 180◦ out of phase with respect to the elastic force. 95 2. The viscous damping force is proportional to ω and leads the elastic force by 90◦ . Thus, in the free surface problem, we could interpret the free surface matrix HFS as an inertial effect (in a mechanical analogy) with a “surface mass matrix” M given by 1 L M= HFS = 2 −ω 6g  2 1 1 2  , (8.63) where the diagonal “masses” are positive. Similarly, we could interpret the radiation boundary matrix HRB as a “damping” effect with the “boundary damping matrix” B given by   1 αL 2 1 B = HRB = , (8.64) iω 6ω 1 2 where the diagonal dampers are positive. This “damping” matrix is frequency-dependent. The free surface problem is a degenerate equivalent to the mechanical problem, since the mass M occurs only on the free surface rather than at every point in the domain. In fact, the ideal fluid, for which ∇2 φ = 0, behaves like a degenerate mechanical system, because the ideal fluid possesses the counterpart to the elastic forces but not the inertial forces. This degeneracy is a consequence of the incompressibility of the ideal fluid. A compressible fluid (such as occurs in acoustics) has the analogous mass effects everywhere. 96 Bibliography [1] K.J. Bathe. Finite Element Procedures. Prentice-Hall, Inc., Englewood Cliffs, NJ, 1996. [2] J.W. Brown and R.V. Churchill. Fourier Series and Boundary Value Problems. McGraw-Hill, Inc., New York, seventh edition, 2006. [3] G.R. Buchanan. Schaum’s Outline of Theory and Problems of Finite Element Analysis. Schaum’s Outline Series. McGraw-Hill, Inc., New York, 1995. [4] D.S. Burnett. Finite Element Analysis: From Concepts to Applications. Addison-Wesley Publishing Company, Inc., Reading, Mass., 1987. [5] R.W. Clough and J. Penzien. Dynamics of Structures. McGraw-Hill, Inc., New York, second edition, 1993. [6] R.D. Cook, D.S. Malkus, M.E. Plesha, and R.J. Witt. Concepts and Applications of Finite Element Analysis. John Wiley and Sons, Inc., New York, fourth edition, 2001. [7] K.H. Huebner, D.L. Dewhirst, D.E. Smith, and T.G. Byrom. The Finite Element Method for Engineers. John Wiley and Sons, Inc., New York, fourth edition, 2001. [8] J.H. Mathews. Numerical Methods for Mathematics, Science, and Engineering. Prentice-Hall, Inc., Englewood Cliffs, NJ, second edition, 1992. [9] K.W. Morton and D.F. Mayers. Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge, 1994. [10] J.S. Przemieniecki. Theory of Matrix Structural Analysis. McGraw-Hill, Inc., New York, 1968 (also, Dover, 1985). [11] J.N. Reddy. An Introduction to the Finite Element Method. McGraw-Hill, Inc., New York, third edition, 2006. [12] F. Scheid. Schaum’s Outline of Theory and Problems of Numerical Analysis. Schaum’s Outline Series. McGraw-Hill, Inc., New York, second edition, 1989. [13] G.D. Smith. Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press, Oxford, England, third edition, 1985. [14] I.M. Smith and D.V. Griffiths. Programming the Finite Element Method. John Wiley and Sons, Inc., New York, fourth edition, 2004. [15] M.R. Spiegel. Schaum’s Outline of Theory and Problems of Advanced Mathematics for Engineers and Scientists. Schaum’s Outline Series. McGraw-Hill, Inc., New York, 1971. [16] J.S. Vandergraft. Introduction to Numerical Calculations. Academic Press, Inc., New York, second edition, 1983. 97 [17] O.C. Zienkiewicz and R.L. Taylor. The Finite Element Method for Solid and Structural Mechanics. Elsevier Butterworth-Heinemann, Oxford, England, sixth edition, 2005. [18] O.C. Zienkiewicz, R.L. Taylor, and P. Nithiarasu. The Finite Element Method for Fluid Dynamics. Elsevier Butterworth-Heinemann, Oxford, England, sixth edition, 2005. [19] O.C. Zienkiewicz, R.L. Taylor, and J.Z. Zhu. The Finite Element Method: Its Basis and Fundamentals. Elsevier Butterworth-Heinemann, Oxford, England, sixth edition, 2005. 98 Index acoustics, 13 domain of influence, 23 back-solving, 10 banded matrix, 33 beams in flexure, 44 Bernoulli equation, 86, 89 big O notation, 4 boundary conditions, 1 derivative, 20, 33 essential, 76 natural, 65, 76, 88 Neumann, 20 nonreflecting, 27 radiation, 92 boundary value problem(s), 1, 8 brachistochrone, 61 electrostatics, 12 equation(s) acoustics, 13 Bernoulli, 86, 89 classification, 14 elliptic, 14, 31 Euler-Lagrange, 59, 63 heat, 13 Helmholtz, 13 hyperbolic, 14, 21 Laplace, 11 mathematical physics, 11 nondimensional, 15 ordinary differential, 1 parabolic, 14, 16 partial differential, 11 Poisson, 12, 66, 73 potential, 11 sparse systems, 32 systems, 5 wave, 12, 22 error global, 3 local, 3 rounding, 3 truncation, 3 essential boundary condition, 76 Euler beam, 45 Euler’s method, 2 truncation error, 3 Euler-Lagrange equation, 59, 63 explicit method, 4, 16 calculus of variations, 57 CFL condition, 26 change of basis, 49 compatibility, 81 complex amplitude, 90 complex numbers, 90 conic sections, 14 constant strain triangle, 48 constrained minimization, 63 constraints, 36 continuum problems direct approach, 45 Courant condition, 26 Crank-Nicolson method, 10, 18 stencil, 20 CST, 48 cyclic permutation, 71 cycloid, 62 d’Alembert solution, 21, 22 del operator, 11 determinant of matrix, 52 discriminant, 14 displacement vector, 35 divergence theorem, 67 domain of dependence, 25 finite difference(s), 6 backward, 6 central, 6 forward, 2, 6 Neumann boundary conditions, 33 relaxation, 33 Fourier’s law, 13 99 free surface flows, 89 functional, 58 Galerkin’s method, 83 Gaussian elimination, 10, 80 gravitational potential, 11 heat conduction, 12 Helmholtz equation, 13 Hooke’s law, 48 implicit method, 4 incompressible fluid flow, 11 index notation, 67 initial conditions, 1 initial value problem(s), 1 interpretation of functional, 79 Kronecker delta, 50, 57, 68 Laplacian operator, 11 large spring method, 43 Leibnitz’s rule, 22 magnetostatics, 12 mass matrix, 35 mass-spring system, 1, 34 matrix assembly, 36 matrix partitioning, 42 mechanical analogy, 95 method of weighted residuals, 83 natural boundary condition, 65, 76, 88 Neumann boundary condition finite differences, 20, 33 Newton’s second law, 34 nonconforming element, 82 nondimensional form, 15 orthogonal coordinate transformation, 52 orthogonal matrix, 52 orthonormal basis, 50 phantom points, 20 phasors, 90 addition, 91 pin-jointed frame, 41 pivoting, 80 Poisson equation, 12 positive definite matrix, 80 potential energy, 80 potential fluid flow, 85 radiation boundary condition, 92 relaxation, 33 rod element, 38 rotation matrix, 52 rounding error, 3 Runge-Kutta methods, 4 separation of variables, 18 shape functions, 70 shooting methods, 10 solution procedure, 38 sparse system of equations, 32 speed of propagation, 12 stable solution, 18 stencil, 16, 20, 26 stiffness matrix, 35, 76 properties, 38 strain energy, 80 summation convention, 50, 67 symmetry, 87 Taylor’s theorem, 3, 4 tensors, 53 examples, 54 isotropic, 57 torsion, 12 transverse shear, 45 triangular element, 46 tridiagonal system, 8, 10 truncation error, 3 truss structure, 40 unit vectors, 50 unstable solution, 18 velocity potential, 11, 85 vibrations bar, 12 membrane, 12 string, 12 virtual work, 48 100 warping function, 12 wave equation, 12, 22 wave maker problem, 91 matrices, 94 wave speed, 23 101