Computational Fluid Dynamics Programming Project

Question Description

Help me study for my Python class. I’m stuck and don’t understand.

From Professor's email:

Hi class,

I will be grading Homework #2 on Sunday, March 29th. Since I won't look at the submissions until that time, I added a grace period of an extra day to the submission deadline on the course webpage.

1. Please make sure the file you upload contains plots of your results.
2. Attach all of your code.
3. In addition, to reviewing your results (plots). I also look at your code. I try to understand what you are doing and if something looks strange, I will try to run it. It is helpful for me if you add some documentation to your code.

I know many of you work together on homework. I encourage this, but I also want to make sure each of you individually understand the work you are submitting. With the code I have already published and with any code potentially shared amongst you during discussions of the homework, there exists is a gray area of what is plagiarism and what is not.

The key point here is that you need to show me in your own way that you understand how to solve these problems. Some things to remember:

1. Add comments in your code that tell me in your owns words what you are doing. I don't need paragraphs, but when I review code, I try to follow your line of thought. I want to know what your logic is.

2. Be an editor. Just like you would proof read a final term paper, you should also edit your code. Think to yourself, is that really the best way to do that? Can I improve this line of code to make it more clear to anyone who might use my code what it is doing? Do I really understand this bit of code, or is it just some snippet of code I copied from stackoverflow? Can I make it more concise? Should I use a function here instead?

3. Go beyond just generating the plots I asked for, show me you really understand the steps involved. One way to do this is to test each component of your code. It is a great habit to get into, and you probably already have done a decent amount of debugging why you solved the problem. Why not formalize it? For instance, you could test that the components of your banded matrix are indeed what you expected them to be or that your initial condition looks right.

Unformatted Attachment Preview

