Computational Fluid Dynamics Chapter 11Chapter 1 Simple initial-value problemsIn this chapter, simple, ordinary differential equations for initial-value problems will be solvedto understand the basic concept of numerical methods. Two practical examples will be used todemonstrate the principle of numerical methods. In section 1.1, the solution of a moving box along aflat plate under an external … Continue reading “Fluid Dynamics | My Assignment Tutor”
Computational Fluid Dynamics Chapter 11Chapter 1 Simple initial-value problemsIn this chapter, simple, ordinary differential equations for initial-value problems will be solvedto understand the basic concept of numerical methods. Two practical examples will be used todemonstrate the principle of numerical methods. In section 1.1, the solution of a moving box along aflat plate under an external force will be discussed. In section 1.2, the vibration of a circular cylinderin a fluid with an initial displacement will be solved. The Euler method and the Fourth-order LungeKutta method will be introduced.2.1 Acceleration of a box under an external forceWe consider a very simple case of a box moving horizontally on a flat plate under an externalforce F as shown in Figure 1. Now we will use a very simple numerical method to calculate thedisplacement (position) of the box at any time t.Figure 1 A box moves along a flat plate under a force FIf the box’s mass is m and it is initially static. The acceleration of the box (a) isFma . Thevelocity V is the first-order derivative of x with respect to the time, t, asxtV d(1)Initially, both the displacement x and the velocity V are zero, and they increase with time. Theacceleration is defined as the increase of the velocity in unit time, i.e.aV(2) dtddThe solution of Eq. (2) isV (t) at V0 (3)where V0 is the velocity at t=0, i.e. the initial velocity. If initially the box is static, V0=0. In practice,what we want to know is the position of the box at any time t. If we use x to stand for thedisplacement of the box, the derivative of x with respective to time is the velocity, i.e.atxtdd(4)The problem described above is a very simple problem that we can solve theoretically. Inpractical engineering, many problems do not have theoretical solutions. Then they have to be solved xF Initial condition (at t=0), x=0, V=0Computational Fluid Dynamics Chapter 12by numerical method. Numerical methods use numerical approximation to solve a problem. Equation(4) is the governing equation for calculating x.Now we solve Eq. (4) Using the numerical method to see how the numerical method works. Tosolve Equation (4), we need to know the initial displacement of the box. The problem is an initialvalue problem. Initially, we know the displacement of the box is zero, i.e.0x t 0 (5)Equation (5) is the initial condition of the problem, i.e. Equation (4). At t=0, the velocity is zero.If we increase the time gradually by a very small interval t and use xnt to stand for thedisplacement of the box at nt . We assume the velocity at a very small time interval from t nt to(n 1)t does not change and is the same as the velocity at t nt . If the time interval t is smallenough, this assumption does not lead to too much error to the solution.In CFD, the time interval t is usually called computational time step, or time step.Taking a=0.5 m/s2, and the time interval Δt=0.1 s as an example, the procedure for calculatingthe displacements x1t , x2t , , xnt , are (Note: the distance that the box moves from t nt to(n 1)t is V ntt1. t 0 to t , 0 V , displacement x1t 0 0 t 02. t t to 2 , t V 1t 0.5 0.1 0.05 m/s , x2t x1t V 1t t 0 0.05 0.1 0.005 m3. t 2 to t 3 , t V 2t 0.5 0.2 0.1 m/s , x3t x2t V 2t t 0 0.1 0.1 0.015 m4. t nt to (n 1)t , V nt 0.5 nt , x(n1)t xnt V nt tContinuing the above steps, we can find the displacement at any time. It should be noted that theabove method is an approximation method, because it is assumed that the velocity does not changefrom t nt to (n 1)t .Figure 2 Solutions to the problem of a box moving on a flat plate with an acceleration of 0.5 m/s200.050.10.150.20.250 0.2 0.4 0.6 0.8 1x (m)t (s)TheoreticalΔt=0.1Δt=0.05Δt=0.02Computational Fluid Dynamics Chapter 13The theoretical solution to Eq. (4) is212x at (6)To see how accurate the numerical results are, Figure 2 shows the comparison between thenumerical solutions and the theoretical solution for the problem defined in Figure 1, with anacceleration of a=0.5 m/s2. It can be seen the smaller the time interval t , the closer between thenumerical results and the theoretical solution.The example discussed above has theoretical solution. However, many engineering problemscannot be solved by the theoretical method because of the complexity of the problems. In these cases,they have to be solved by numerical method.Above method for solving a differential equation is called Euler method. If we have adifferential equationf x tdtdx , (7)with the initial conditionx t 0 x0 (8)If we choose a time step Δt. Based on the initial value, we can use the Euler method to predict x at Δt,2Δt, 3Δt, … usingEuler Method:xt x0t f 0ttx2t x1t f 1ttx(n1)t xnt f ntt (9)The Euler method is a very simple and basic method. To increase the accuracy, we can use moreaccurate methods, for example the Fourth order Runge-Kutta method shown below.Fourth-order Runge-Kutta method:( 1) 1 2 2 2 3 46x n t xnt t k k k k (10)where( , )( / 2, / 2)( / 2, / 2)( , )4 33 22 11k f n t t x k tk f n t t x k tk f n t t x k tk f n t xn tn tn tn t Figure 3 shows the comparison between the Euler method and the Fourth-order Lunge-Kuttamethod for the differential equation dx/dt=0.5t. It can be seen that, using the same time step ofΔt=0.1 s, the Four-order Lunge-Kutta solution is the same as the theoretical results, while the Eulermethod solution has errors. For higher order governing equations, the Lunge-Kutta solution may notgive exact solution, but its solution will still be much more accurate than the Euler method solution.Computational Fluid Dynamics Chapter 14Figure 3 Comparison between the solutions of the Euler method and the Fourth-order Lunge-Kuttamethod for dx/dt=0.5t1.2 Vibrating of a cylinder in a still fluidA second example is the vibration of an elastically mounted circular cylinder in still fluid asshown in Figure 4. At any instant, the displacement of the cylinder is represented by Y, and thevelocity and the acceleration of the cylinder are:YtVdd (6)and2 2Ya (7)When the cylinder is moving in the fluid, it experiences forces from the spring, the damper andthe fluid. The force of the spring is:Spring force = kY(8) ddtwhere K is the spring constant. The negative sign means that the force is in the opposite direction ofY. The force of the damper is proportional to the velocity as:damping force =Ytcdd (9)where c is the damping constant, or the damping coefficient. The equation of motion iskY FYtctYm dddd22(10)where F is the external force. For a cylinder vibrating in a fluid, the force term is the force from thefluid. The force from the fluid can be calculated by the empirical formula as: YdYdC AtYC mA D p12ddForce from fluid22 d tdtd (11)00.050.10.150.20.250 0.2 0.4 0.6 0.8 1x (m)t (s)TheoreticalΔt=0.1, Euler methodΔt=0.1, Fourth‐orderRLunge‐Kutta methodComputational Fluid Dynamics Chapter 15(a) Three-dimensional view (b) Cross-sectionFigure 4 Elastically mounted circular cylinder in still fluidThe first term of Equation (11), which is proportional to the acceleration, is called inertia force,and the second term, which is proportional to the velocity square, is called the drag force. Theparameters in Equation (11) are defined as: ρ = fluid density CA=added mass coefficient, in this example CA=1md =displaced fluid mass, m D L42d ρ=density of the fluid D=diameter of the cylinder L=length of the cylinder CD=drag coefficient, in this example CD=1.2 Ap=DL, projected area of the cylinder in the flow directionIt should be noted that the force from the fluid is calculated from an empirical formula and theforce coefficients in Equation (11) are also empirical numbers that have uncertainty. There is notheoretical method to calculate the force of the fluid.Based on the Newton’s second law of motion, the total force equals to the mass multiplied bythe acceleration of the cylinder, we have YdYdtdtd C AtYC mYtkY ctYmA D p 12dddddd222 d2 (12)By rearranging the terms in Equation (12) we get the equation of motion in a fluid as: 0dddd12dd22d m CAm tY c CD Ap Yt Yt kY (13) YtVddY Balance positionV, aYtYtC Adtd YC VYtckYA f D dddd12Force from fluidddforce from damperForce from Spring22 Computational Fluid Dynamics Chapter 16Comparing Equation (13) with Equation (11), it can be seen that the fluid leads to increase in theeffective mass and the damping coefficient. CAmd is called the added mass. The effective mass isincreased because when the cylinder moves, it also drives its surrounding fluid to move.Now we solve the equation of motion, Equation (13) using the fourth-order Runge-Kutta method.The second-order equation (13) can be rewritten into two first-order equations:vYtdd(14)( , )ddf Y vvt(15)whered12( , , )m C mc C A v v kYf Y v tp AD (16)The Fourth-order Runge-Kutta formulae for Equations (14) and (15) are( 1) 1 2 2 2 3 46 Y Y Y YY n t Y nt t k k k k (17)( 1) 1 2 2 2 3 46 v v v vv n t vnt t k k k k (18)wheren tkY1 v (19-1)kv1 f Y nt ,vnt (19-2)2 12 vn tkY v t k (19-3) 2 1 12,2 vn tYn tkv f Y t k v t k (19-4)3 22 vn tkY v t k (19-5) 2 3 2 2, 2n tYn tt k vvt kkv f Y4 v3kY v tkn tkv f Y tkv tkYv (19-6)(19-7)4 nt 3,nt 3 (19-8)Now we solve an example problem of free vibration of a cylinder in fluid fluid using the fourthorder Runge-Kutta method. The parameters are listed in table 1. Initially the cylinder’s displacementis 0.5 m and its velocity is zero. The results of the example are shown in Figure 5, where thedisplacement-versus-time curves are plotted. It can be seen that the vibration amplitude in thevacuum does not change with time, because there are not any resistance from the fluid and thedamping coefficient c is zero. However, when the cylinder is placed in water, the amplitude reducedvery quickly due to the strong resistance from the water. If the fluid is air, the amplitude also reduces,Computational Fluid Dynamics Chapter 17but with very small rate. This is because the density of the air is much smaller than that of the water,resulting much smaller resistance force on the cylinder.Table 1 parameters for vibration of a cylinder in fluid DiameterD (m)0.1LengthL (m)1Density of the cylinderρcylinder (kg/m3)2000Density of the fluidρ (kg/m3)0 (vacuum), 1.204 (air),1000 (water)Spring stiffnessk (N/m)160Damping coefficientc (Ns/m)0Drag coefficientCD1.2Added mass coefficientCA1 Figure 5 Solution of the vibration of a cylinder in a fluid using the parameters listed in Table 1.ExerciseDevelop a MATLAB program to solve the problem of an elastically mounted circular cylindervibrating in a still water with an initial displacement of 0.2 m. Initial displacement, Y=0.2 m at t=01 m/sPeriod of the oscillatory flow, T10 sDiameter of the cylinder, D0.1 mLength of the cylinder, L1 mMass of the cylinder, m6 kgDensity of the water, ρ1024 kg/m3Stiffness of the spring (total), K200 N/mDamping coefficient (total), c10 N·m/s Solution‐0.6‐0.4‐0.200.20.40.60 5 10 15 20 25 30x (m)t (s)Vacuum: c=0, ρ=0 Air: c=0, ρ=1.204 kg/m3 Water: c=0, ρ=1000 kg/m3Computational Fluid Dynamics Chapter 18The equation of motion and the Fourth-order Runge-Kutta method has been discussed in Chapter 1.Following the MATLAB program.pi = acos(-1) % pi = 3.14159D = 0.1 % cylinder diamterL = 1 % cylinder lengthm = 6 % cylinder massro = 1024 % fluid densityK = 20 % spring stiffnessc = 0 % damping constantCA = 1 % added mass coefficientCD = 1.2 % drag coefficientmd = ro*pi*(D*D)/4*L % added massAp = D*L % projecting areadt = 1/60 % time stepndt = 1000 %total step to be calculatedX(1:ndt+1) = 0 % displacement X from step 0 to step ndtv(1:ndt+1) = 0 % velocity from step 0 to step ndttime(1:ndt+1) = (0:ndt)*dt % time from step 0 to step ndtX(1) = 0.2 % initially displacement is 0.2 mfor n=1:ndtt = time(n)% to calculate kx1 and kv1vtemp = v(n)xtemp = X(n)t2 = 0.5*ro*CD*Ap*abs(vtemp)*(-vtemp) % drag forcet3 = -K*xtemp-c*vtemp % spring and dammping forceskx1 = v(n)kv1 = (t2+t3)/(m+CA*md)% to calculate kx2 and kv2t = time(n)+dt/2vtemp = v(n)+0.5*dt*kv1xtemp = X(n)+0.5*dt*kx1t2 = 0.5*ro*CD*Ap*abs(vtemp)*(-vtemp)t3 = -K*xtemp-c*vtempkx2 = v(n)+0.5*dt*kv1kv2 = (t2+t3)/(m+CA*md)% to calculate kx3 and kv3t = time(n)+dt/2vtemp = v(n)+0.5*dt*kv2xtemp = X(n)+0.5*dt*kx2t2 = 0.5*ro*CD*Ap*abs(vtemp)*(-vtemp)t3 = -K*xtemp-c*vtempkx3 = v(n)+0.5*dt*kv2kv3 = (t2+t3)/(m+CA*md)% to calculate kx4 and kv4t = time(n)+dtvtemp = v(n)+0.5*dt*kv3xtemp = X(n)+0.5*dt*kx3t2 = 0.5*ro*CD*Ap*abs(vtemp)*(-vtemp)Computational Fluid Dynamics Chapter 19t3 = -K*xtemp-c*vtempkx4 = v(n)+dt*kv3kv4 = (t2+t3)/(m+CA*md)% next step valueX(n+1) = X(n)+(dt/6)*(kx1+2*kx2+2*kx3+kx4)v(n+1) = v(n)+(dt/6)*(kv1+2*kv2+2*kv3+kv4)end% write the results into a file “output.txt”fileID = fopen(‘output.txt’,’w’);for n=1:ndt+1fprintf(fileID,’%12.6f %12.6f %12.6frn’,time(n),X(n),v(n));endfclose(fileID);% plot the pictureplot(time(1:ndt),X(1:ndt))The above program plots a displacement (vertical axis) versus time (horizontal axis) shownbelow. It also writes the results into a file “output.txt”. This file can be used to generate diagramsusing other software such as Miscroft Excel.Figure 6 The plot generated using the program in the example. The horizontal axis is time and thevertical axis is displacement of the cylinder.0 2 4 6 8 10 12 14 16 18-0.1-0.0500.050.10.150.20.25
