Homework # 3 ME526/NSE526 Due: November 5 1. You are encouraged to work in a group of up to 2 students and submit one solution per group. 2. Your solution must be clearly legible. Illegible work may not be graded and returned without any points. Although not necessary, you may type your work. 3. All problems must be solved. However, all problems may not be graded. A random sample of problems will be selected for grading. 4. If you are required to write a computer program, attach your code with several comment statements on the code wherever possible. 1 1. The following scheme is proposed to solve y′ = f(y): y∗ = yn + γ1hf(yn) (1) y∗∗ = y∗ + γ2hf(y∗) + ω2hf(yn) (2) yn+1 = y ∗∗ + γ3hf(y∗∗) + ω3hf(y∗) (3) where γ1 = 8/15, γ2 = 5/12, γ3 = 3/4, ω2 = −17/60, and ω3 = −5/12 with h being the time step. (a) What is the order of accuracy of this method? Prove the order by applying the method to y′ = λy and obtaining a finite difference approximation that can be written in the yn+1 = σyn. (b) Draw a stability diagram in the (hλR, hλI) plane for this method applied to the model problem y′ = λy. What is the maximum step size for λ pure imaginary and λ negative real? 2. Consider motion of a simple pendulum (see Fig 1) consisting of a mass m attached to a string of length ` and immersed in a viscous fluid adding some damping resistance. The governing equation is given by θ′′ + cθ′ + g ` sin(θ) = 0, (4) where g is the gravitational acceleration, and θ the angle made by the string with the vertical axis, c is the damping coefficient. For small angles θ, the linearized equation of motion is: θ′′ + cθ′ + g ` θ = 0. (5) Let g = 9.81 m/sec2, and ` = 0.6 m, c = 4 sec−1. Assume that the pendulum starts from rest (i.e. θ′(t = 0) = 0) with θ(t = 0) = 100. Figure 1: Schematic of a simple pendulum of mass m and length ` immersed in a viscous fluid. (a) Evaluate the exact solution for the linearized system (equation 5). 2 (b) Consider the following numerical schemes to advance the solution over the time interval 0 ≤ t ≤ 6. For each scheme clearly indicate the system of discretized equations (to advance from t = tn to t = tn+1) along with the proper initial conditions to be used. i. The Forward Euler, ii. The Second Order Runge-Kutta, iii. The Fourth Order Runge-Kutta. (c) Perform linear stability analysis and discuss quantitatively and in detail the stability limits for each scheme. (d) Evaluate the solution to the linearized problem 1, θ(t), using all the schemes for 0 ≤ t ≤ 6. Use ∆t = 0.15, 0.5, 1 s. For each case, and on separate figures, plot the solution over the time period together with the exact solution to the linearized system2 Use subplot to plot all figures (that is total of 9 figures, 3 for each scheme and 3 different time-steps) on a single page. Clearly label all plots. (e) Consider the non-linear equation 4. For an undamped system (i.e. c = 0), solve the differential equations with θ(t = 0) = 600 with RK4. What steps have you taken to be certain of the accuracy of your results? Why should your results be believable? Plot the solution for 0 ≤ t ≤ 6. Compare your solution with exact solution of the linearized undamped governing equation (i.e. solution to equation 5 with c = 0). Comment on what is the effect of the non-linearity on the solution. For the numerical approach, find out any limit on the time-step to obtain stable solution for the non-linear problem using RK4 (This should be done by increasing the time-steps in small increments and finding a time-step beyond which unphysical solution is obtained. List down the time-step at which this happens). Compare this time-step obtained from numerical experiments to the time-step restriction provided by the linear-stability theory for the undamped problem. How do they compare? Comment on your findings. 3. Consider chemical reactions during food digestion in our bodies. An enzyme (E) combines with a substrate (S) to form a complex (ES). The complex can dissociate back to E and S or it can proceed to form product P. The time history of the reactions E + S k1−⇀↽− k2 ES k3→ E + P, (6) is gorverned by the following rate equations The above reaction are governed by the following 1You will have to write a computer program to do this with separate routines for each scheme. 2In your plot θ should be in degrees. Make sure you use separate line-styles or symbols to distinguish the plots, also clearly label x and y axes, and use the same plotting format for all figures (e.g. the axis range etc.). If your answer blows up, plot the answer up to a certain range of θ axis to allow comparison with other schemes. 3 rate equations: dCS dt = −k1CSCE + k2CES (7) dCE dt = −k1CSCE + (k2 + k3)CES (8) dCES dt = k1CSCE − (k2 + k3)CES (9) dCP dt = k3CES (10) where Cis are concentrations and k1, k2 and k3 are reaction rate constants given as k1 = 2× 103; k2 = 1× 10−3; k3 = 10.0. (11) Initially, CS(0) ∼ 1.0, CE(0) = 5×10−5, CES(0) = 0, and CP = 0. Using C = (C1, C2, C3, C4) = (CS, CE, CES, CP), Solve the above system numerically to steady state (let t = 2500s repre- sent steady state) using the following schemes. (a) RK4. Do you think ∆t = 0.0005 will give you a stable solution? Justify your answer based on Stability Analysis at t = 0. Use this time-step to obtain the steady state solution. Note the computing time required to obtain steady state. Compare steady state solution to the exact. Plot all concentrations on the same plot (use axes ranges from 10−10 to 1.5× 100 for concentrations and 10−5 to 3000 for time. Provide solutions in a log-log plot as well as semi-log plot (where the concentrations are on a log scale and time in a linear scale). (b) Setup the problem with a linearized trapezoidal scheme. To do this, first write down the finite-difference equations for the standard trapezoidal scheme. Since this is an implicit method, you will have non-linear terms for the unknown concentrations at time level tn+1. Apply the linearization technique to the unknown non-linear terms in your finite-difference approximation (fda) and re-write the fda indicating the time level at which each term is evaluated. Your final answer should be in the form: MY = NX+b where M, N are matrices, Y , X, and b are some vectors. Clearly indicate the elements of each matrix and vector and also the time level at which each element is evaluated. Would a time step of 0.001 provide an accurate solution at steady state? Solve and check. What advantages would such a linearized central differencing scheme would have over say the fourth-order RK method? (c) Compare the computing time required for the three approaches above. 4
欢迎咨询51作业君