Download Exercises Computational Fluid Dynamics A.E.P. Veldman ... and more Schemes and Mind Maps Fluid Dynamics in PDF only on Docsity!
Exercises Computational Fluid Dynamics
A.E.P. Veldman (University of Groningen)
General remarks
During the course several computer exercises have to be made. These exercises are to be solved in ‘pairs’ of one or two students. A (short) written report of the findings has to be handed in the following week; discussion of the homework with other ‘pairs’ of CFD students is not allowed during that period.
All software required for the exercises is available by downloading the files from the Nestor pages for this course:
Computational Fluid Dynamics > Course Information
An alternative is to download them through
http://www.math.rug.nl/∼veldman/Colleges/CFD/CFD-Exercises.html
Exercise 1 – Artificial diffusion
Files and commands Required files: cfd 1.m, cfd upw.m, ns exact.m The program is started in Matlab with the command cfd 1.
Description
Consider the inhomogeneous convection-diffusion equation
dφ dx − k
d^2 φ dx^2 = S on 0 ≤ x ≤ 1 , with φ(0) = 0, φ(1) = 1.
The right-hand side is given by
S(x) =
{ 2 − 5 x , 0 ≤ x ≤ 0. 4 , 0 , 0. 4 < x ≤ 1.
In this exercise an upwind discretization of the above equation is applied. The discrete solution for three values of the diffusion coefficient k is plotted in one picture. Also the exact solution for one value of k can be included in this picture.
Questions
- Firstly, study the discrete solution for a grid with N = 20 grid points, at the following values of k: k = 0.1, 0.01 and 0.001. These values correspond with mesh-P´eclet numbers
P = 0.5, 5 and 50, respectively. One would expect (or at least hope) that these solutions tend to the (analytical) solution for k = 0. Hence, compare the discrete solutions to the exact solution at k = 0 (ignore the Matlab warnings concerning division by zero). What do you observe; is there any such relation?
- Next, compare the same three discrete solutions at k = 0.1, 0.01 and 0.001 (again for N = 20) to the exact solution at k = 0.025. Now, there seems to be some relation between the upwind limit for k → 0 and the exact solution at k = 0.025. Explain this relation.
- Finally, study once more the upwind solutions at k = 0.1, 0.01 and 0.001, but now on a grid with N = 40 grid points. Which exact solution, i.e. at which value of k, corresponds with the upwind limit for k → 0?
- Finally, consider the results on the stretched grid. Compare the upwind method and both generalizations (A and B) of the central method (cf. Section 1.6). Describe the difference between the discrete solutions; why has the upwind solution not improved?
- On the stretched grid, determine empirically the value of k for which one of the eigen- values of the coefficient matrix of Method B crosses the imaginary axis. Furthermore, try to find a value of k for which the coefficient matrix is approximately singular. Watch how the solution of Method B behaves at this singularity.
Questions to be solved by pencil-and-paper
- Prove that λ ≥ max( 12 − (^) uhk , 0) is a sufficient condition for the upwind-biased lambda schemes from Eq. (1.15) on page 12 to be wiggle-free. In particular show that the QUICK method is wiggle-free for P ≤ 8 /3. Hint: Try fundamental solutions of the form ri, and monitor the sign of r.
- Discretize the convection-diffusion equation
dφ dx − k
d^2 φ dx^2 = 0 on 0 ≤ x ≤ 1 ,
with Dirichlet boundary conditions
φ(0) = 0, φ(1) = 1.
Use a nonuniform grid with only one interior grid point at x = 1 − 2 k; k remains an unspecified variable here. Use the discretization methods A and B. In both cases, A and B, solve the discrete system (which in fact is only one equation). Sketch both solutions using linear interpolation between the grid points. Compute the exact solution of the convection-diffusion equation (at x = 1 − 2 k) as well. Which discrete solution is better?
Exercise 3 – The JOR and SOR method
Description
After discretization of a (system of) partial differential equation(s) one obtains a large system of algebraic equations, Ax = b. Solving such a large system with a direct method is expensive with respect to both memory and CPU time. Therefore, such systems are often solved with an iterative method.
As in Exercise 1 and 2, we consider the inhomogeneous convection-diffusion equation
dφ dx
− k d^2 φ dx^2
= S on 0 ≤ x ≤ 1 , with φ(0) = 0, φ(1) = 1.
The right-hand side is given by
S(x) =
{ 2 − 5 x , 0 ≤ x ≤ 0. 4 , 0 , 0. 4 < x ≤ 1.
The Matlab file cfd 3.m discretizes this equation with a number of finite-difference (-volume) methods, and solves the resulting linear systems iteratively.
Two types of grid are used: uniform grids and non-uniform (stretched) grids. The stretched grid possesses N/2 equal grid cells between 0 and 1 − 5 k, and N/2 equal grid cells between 1 − 5 k and 1. Several types of discretization are employed; on the uniform grid both an upwind and a central discretization, and on the stretched grid an upwind, a central A and a central B discretization are used.
The resulting linear systems are solved iteratively with the JOR and SOR methods. As stopping criterion in these methods
‖ φn+1^ − φn^ ‖∞≤ 10 −^5
is used, in which n is the counter of the iterative steps; the maximum norm of the difference between the solutions of two consecutive iteration steps must be sufficiently small.
Required files
This exercise requires the files cfd 3.m, plotev.m, plotsol.m, ns iterative.m, ns jor.m, ns sor.m, ns exact.m.
Typing the command cfd 3 in Matlab starts a menu in which the choice of the grid, the discretization method, the iterative method and the value of the relaxation factor ω can be given. Then, with the menu a program can be started which solves the convection- diffusion equation in the prescribed way. The results produced by the program consist of the eigenvalues of the Jacobi matrix, a plot of the iteration error (this error is written to the screen as well), the total number of iterations and a plot showing the continuous solution, the discrete solution computed with a direct solver and the discrete solution computed with the iterative method (provided the iterative method did converge).
stretched Jacobi method Gauss-Seidel method grid convergent? # iterations convergent? # iterations upwind central A central B
each value of ω the number of iterations needed by the JOR method. Determine in this way the optimal and maximal relaxation factor, ωopt and ωmax for the JOR method. Do the same for the SOR method. Put the results in the table below.
stretched JOR method SOR method grid ωopt # it. ωmax # it. ωopt # it. ωmax # it. upwind central A central B
In case of the central B discretization compare the discrete solution obtained with the JOR or SOR method with the discrete solution obtained with the exact solver. Do the JOR and SOR method really converge? Use the eigenvalues of the Jacobi matrix, computed by cfd 3, to explain this result. Which systems are solved better by the JOR and SOR method, those resulting from the upwind discretization, those resulting from the central A discretization or those resulting from the central B discretization?
Questions to be solved by pencil-and-paper
- In Section A.3.1 for a central discretization on the uniform grid the Jacobi matrix and its eigenvalues are derived. The Jacobi matrix reads
diag[
P
P
],
(P is the mesh-P´eclet number) and its eigenvalues μJ are given by
μJ =
√ 4 − P 2 cos( iπ N ), i = 1,... , N − 1 ,
where N is the number of grid cells. Derive in a similar way for the upwind discretization on the uniform grid the eigenvalues μJ of the Jacobi matrix. Show that these eigenvalues are real and between -1 and 1 for every value of P. Hint: the eigenvalues of the tri-diagonal matrix diag(a, 0 , c) are given by μ = 2
ac cos(iπ/N ), i = 1,... , N − 1.
- For a central discretization on the uniform grid the eigenvalues of the Jacobi matrix are pure imaginary when P > 2. Let the largest pure imaginary eigenvalue be given
by μJ = iμI , with μI real. The optimal and maximal relaxation factor for the JOR method are given by (cf. Section A.3.1)
ωJOR,opt =
1 + μ^2 I , ωJOR,max =
1 + μ^2 I
For an upwind discretization on the uniform grid the eigenvalues of the Jacobi matrix are real (see Question 5). Let the largest real eigenvalues be given by μJ = μR and μJ = −μR, with μR real. Derive, as a function of μR, the optimal and maximal relax- ation factor for the JOR method. Hint: the eigenvalues μJOR of the iteration matrix of the JOR method are related to the eigenvalues μJ of the Jacobi matrix by μJOR = 1 + ω(μJ − 1), and the JOR method converges when the eigenvalues μJOR lie within the unit circle. Furthermore, for the SOR method give the optimal and maximal relaxation factor (as a function of μI and μR respectively) when using the central and the upwind discretiza- tion.
- Use Questions 5 and 6 (remember that k = 0.01 and N = 10) to determine for the JOR method the values of the optimal and maximal relaxation factor when solving the systems resulting from the discretization on the uniform grid with respectively the upwind and central method. Do the same for the SOR method. Compare these theoretical results with the numerical results from Question 2.
- Consider the upwind discretization on the stretched grid. Use Question 6 (remember that k = 0.01 and N = 10) and the numerical values of the eigenvalues μJ of the Jacobi matrix, computed as a byproduct by cfd 3 in Question 4, to determine for the JOR method the theoretical values of the optimal and maximal relaxation factor. Do the same for the SOR method. Compare these theoretical results with the numerical results for the upwind discretization on the stretched grid from Question 4. Why is it difficult to compute the theoretical optimal and maximal relaxation factor when solving the systems resulting from a discretization on the stretched grid with a central A or central B method?
time step δt = 0. 05 δt = 0. 5 stable? stable? damping ω = 0. 5 ω = 0. 6
- Finally, again for N = 10, investigate the fully implicit method (ω = 1). Use the same time steps as in Question 2. Observe whether or not the solution shows oscillations.
time step δt = 0. 05 δt = 0. 5 stable? oscillations? stable? oscillations? ω = 1. 0
Questions to be solved by pencil-and-paper
First consider the fully explicit method (ω = 0). Carry out a Fourier analysis of this method. Determine the maximum allowable time step for N = 10 and N = 20. Compare these stability limits with the empirical observations in Question 1.
Carry out a Fourier stability analysis of the generalized Crank-Nicolson method, and determine the amplification factor. Show that the generalized Crank-Nicolson method is unconditionally stable for ω ≥ 0 .5. Investigate how the amplification factor behaves when δt → ∞. Now explain the oscillations visible in Question 2. How are these oscillations influenced when ω is increased from 0.5 to 0.6? Give a possible explanation why the observed amplification factors in Question 2 are somewhat different from the theoretical amplification factors computed in this Question.
Derive theoretically for which values of the Crank-Nicolson parameter ω the solution is wiggle-free. As a special case, explain why the fully implicit method (ω = 1) does not show oscillations (see Question 3). Hint: Use the concept of a positive operator.
Exercise 5 – Navier-Stokes solver
General description
After semi-discretization with an explicit time-discretization method with two time levels, the Navier-Stokes equations can be written as
div un+1^ = 0, un+1^ − un δt
in which R = −(u.grad)u + ν div grad u consists of the convection and diffusion terms. Rewriting and combining these two equations gives
div grad pn+1^ = div ( un δt
- Rn), (1a) un+1^ = un^ + δt Rn^ − δt grad pn+1. (1b)
The fortran program cfd 5.f solves the Navier-Stokes equations.
Files and commands
This exercise requires the files cfd 5.f, cfd.in, liqdmenu.m, liqdplot.m, liqprint.m, xplot.m.
The program cfd 5.f has to be compiled with the command gfortran -o cfd 5 -O cfd 5.f (any other Fortran compiler, e.g. ifort, will also do). Then, the program can be run with the command ./cfd 5. This program needs the input file cfd.in, this file has to be adjusted to perform the different computations. The input file, with the parameter setting of Question 0, is included at the end of this exercise. At the end of the input file an explanation of the different parameters is given.
The program writes output both to the screen and to files. The output to the screen first gives the values of the parameters describing the problem, and then gives each (user specified) time interval the time, div u, the number of Poisson iterations and the horizontal velocities at three positions at the vertical center line of the domain (at the bottom, halfway, and at the top). Furthermore, output is written to the files uvpf#.dat and config.dat. The output file config.dat contains a description of the grid. The output file uvpf#.dat contains the horizontal and vertical velocity and pressure in each grid cell at a certain time. In the input file the user can specify whether the file uvpf#.dat is made each (user specified) time interval or only at the end of the computation.
The output files uvpf#.dat and config.dat are needed for the post-processing in Matlab. Typing the command liqdmenu in Matlab starts a menu with which the post-processing can be controlled. With the upper part of the menu surface or contour plots of the (different) velocities, pressure, vorticity or streamfunction can be made and with the lower part of the menu a movie of these plots at consecutive time steps can be made. With the middle part of the menu plots of the (different) velocities or the pressure along horizontal or vertical sections of the domain can be made.
written to the screen and the velocity field is written to the file uvpf0001.dat at the end of the simulation.
Question 0: ’exact’ solution
First the ’exact’ solution is computed, i.e. the discrete steady solution on the given grid. Hereto set in the input file Divuv=1 (i.e. retain the term div un^ in the Poisson equation), StrtP=1 (i.e. use the previously found pressure as initial guess for the iterations in the new time step) and ε = 10−^3. Note that the solution has converged in time, and write the results at t = 25 in the table below.
time div u # it. u bottom u mid u top 25
Questions without div un
First, substitute div un^ = 0 in the Poisson equation (i.e. set Divuv=0 in the input file). Hence, solve the Navier-Stokes equations in the form (2a) and (2b). In each of the following questions, put the output of the program at t = 20, 21 ,... , 25 in a table of the format given below.
time div u # it. u bottom u mid u top 20 21 22 23 24 25
- Use for the Poisson solver an accuracy ε = 10−^1 and use as initial solution for the Poisson solver p = 0, i.e. start the Poisson solver from scratch.
a) Look what happens with the term div u in time. Furthermore, make a plot of the solution; e.g. use the post-processing menu in Matlab to view the flow, use the upper part of the menu to make a surface plot of the streamfunction. Also have a look at the absolute velocity (crange [0, 0.1]) and the vertical velocity. Repeat the simulation with TFin=100. Does the velocity field converge? b) Now we will consider what happens with the divergence when δt → 0: will this improve the solution? Set (temporarily) DelT=0.0025 and CFLMax=0.05 in the input file. Repeat the simulation (with TFin=25). Look what happens with the term div u in time. Compare this with the behaviour of div u in Question 1a. Explain the results.
Use for the Poisson solver an accuracy ε = 10−^3 , and use as initial solution for the Poisson solver p = 0, i.e. start the Poisson solver from scratch. Look what happens with the term div u in time. Furthermore, make various plots of the solution as in the previous question. The solution has improved, but is it really accurate?
Use for the Poisson solver an accuracy ε = 10−^3 , and use as initial solution for the Poisson solver the solution from the previous time step. Look once more what happens with the term div u in time. Are you satisfied now?
Use for the Poisson solver an accuracy ε = 10−^1 , and use as initial solution for the Poisson solver the solution from the previous time step. Once again, describe the behaviour of div u in time.
Questions with div un
Next, retain the term div un^ in the Poisson equation (i.e. set Divuv=1 in the input file). Hence, solve the Navier-Stokes equations in the form (1a) and (1b).
- Use for the Poisson solver an accuracy ε = 10−^1 and use as initial solution for the Poisson solver p = 0, i.e. start the Poisson solver from scratch.
a) Look what happens with the term div u in time. How do you like the solution? b) Now we will consider what happens with the divergence when we let δt → 0. Set (temporarily) DelT=0.0025 and CFLMax=0.05 in the input file. Repeat the simulation. Look what happens with the term div u in time. Compare the results with those from Question 5a, and explain the differences.
Use for the Poisson solver an accuracy ε = 10−^3 , and use as initial solution for the Poisson solver p = 0, i.e. start the Poisson solver from scratch (set StrtP=0 in the input file). Look what happens with the term div u in time. Compare the results with those from Question 5a, and explain the differences.
Use for the Poisson solver an accuracy ε = 10−^3 , and use as initial solution for the Poisson solver the solution from the previous time step (i.e. set StrtP=1 in the input file). Monitor div u in time, and make plots of the solution. Explain the results.
Use for the Poisson solver an accuracy ε = 10−^1 , and use as initial solution for the Poisson solver the solution from the previous time step. Monitor div u and plot the solution. How do you like the results?
Summary We have seen that the various approaches give quite different results. In the last question of this part you will have to decide which approaches give acceptable results.
- In practical applications one is often satisfied when the final solution has an accuracy of two digits. Keeping this in mind, compare the results from the Questions 1 to 8 and give your opinion about which approaches give acceptable results?
In every question below, plot the flow at the first five seconds of the simulation and at the end of the simulation, i.e. at t = 1, 2 ,... , 5 and t = 15. Use the post-processing menu in Matlab to plot the flow. Use the upper part of the menu to make a surface plot of the abso- lute velocity. Furthermore, use the middle part of the menu to make a plot of the horizontal velocity along the horizontal center line (i.e. j-plane=10). By using the ’hold’ button in the menu the plots of the horizontal velocity along the horizontal center line at t = 1, 2 ,... , 5 and t = 15 can be put in one figure. A boundary layer and irregularities of the flow can easily be seen in the second plot.
Questions with prescribed outflow
Firstly, prescribe a horizontal velocity at the right (outflow) boundary of the channel as well: u = 1 and v = 0 (i.e. set right=8 in the input file).
Simulate the flow for ν = 0.01. Make plots of the flow as described above. Describe and explain the development in time of the horizontal velocity along the horizontal center line.
Simulate the flow for ν = 0.005. Make plots of the flow as described above. Compare the flow with the one computed in Question 10, and explain the differences.
Simulate the flow for ν = 0.001. Make plots of the flow as described above. Com- pare the flow with those computed in Question 10 and 11, and explain the differences. Furthermore, make a surface plot of the absolute velocity in which the vector plot is included (select in the menu ’arrows automatic’ instead of ’no velocity vectors’), and zoom first in on the x-interval [0,1] and then on the x-interval [9,10] (i.e. take ’manual x-axis’ instead of ’automatic zoom’ and select the interval). Furthermore, make a plot of the horizontal velocity at the vertical section at the i-plane 32. What happens with the velocity at the upper and lower boundary?
Questions with free outflow
Secondly, allow at the right boundary of the channel a free outflow condition (i.e. set right= in the input file).
Simulate the flow for ν = 0.01. Make plots of the flow as described above. Compare the flow with the one computed in Question 10, and explain the differences. What is the limit value of the horizontal velocity at the center line y = 0.5. When the flow is fully developed the (horizontal) velocity profile is a parabola. Compute this parabola theoretically. What is the theoretical value of the horizontal velocity of the fully developed flow at the center line y = 0.5. Compare this result with the numerical value. Also make a cross sectional plot of the horizontal velocity near the end of the channel to check the parabolic profile.
Simulate the flow for ν = 0.005. Make plots of the flow as described above. Compare the flow with those computed in Question 11 and 13, and explain the differences. What
is the limit value of the horizontal velocity at the center line y = 0.5. Compare this numerical value with the theoretical value computed in Question 13.
- Simulate the flow for ν = 0.001. Make plots of the flow as described above. Compare the flow with those computed in Question 12, 13 and 14 and explain the differences. What is the value of the horizontal velocity at the center line y = 0.5 at the end of the channel. Explain why this numerical value of the horizontal velocity at the center line y = 0.5 differs from the theoretical value computed in Question 13. Also make a cross sectional plot of the horizontal velocity near the end of the channel to see the shape of the velocity profile: is it a parabola? Repeat the simulation with a channel of respectively length 20 and length 50 (set first Xmax=20 and then Xmax=50), use TFin=50 in order to reach a steady-state solution; also increase Imax (although it has little influence). Set uvp=0 in order to save memory. Make only at the end of the simulation (i.e. at t = 50) a plot of the horizontal velocity along the horizontal center line (i.e. j-plane=10). What is the limit value of the hor- izontal velocity at the center line y = 0.5. Again make a cross sectional plot of the horizontal velocity near the end of the channel to see the shape of the velocity profile: is it a parabola now? Find a physical explanation of the results.
BOUNDARY CONDITIONS AND INFLOW CHARACTERISTICS left,right,top,bottom : type of boundary condition 1=slip 2=no-slip 7=free outflow : as boundary condition the pressure is set at 0 8=inflow : the velocities are prescribed UIn,VIn : in case of inflow conditions, prescribed velocities the sign of UIn/VIn corresponds to the x/y-direction (hence it does not necessarily correspond to the inward pointing direction)
DISCRETIZATION Alpha : upwind parameter, discretization of convection term 0 = central 1 = upwind Divuv : treatment of divergence 0 = substitute div u^n=0 in Poisson equation 1 = maintain div u^n in Poisson equation
POISSON ITERATION PARAMETERS Epsi : convergence criterion of Poisson solver ItMax : maximal allowed number of iterations of Poisson solver Strtp : initial field of Poisson solver 0 = zero field p= 1 = solution from previous time step, p=poud;
TIME STEP TFin : endtime of simulation DelT : timestep (is adjusted during the simulation in order to keep the CFL number smaller than CFLMax) CFLMax : maximum allowed CFL number
** PRINT/PLOT CONTROL** PrtDt : time between 2 small printouts (to the screen and a file) uvp : velocity and pressure in Matlab format 0=only at end of computation (preferable) 1=at every small printout (takes a lot of memory)
