The assignment consists of hand calculations and problems that need to be worked on through MATLAB. There are 2 problems, each with a variety of subproblems.
Problem 1: Numerical Differentiation (50 points) In this problem, we will use numerical differentiation and Euler-Bernoulli beam theory to approximate the elastic curve of a clean airplane wing under elliptical lift distribution. The wing has length L and, for simplicity, is assumed to be a uniform cantilevered beam as shown in Figure 1. Recall that the governing ODE for the deflection v of a homogeneous beam with constant cross-section is given by EIv(4)(x) = w(x) , (1) where • E is the beam’s Young’s modulus, • I is the area moment of inertia of the cross section, and • w is the distributed transverse load with positive direction as shown in the figure. Since the wing is fixed at the left end and free at the right, the deflection must satisfy the following boundary conditions: • Zero displacement at x = 0 =) v(0) = 0 , (2a) • Zero rotation at x = 0 =) v0(0) = 0 , (2b) • Zero bending moment at x = L =) v00(L) = 0 , (2c) • Zero shear force at x = L =) v000(L) = 0 . (2d) 1. (40 points) Problem Set-up and Implementation (a) (6 points) Assuming that v has as many derivatives as needed, use the data v(x � 2h) , v(x � h) , v(x) , v(x + h) , and v(x + 2h) to show that the 4th derivative of v at the point x can be approximated as v(4)(x) ⇡ v(x� 2h)� 4v(x� h) + 6v(x)� 4v(x+ h) + v(x+ 2h) h4 . (3) (b) (1 point) What is the order of the approximation? We set up a uniform grid along the beam, i.e. on [0 , L], by letting h = L n and xj = jh , j = 0 , 1 , . . . , n , where n+ 1 is the total number of grid points (including the endpoints). Furthermore, we let vj be the approximation of the deflection at the mesh point xj , i.e., vj ⇡ v(xj) , and let wj = w(xj) h4 EI . 2 (c) (2 points) Use the approximation of Equation (3) to evaluate Equation (1) at xj, thus relating wj and the approximations vj�2 , vj�1 , vj , vj+1 , vj+2 of the beam deflection at the grid points xj�2 , xj�1 , xj , xj+1 , xj+2 respectively. Note that the above result is problematic near the endpoints as it requires knowledge of data beyond the beam. For example, using the result of part (c) with j = 1 would require knowledge of v�1 which does not exist within the grid domain. We consider this below. (d) (1.5 points) For j = 1 , use the result of part (c) to obtain an equation relating v1 , v2 , v3 , v�1 , and w1 . Hint: Remember that v(0) = 0 from the boundary condition of Equation (2a). (e) (1 point) Recall the central difference formula for the 1st derivative: v0(x) ⇡ v(x+ h)� v(x� h) 2h . (4) State the order of the central difference scheme of Equation (4). (f) (2 points) Use Equation (4) and the boundary condition of Equation (2b) to express v�1 in terms of v1 . (g) (2 points) Use the results of parts (d) and (f) to obtain an equation relating v1 , v2 , v3 , and w1 . (h) (1.5 points) Next, for j = n � 1 , use the result of part (c) to obtain an equation relating vn�3 , vn�2 , vn�1 , vn , vn+1, and wn�1 . (i) (1 point) Likewise, recall the central difference formula for the 2nd derivative: v00(x) ⇡ v(x� h)� 2 v(x) + v(x+ h) h2 . (5) State the order of this central difference scheme. (j) (2 points) Use Equation (5) and the boundary condition of Equation (2c) to express vn+1 in terms of vn�1 and vn . (k) (2 points) Use the results of parts (h) and (j) to obtain an equation relating vn�3 , vn�2 , vn�1 , vn , and wn�1 . (l) (1.5 points) Finally, for j = n , use the result of part (c) to obtain an equation relating vn�2 , vn�1 , vn , vn+1 , vn+2, and wn . (m) (2 points) It can be shown that a second order approximation for the 3rd derivative is v000(x) ⇡ �v(x� 2h) + 2v(x� h)� 2v(x+ h) + v(x+ 2h) 2h3 . (6) Use Equation (6) and the boundary condition of Equation (2d) to express vn+2 in terms of vn�2 , vn�1 and vn+1 . 3 (n) (2.5 points) Insert your results of parts (j) and (m) into that of part (l) to obtain an equation relating vn�2 , vn�1 , vn , and wn . (o) (6 points) Combine the results of parts (c), (g), (k), and (n) to obtain the system of equations to be solved for the unknown deflections v1 , v2 , . . . , vn . Specifically, letting v = (v1 , v2 , . . . , vn)T and w = (w1 , w2 , . . . , wn)T , determine the matrix A 2 Rn⇥n such that Av = w . In what follows, use values for a beam made out of steel: E = 2.10⇥1011 N/m2, I = 9.1⇥10�6 m4, L = 5m, w(x) = w0 r 1� x2 L2 with w0 = 1284 N/m . (p) (6 points) Use Matlab® to solve for and plot the approximate deflection of the beam for n = 24 . Call your script Problem1.m. Hint: Remember to include the point x0 = 0 . On the same figure, also plot the exact solution given by v(x) = w0 3EI ( 1 16 " 1 15 (16L4 + 83L2x2 + 6x4) r 1� x2 L2 +x(3L3 + 4Lx2) arcsin ⇣x L ⌘i � ⇡L 8 x3 � L4 15 � 2. (10 points) Error Analysis For a given n , the infinity norm of the approximation error is computed as en = max 0jn ��vj � v(xj) �� = max ���v0 � v(x0) ��, ��v1 � v(x1) ��, . . . , ��vn � v(xn) �� , (7) where, recall, v(xj) is the exact solution at xj and vj is the approximate deflection at xj . If p is the order of the approximation, then it is expected that the infinity norm of the error would obey en = O(h p) = O(n�p) as n approaches infinity. In other words, as n ! 1 , we can write en ⇡ C · n �p for some constant C . (8) (a) (4 points) Using Equation (7), add to your script from Part 1(p) to tabulate the infinity norm of the approximation error for the following values of n : 24 , 25 , 26 , 27 , 28 , 29 , 210. (b) (3 points) Plot the result of part (a) on a log-log plot and use Matlab® ’s polyfit function to determine the slope of the best-fit line through the points (log n , log en) . (c) (3 points) Noticing from Equation (8) that log en ⇡ �p log n + logC , deduce the order p of the approximation. Is this what you expected in light of the order(s) of the approximations given in Equations (3), (4), (5). Why or why not? 4 L
AAACbHicbVHNSsNAEN7G//rXqjcRFoviqSQq6LHoxYMHBVOFNpTNdmKXbrJhd9JSQp/Aqz6cL+EzuKk9GOvAsh/f/H0zE6ZSGHTdz4qztLyyura+Ud3c2t7ZrdX32kZlmoPPlVT6JWQGpEjAR4ESXlINLA4lPIfD28L/PAJthEqecJJCELPXRESCM7TU432v1nCb7szoIvDmoEHm9tCrV/xuX/EshgS5ZMZ0PDfFIGcaBZcwrXYzAynjQ/YKHQsTFoMJ8pnSKT2xTJ9GStuXIJ2xvzNyFhsziUMbGTMcmL++gvzP18kwug5ykaQZQsJ/GkWZpKhoMTbtCw0c5cQCxrWwWikfMM042uWUuoRxaYYczMhm2orFZ4WPBQ5s2ZQyrdW4FFrIQqWkKVcIQyX75cXMBuA6yCGzn0hxag/h/V37ImifN72Lpvt42WjdzE+yTg7JMTkjHrkiLXJHHohPOAHyRt7JR+XLOXAOnaOfUKcyz9knJXNOvwEa7L76 v
AAACbHicbVHLTgIxFC3jC/EF6s6YNBKNKzKjJrokunEJiTwSmJhOuUBjZzpp70DIhC9wqx/nT/gNdoCFI96k6cm5r3PvDWIpDLruV8HZ2Nza3inulvb2Dw6PypXjtlGJ5tDiSirdDZgBKSJooUAJ3VgDCwMJneDtKfN3JqCNUNELzmLwQzaKxFBwhpZqTl7LVbfmLoyuA28FqmRljddKodUfKJ6EECGXzJie58bop0yj4BLmpX5iIGb8jY2gZ2HEQjB+ulA6p5eWGdCh0vZFSBfs74yUhcbMwsBGhgzH5q8vI//z9RIcPvipiOIEIeLLRsNEUlQ0G5sOhAaOcmYB41pYrZSPmWYc7XJyXYIwN0MKZmIzbcXss8KnAse2bEyZ1mqaC81koVLS5CsEgZKD/GIWA3Dtp5DYT8Q4t4fw/q59HbRvat5tzW3eVeuPq5MUyRm5INfEI/ekTp5Jg7QIJ0DeyQf5LHw7p86Zc74MdQqrnBOSM+fqB3A8vyQ= x
AAACbHicbVHNSsNAEN7E//rXqjcRFoviqSQq6FH04lHB2EIbZLOd2qWbbNidWEvoE3jVh/MlfAY3bQ+mdWDZj2/+vpmJUikMet634y4tr6yurW9UNre2d3artb1nozLNIeBKKt2KmAEpEghQoIRWqoHFkYRmNLgr/M030Eao5AlHKYQxe01ET3CGlnp8f6nWvYY3MboI/Bmok5k9vNScoNNVPIshQS6ZMW3fSzHMmUbBJYwrncxAyviAvULbwoTFYMJ8onRMTyzTpT2l7UuQTti/GTmLjRnFkY2MGfbNvK8g//O1M+xdh7lI0gwh4dNGvUxSVLQYm3aFBo5yZAHjWlitlPeZZhztckpdorg0Qw7mzWbaisVnhQ8F9m3ZlDKt1bAUWshCpaQpV4giJbvlxUwG4DrMIbOfSHFsD+HPr30RPJ83/IuG93hZv7mdnWSdHJJjckZ8ckVuyD15IAHhBMgH+SRfzo974B66R9NQ15nl7JOSuae/dEy/Jg== AAAB7HicbVBNSwMxEJ2tX7V+VT16CRahXsquFOqx6MVjBWsL7VKyabYNTbJLklXL0r/gxYOCePUHefPfmG33oK0PBh7vzTAzL4g508Z1v53C2vrG5lZxu7Szu7d/UD48utdRoghtk4hHqhtgTTmTtG2Y4bQbK4pFwGknmFxnfueBKs0ieWemMfUFHkkWMoJNJj1Wn84H5Ypbc+dAq8TLSQVytAblr/4wIomg0hCOte55bmz8FCvDCKezUj/RNMZkgke0Z6nEgmo/nd86Q2dWGaIwUrakQXP190SKhdZTEdhOgc1YL3uZ+J/XS0x46adMxomhkiwWhQlHJkLZ42jIFCWGTy3BRDF7KyJjrDAxNp6SDcFbfnmVdC5qXr3mebf1SvMqz6MIJ3AKVfCgAU24gRa0gcAYnuEV3hzhvDjvzseiteDkM8fwB87nDx1Bjhs= w(x) Figure 1: Simplified (clean) airplane wing subjected to an elliptical lift distribution. Problem 2: Partial Differential Equations (50 points) As we saw in class, the 1-D heat equation can be used to model the diffusion of a property of a substance (e.g. incompressible fluid) in a one-dimensional medium (rod, bar. etc.) If in addition we consider the transport of the substance in the medium, we obtain the 1-D advection-diffusion equation given by ut + c ux = Duxx , (9) where u(x, t) describes the state of the fluid (e.g.: temperature, density, energy, etc.), D > 0 is the diffusion coefficient (assumed to be constant), and c > 0 is the advection velocity (which we will also assume to be constant). Furthermore, the initial state of the fluid is given by u(x, 0) = f(x) for all x . (10) Let wi,j be the numerical approximation to the true solution u(xi, tj) at the point (xi, tj) on a grid of (spatial) length L , where, taking h and k as the spatial and temporal step sizes respectively, we have • xi = ih , for i = 0, 1, 2, 3, . . . ,m+ 1 , with x0 = 0 and xm+1 = (m+ 1)h = L corresponding to the left and right boundaries of the grid respectively, and • tj = jk , for j = 0, 1, 2, 3, . . . Because Equation (9) has spatial derivatives, we must also specify one boundary conditions for the variable x . Here, we will do so through the use of periodic boundary conditions. Specifically, we will assume that the grid repeats indefinitely so that we may write u(0, t) = u(L, t) for t � 0 . (11) Physically, this means that whatever flows out of the boundary x = L must flow back into the grid at the boundary x = 0 . 5 1. (17.5 points) Forward in Time – Central in Space (FTCS) Scheme In this part, we use the forward Euler in time and central difference in space for the 1st and 2nd derivatives to discretize Equation (9) over the grid. (a) (3 points) Using the FTCS finite difference scheme, show that wi,j satisfies wi,j+1 = wi,j � � 2 (wi+1,j � wi�1,j) + � (wi�1,j � 2wi,j + wi+1,j) , or, equivalently, wi,j+1 = ✓ � + � 2 ◆ wi�1,j + (1� 2�)wi,j + ✓ � � � 2 ◆ wi+1,j , (12) where � ⌘ c k h and � ⌘ Dk h2 . (b) (1.5 points) Draw the stencil for the FTCS scheme of Equation (12). Be sure to label it as done in class. Is the scheme implicit or explicit? Why? (c) (3 points) By considering i = 0, 1, 2, . . . ,m � 1,m,m + 1 , show that Equation (12) can be represented in matrix form as