Computational Fluid Dynamics: Homework #2 Due on March 16, 2020 at 5:00 pm 1 Computational Fluid Dynamics : Homework #2 Problem 1 Problem 1 [5 Points] Solve the initial value problem for the model hyperbolic equation (the one-dimensional linear advection equation), described by the partial differential equation ∂u ∂u +c = 0, ∂t ∂x (1) using the five explicit methods provided below. Bonus points will be given for correct solutions obtained from the MacCormack and Warming-Beam methods. For each numerical method, use N = 41 points to span the interval 0 ≤ x ≤ Lx with a Lx = 2. Assume the mesh is uniform in the x-direction, where the location at each point is given by xi = (i − 1)∆x with ∆x = Lx /(N − 1) For each method use the stability condition, c∆t = 0.9 ∆x to compute the time step, ∆t. Use the value c = 1. Run each method for 10 time steps and compare the results to the analytical solution, u(x, t) = u0 (x − ct) For reference, the derivation for the analytical solution to the linear convection problem using the method of characteristics is provided below. Use the following initial condition, ( 1.0 x ≤ 0.5 u0 (x, 0) = 0.5 x > 0.5 Since the wave is moving to the right (c > 0), the boundary condition at x = 0 is fixed by the initial condition, u(0, t) = u0 (0, t). At the right boundary x = 2, always use the backward difference formula, (D− /∆x)·, for the spatial derivative to avoid needing to specify a ”ghost” node at that boundary. Explicit Numerical Methods • Explicit Backward [1 Point] D− · n u ∆x i ∆x Stability Criteria : c > 0, ∆t ≤ |c| Order of Accuracy : O (∆t, ∆x) Method : un+1 = uni − c∆t i (2) • Explicit Forward [1 Point] D+ · n u ∆x i ∆x Stability Criteria : c > 0, ∆t ≤ |c| Order of Accuracy : O (∆t, ∆x) Method : un+1 = uni − c∆t i (3) • Explicit Central [1 Point] Method : un+1 = uni − c∆t i Stability Criteria : Always unstable  Order of Accuracy : O ∆t, ∆x2 Problem 1 continued on next page. . . 2 D0 · n u ∆x i (4) Computational Fluid Dynamics : Homework #2 Problem 1 (continued) • Lax [1 Point] uni−1 + uni+1 D0 · n − c∆t u 2 ∆x i ∆x Stability Criteria : ∆t ≤ |c|  Order of Accuracy : O ∆t, ∆x2 Method : un+1 = i (5) • Lax-Wendroff [1 Point] Method : Stability Criteria : un+1 = uni − c∆t i ∆t ≤ D0 · n 1 2 2 D+ · D− · n u + c ∆t u ∆x i 2 ∆x ∆x i (6) ∆x |c| Order of Accuracy : O ∆t2 , ∆x2  • MacCormack [Bonus: 1 Point] D+ · n Method : un+1 u = uni − c∆t i ∆x i   1 n D− · n+1 n+1 n+1 ui = u + ui − c∆t u 2 i ∆x i ∆x Stability Criteria : ∆t ≤ |c|  2 Order of Accuracy : O ∆t , ∆x2 (7) • Warming-Beam [Bonus: 1 Point] c∆t D− · n ui 2 ∆x   D− · n+1/2 ∆x D− · n n un+1 = u − c∆t u u + i i i ∆x 2 ∆x i 2∆x Stability Criteria : ∆t ≤ |c|  Order of Accuracy : O ∆t2 , ∆x2 Method : n+1/2 ui = uni − (8) Difference Formulas Assuming a uniform mesh size of ∆x, the difference operators, D, used in the difference equations for the explicit numerical methods given in the previous section are defined as • Forward difference, first-order accurate φi+1 − φi D+ · φi = ∆x ∆x (9) • Backward difference, first-order accurate D− · φi − φi−1 φi = ∆x ∆x (10) • Central difference, second-order accurate D0 · φi+1 − φi−1 φi = ∆x 2∆x Problem 1 continued on next page. . . 3 (11) Computational Fluid Dynamics : Homework #2 Problem 1 (continued) Analytical Solution to a Hyperbolic PDE The solution to Eq. 1 can be determined from the initial condition u(x, 0) = u0 (x) and the boundary condition, u(0, t) = u0 (0, t) = ul . If c ≥ 0, then we only need the boundary condition at x = 0. Also, note that the boundary condition is informed by the initial condition; hence we are providing solutions to the initial value problem. The method of characteristics can be used to find an analytical solution to the above equation. Along a characteristic line, we have dx = c. (12) dt As a result, along a particular curve x = x(t), we have i ∂u dx ∂u dh ∂u ∂u u(x(t), t) = + = +c . (13) dt ∂t dt ∂x ∂t ∂x The value u is therefore constant along the characteristics (the material derivative du/dt is zero, thus u is constant along lines x − t lines). We also know that these characteristic lines are straight, since the curvature (second-derivative) of the characteristics is zero,   d dx dc d2 x = =0 (14) = 2 dt dt dt dt Knowing that the value of u is constant along a characteristic, we can define a series of parametric lines that only depend on the initial value of u at a point emanating from (η, 0). Direct integration of dx = c. dt (15) x = ct + ξ (16) ξ = x − ct (17) results in the characteristic line or where the solution of u(x, t) of the initial value problem is given by the u(x, t) = u(ξ, 0) = u0 (ξ) = u0 (x − ct) (18) Problem 2 [2 Points] Solve the initial value problem for the model hyperbolic equation, ∂u ∂u +c =0 ∂t ∂x using the implicit central difference method, given as Method : un+1 = uni − c∆t i Stability Criteria : Always stable  Order of Accuracy : O ∆t, ∆x2 (19) D0 · n+1 u ∆x i (20) (21) (22) Use the same mesh as used in Problem 1, but for c∆t ∆x use the values 0.9, 1.8, 2.0. As in Problem 1, compare the computed solutions to the analytical solution after 10 time steps. Note: The value c∆t ∆x is often referred to as the CFL number, in reference to Richard Courant, Kurt Friedrichs, and Hans Lewy, who described the stability of hyperbolic PDEs in a 1928 paper. While the criteria is a direct result of von Neumann stability analysis (deduced many years later), Courant, Friedrichs, and Lewy derived the result more intuitively. They reasoned that during a single discrete time step, a wave (characteristic) moving across a mesh must not move further than the adjacent grid point during that time step. 4 Computational Fluid Dynamics : Homework #2 Problem 3 Problem 3 [3 Points] Solve the boundary value problem for the model parabolic equation, described by the partial differential equation ∂2u ∂u (23) =ν 2 ∂t ∂x using both the full explicit (α = 0) and full implicit (α = 1), central difference methods described as  Dxx ·  (1 − α)uni + αun+1 i ∆x Stability Criteria : α = 1 : Always stable Method : un+1 = uni − ν∆t i α = 0 : ∆t ≤ ∆x2 /(2k)  Order of Accuracy : O ∆t, ∆x2 (24) (25) where the second-order spatial difference equation is • Central difference, second-order accurate – Uniform mesh Dxx · φi+1 − 2φi + φi−1 φi = ∆x ∆x2 (26) φi+1 −φi φi −φi−1 xi+1 −xi − xi −xi−1 1 2 (xi+1 − xi−1 ) (27) – Non-uniform mesh Dxx · φi = ∆x Use N = 41 points to span the interval 0 ≤ x ≤ Lx with Lx = 1. Assume the mesh is uniform in the x-direction, where the location at each point is given by xi = (i − 1)∆x with ∆x = Lx /(N − 1) (28) The parabolic PDE is a boundary value problem and requires us to specify a boundary condition at each boundary in the domain. Assume ν = 1 and use a Dirichlet boundary condition, i.e., fixed value of u at each end of the domain, x = 0 and x = Lx = 1. The initial condition and boundary conditions are given as Left boundary condition x = 0 : u(0, t) = 0 Right boundary condition x = Lx : u(Lx , t) = 1 Initial condition t = 0 : u(x, 0) = 0 (29) (30) (31) For the full explicit (α = 0) method, use the stability constraint to compute the time step, ∆t, and run for 100 time steps. For the full implicit (α = 1) method, use two values of ∆t. For one simulation use ∆t = 1/8 for 10 time steps, and for a second simulation use ∆t = 109 for 1 time step. Plot the results. The steady-state solution of u(x) should be a linear function bounded by 0 at x = 0 and 1 at x = Lx . Problem 4 [5 Points] Use the implicit central difference method as described Eqs. 24 to solve heat equation, ∂T k ∂2T ∂2T = =ν 2, 2 ∂t ρcp ∂x ∂x (32) where ν = k/(ρcp ) is the thermal diffusivity. In the previous problem, we essentially solved a homogenous, non-dimensional form of the heat equation. In this problem, let us find the solution to a more practical Problem 4 continued on next page. . . 5 Computational Fluid Dynamics : Homework #2 Problem 4 (continued) problem. Use the following values, ρ = 8000 kg/m3 , cp = 500 J/kg · K, and k = 10 W/m · K. At the left boundary, assume we have been given a known constant heat flux, q0 , and at theright boundary, we will assume the surface is adiabatic, meaning the heat flux across that boundary is zero. The initial and boundary conditions as Left BC - Specified Heat Flux, x = 0 : qw (0) = q0 = −7.5 × 105 W/m2 Right BC - Adiabatic, x = Lx : qw (Lx ) = 0.0 Initial condition, t = 0 : T (x, 0) = T0 = 300K (33) (34) (35) Note that the heat flux is defined as negative. Following sign convention, this is the heat flux into the material. We can find an analytical solution to the heat condition, which is given by the equation, 1 νt x 1 T (x, t) − T0  = 2 + − + Lx Lx 3 Lx 2 |q0 | k  x Lx 2 −2     ∞ X nπx 1 2 2 νt exp −n π cos n2 L2x Lx n=1 (36) Compare your numerical solutions to the analytical solution at t = 40s. To evaluate the Fourier series in the analytical solution, you can truncate the series at n = 1000. Note: For long times the scenario described by the boundary conditions above is unphysical, since as time continues (t → ∞) more and more energy is added to the system, q0 is constant. In practice, q0 is not constant. At high temperatures, black-body radiation at this surface must be taken into account. This would imply that the boundary condition is a function of temperature. A mixed boundary condition like this is called a Robin boundary condition. Bonus: 3 points Find a solution to the heat conduction equation as described above, but instead of using a uniform mesh, use a non-uniform mesh that is hyperbolically spaced. Compare your results to those obtained from using a uniform mesh. The mesh sizing can be computed using the following procedure. Given the length of the domain as Lx , the number of points in the mesh as N , the location of the first node index as x1 , and the minimum spacing as ∆xmin = x2 −x1 , then the x-location of each node in the mesh is given by the hyperbolic formula i−1 xi = x1 + Lx e N −1 − 1 eκ − 1 (37) where the value of ∆xmin can be used to determine stretching factor κ. Using Newton’s method (a rootfinding algorithm), we can iteratively solve for the value of κ using the equation, 1 x2 − x1 = ∆xmin = Lx e N −1 − 1 eκ − 1 (38) Once κ is determined, the mesh can be defined using Eq. 37. Use the same implicit central difference equation, however, the difference formula for the second derivative should account for the non-uniform mesh size. 6 ...
Purchase answer to see full attachment
Student has agreed that all tutoring, explanations, and answers provided by the tutor will be used to help in the learning process and in accordance with Studypool's honor code & terms of service.

Final Answer



I was on a very tight deadline but thanks to Studypool I was able to deliver my assignment on time.

The tutor was pretty knowledgeable, efficient and polite. Great service!

Heard about Studypool for a while and finally tried it. Glad I did caus this was really helpful.