Studying the Viscous Flow Around a Cylinder Using OpenFoam

Studying the Viscous Flow Around a Cylinder Using OpenFoam –2– ABSTRACT We model the incompressible, viscous flow of a fluid around a cylinder using the OpenF oam toolbox. We develop a mesh with varying resolutions to properly capture the flow around the cylinder. We run our simulations for varying Reynolds number (R) cases between 32≤R≤161, and see an increase in shedding within the flow for increasing Reynolds number. We find decent qualitative agreement with observations, though quantitatively the simulated Strouhal number is not in good agreement with the expected Strouhal number. Given the conditions of our simulations, OpenF oam is able to accurately model the change in Strouhal number as Reynolds number changes. Additionally, OpenF oam is confirmed to follow the law of similarity. 1. Introduction1 Viscous flow around a cylinder has been an oft studied topic in the fluid community. One of the first well-known cases of studying flow around a cylinder was carried out by Hiemenz (1911), who took quantitative measurements of the pressure around a cylinder. Later, Homann (1936) used oil to visualize the fluid flow around cylinders. Tests of flow around a cylinder, such as via simulations, have piqued the imagination of many scientists and engineers, and thus have continued to this day. An incompressible, viscous fluid is used to model the flow around a cylinder. The equations for an incompressible viscous fluid are governed by 1 Introduction adopted from Secs. 5.4 and 5.5 of “The Physics of Fluids and Plasma: An Introduction for Astrophysicists” by Arnab Rai Choudhuri (1998). –3– ω = ∇×v, (1) ∇·v = 0, (2) ∂ω = ∇×(v×ω) + ν∇2 ω, ∂t (3) where eq. 1 is the vorticity, eq. 2 is the condition of incompressibility, and eq. 3 is Navier-Stokes with vorticity. Here, v is the velocity of the fluid and ν is the kinematic viscosity of the fluid. For flow around a geometric object, with unit length L and unit fluid velocity V , we can rewrite the variables of eq. 3 in scaled units as x = x0 L, v = v0 V, t = t0 L V , ω = ω0 . V L (4) With these scaled variables, eq. 3 can be re-written as ∂ω 0 1 02 0 0 0 0 = ∇ ×(v ×ω ) + ∇ ω, ∂t0 R (5) where ∇0 = L∇ is the gradient in space for the scaled length and R is the Reynolds number, given by R= LV . ν (6) If the Reynolds number were the same for two cases with different flows around similarly shaped objects, then the flow patterns should be the same (assuming the boundary –4– Fig. 1.—: Viscous flow around a cylinder for high Reynolds number. Originally from Homann (1936), but reproduced from Batchelor (1967). conditions are identical). Based on this observation, Reynolds (1883) developed the law of similarity. Stokes (1851) studied the effect of low Reynolds number. For the case that R << 1, eq. 4 reduces to 1 02 0 ∇ ω = 0. R (7) A result of this analysis was that Stokes found the drag force on a sphere moving through a viscous fluid. The drag force for this low Reynolds number case is FD = 6πµaU, (8) where a is the radius of the sphere, µ is the coefficient of viscosity and U is the velocity of the sphere through the fluid. This force acts in the direction opposite to the motion of the sphere, and it is known as Stokes’s Law. However, this law is not applicable in the high Reynolds number situation. –5– For high Reynolds number, where R >> 1, the drag force changes. In this case, Stokes’s Law is given by 1 FD = CD πa2 ρU 2 , 2 (9) where CD is a dimensionless drag coefficient defined as CD = 12µ 12 . aρU R (10) Thus, in this case, the coefficient will decrease with increasing Reynolds number. Figure 1 shows the flow past a cylinder for different Reynolds numbers. This figure, originally taken from Homann (1936), but reproduced from Batchelor (1967), shows the appearance of vortices behind the cylinder with high Reynolds number. Once the Reynolds number reaches around R = 60, the vortices shed into a string of vortices behind the cylinder. The properties of this string was studied by von Karman (1911), and was named the Karman vortex sheet. Not shown in Fig. 1 is that for Reynolds number on the order of 104 , a turbulent wake will develop behind the cylinder. Another element of the viscous flow around an object, like a sphere or cylinder, is that v = 0 at the solid boundary. In order to allow for this boundary condition, a boundary layer must exist where the velocity of the flow can quickly go to zero at the surface of the object. Within this layer, ∇2 v must be non-trivial, thus the viscosity term in the Navier-Stokes equation must be important. Regardless of the Reynolds number, viscosity plays a significant role in the layers next to a solid within a fluid. The boundary layer thickness can be approximated as r δ= νL . V (11) –6– From this equation, we can see that the as the fluid moves further down the length of the object, the boundary layer thickness will increase as the square root of the distance. Since viscosity is not negligible in the boundary layer, vortices in the flow can be generated downstream if none existed upstream. For large Reynolds numbers, strong velocity shears can develop, which produce turbulence. Thus, large Reynolds numbers can generate an unstable flow. One way to model the flow around a cylinder is via computational fluid dynamics (CFD) modelers such as OpenF oam (Weller et al. 1998). OpenF oam is an open-source CFD toolbox, developed in 2004, that is able to solve complex fluid flows in systems such as through a cavity or around a cylinder. One particular solver in OpenF oam is called icoF oam. The icoF oam solver uses the incompressible Navier-Stokes equations and solves the equations via the Pressure Implicit with Splitting of Operator (PISO) algorithm. It is the icoF oam solver with the PISO algorithm which is ideal for simulating incompressible flow around a cylinder. While eqs. 1-3 are ideal for modeling the incompressible, viscous flow around cylinder, instead of eq. 3 OpenF oam solves the following form of the momentum equation 1 µ ∂v + v·∇v = − ∇p + ∇2 v. ∂t ρ ρ (12) The main difference between eq. 3 and eq. 12 is that OpenF oam uses the standard momentum equation instead of using the momentum equation for vorticity. Thus, OpenF oam does not use eq. 1 either when it solves the flow around a cylinder. The remaining sections of this paper are as follows: in Sec. 2 we describe the mesh used in our simulation. In Sec. 3 we describe the properties of the model, including the details of the different directories and how to properly setup the physics. In Sec. 4 we present the results of our simulations, and compare these results to observations. In Sec. 5 –7– Fig. 2.—: Basic cartoon of the flow around a cylinder. we summarize our conclusions. 2. Mesh2 To set up the problem of incompressible viscous flow around a cylinder, we use Fig. 2 as a reference for the basic setup. Like Fig. 2, we want a cylinder with r = 0.05 m, with the vertical of the cylinder in the z-direction. Additionally, we wish to have inflow moving in the +x direction from the left boundary (or inlet). Based on the physics and the scale of the problem, such as the size of the boundary layer around the cylinder, certain parts of the domain require higher resolution than others. While in theory the best way to solve the problem would be to create a mesh with equally sized grid points set to the smallest required scale, for the sake of computational simplicity and quickness we use a varied mesh. This mesh is based off of Fig. 3, where we increase our resolution around the cylinder to ensure the physics around the boundary layer is properly captured. In OpenFoam, we must model the cylinder in three dimensions despite only wanting 2 Adopted from Sec. 3.4.1 of “OpenFOAM Guide for Beginners” by Jordi Casacuberta Puig –8– Fig. 3.—: Structure of the mesh around half of the cylinder. a two-dimensional flow. To work around this, we set the z-direction to be one grid cell thick, so we are essentially modeling in two dimensions. We use blockmeshDict in /constant/polyM esh for our mesh creation. In our mesh, we define six boundaries: (1) top, (2) bottom, (3) inlet, (4) outlet, (5) cylinder, and (6) the front and back of the domain. To develop the mesh around the cylinder, we use the arc instruction for a circular geometry. Additionally, we define the top and bottom boundaries via symmetryplane so no rigid walls exist to alter the solution. We define the cylinder via wall, the front and back as empty, and the inlet and outlet as patch. The wall instruction sets up a physical boundary, hence why it is used for the cylinder. The empty instruction is utilized for the front and back planes of a two-dimensional geometry, which is what we have for the front and back boundaries. Lastly, the patch instruction denotes a patch where there is no geometric information about the mesh, and therefore we use this for the inlet and outlet so the fluid can flow through the boundary. The version of blockmeshDict used can be found at the end of this report in the appendix. Using this mesh description, we get the meshes seen in Fig. 4. On the left of Fig. 4 we have the top down view of the mesh. The fluid moves from the -x direction to the +x direction. We note the high degree of resolution around the cylinder. Additionally, –9– Fig. 4.—: Left: Top-down view of the mesh. Right: Inside the mesh, looking in the +x direction. the resolution is lower along the edges of domain (except in the region near the cylinder), as compared to the mesh directly behind the cylinder. This enables better resolution of the flow as it moves around the cylinder, but it does restrict the ability to see smaller fluid elements near the edges. Regardless, as we are pursuing a simplified model that is not computationally expensive, we retain this grid setup for the simulation. Additionally, on the right of Fig. 4 we have the interior to the three-dimensional box that is one grid cell in height. This is evident by the grid on the cylinder, which is on the left in this figure. This perspective is from within the box, looking in the +x direction from the -y direction. It is possible to discern the changing grid structure from this particular figure. 3. Model3 The basics of the model are located in the following three directories: 1) 0, 2) constant, and 3) system. We will discuss the model files located in these directories. The exact versions of the files used can be found at the end of this document in the appendix. 3 Adopted from Secs. 3.4.2-3.4.5 of “OpenFOAM Guide for Beginners” by Jordi Casacu- berta Puig – 10 – 3.1. 0 directory In the 0 directory, there exist two files of note: p and U . These two files define the initial conditions for the fluid flow problem. In p we give the initial pressure conditions for the system. Our internalF ield is set to be uniformly zero. Our boundaries have the conditions that the top and bottom boundaries are defined via symmetryP lane, the front and back is defined as an empty boundary, the cylinder boundary has a zero setting, and the inlet and outlet boundaries are defined via a f reestreamP ressure. The f reestreamP ressure condition sets a fixed pressure value for ingoing flow, but a zero gradient at the outgoing boundary, so this is useful for external flow. The U file gives the initial velocity conditions for the system. The internalF ield is set to have a uniform velocity in the +x direction. The inlet and outlet boundaries also maintain a f reestream condition of a uniform velocity in the +x direction, with the same velocity as that in the internalF ield. The different values for the velocity are given in Table 1. We used 5 m/s for all cases except for the first and last ones. For the first one, OpenF oam was unable to handle a high viscosity in the simulation, so we used a low viscosity and changed the fluid velocity. In the second case of R = 161 we wished to test the law of similarity so we used a different velocity and viscosity from case six. The top and bottom bottom boundaries are given the symmtetryP lane condition. The front and back boundaries are given the empty condition as per usual. Lastly, the cylinder is given a f ixedV alue condition of have zero velocity at the cylinder boundary, hence creating a boundary layer. – 11 – Case R U [km/s] ν [m2 /s] 1 32 0.82 0.002564 2 55 5.00 0.009091 3 65 5.00 0.007692 4 73 5.00 0.006849 5 102 5.00 0.004902 6 161 5.00 0.003106 7 161 4.13 0.002564 Table 1:: Reynolds number cases run in OpenFoam for flow around a cylinder with r = 0.05 m. 3.2. constant directory The constant directory contains the files which define the kinematic viscosity and the Reynolds-averaged Simulation (RAS) model. The first flag in the RAS file, called RASP roperties, is the RAS Model (denoted by RASM odel), which calls the RAS turbulence model. In this case we are using laminar flow, so we set that flag to laminar. Additionally, for the turbulence and printCoef f s flags, we set both to of f . In the viscosity file, called transportP roperties, we set the value of nu (aka the kinematic viscosity) to our desired value. Table 1 shows the different values of ν we used for each case in order to vary the Reynolds number. 3.3. system directory The system directory contains the files which tell the code how to run. We find the controlDict file which tells the code how long to run, with what time step to run, and how – 12 – to write the data. We also find f vSchemes, which tells the code the schemes to use, and f vSolutions, which tells the code what solvers to use. The important aspects of controlDict are as follows. The solver being used for flow around a cylinder is icoF oam. We set the start time to be at 0, since this is where our initial conditions are stored. For the purposes of this project, we use an end time of 2.60 for all of our cases, with a deltaT of 0.00001 and a writeInterval of 1000. For f vSchemes, it is mostly the same as for the cavity simulation, but instead under laplacianSchemes we have the def ault set to none, and we also have laplacian(nu, U ) as Gausslinearorthogonal, and likewise for laplacian((1|A(U )), p). For f vSolutions, there are two main differences from the cavity simulation. First, the U solver is P BiCG, with a preconditioner of DILU , a tolerance of 1e − 05 and a relT ol of 0. Additionally, under P ISO, we have nN onOrthogonalCorrectors set as 3 instead of 0 (as it is in the cavity case), because this helps in getting more physically accurate results. 4. Results The results of our OpenF oam simulations are shown in Fig. 5. Figure 5 shows a comparison between the observed flow around a cylinder originally done by Homann (1936), but then reproduced in Batchelor (1967), and the simulated flow from OpenF oam. For these comparisons, we display cases one through six, and decent agreement is observed. For the R = 32 case, there is a lack of shedding as the fluid propagates past the cylinder. This simulated result is in good agreement with the observed flow. For the R = 55 case, the flow pattern is similar between the simulated and observed cases, though we begin to see evidence of shedding in the simulated case, which is not necessarily observed in the experimental case. As we get to R = 65, we see the flow steepening relative to the earlier – 13 – Fig. 5.—: Left: Results from Fig.1. Right: Results from OpenFoam for cases 1-6. – 14 – Fig. 6.—: Comparison of OpenFoam simulations which both have R = 161, but with different velocities and viscosities. Left: Case 6. Right: Case 7. cases, and more shedding. As the observed case at R = 65 displays shedding, we have better agreement here than for R = 55, but the shedding is still much more significant than should be observed. In the R = 73 through R = 161 cases, we see the flow steepening and the shedding becoming more significant as Reynolds number increases, which agrees with observations. We also wanted to test the law of similarity. For this, two flows with the same Reynolds number should look the same even if their velocities and kinematic viscosities are different. Here we used cases 6 and 7, and the results can be seen in Fig. 6. As may be noted, even with different flow parameters the solutions look exactly the same because they have the same Reynolds numbers. Thus, OpenFoam is able to confirm the law of similarity. While Fig. 5 shows a good qualitative comparison between the observed and simulated results, a quantitative approach is necessary to properly evaluate the effectiveness of OpenF oam. One way to quantitatively describe the flow around a cylinder is via the Strouhal number. The Strouhal number is used to describe the shedding frequency of a fluid and is defined as S= fL , U (13) – 15 – R femp [Hz] Semp fOF [Hz] SOF 55 6.440 0.1288 19.61 0.3922 65 6.995 0.1399 21.31 0.4262 73 7.36 0.1472 22.09 0.4418 102 8.295 0.1659 22.52 0.4504 161 9.340 0.1868 24.26 0.4852 Table 2:: Frequencies and Strouhal numbers for OpenFoam and calculated from an empirical formula from Fey et al. (1998) where f is the shedding frequency. Fey et al. (1998) found an empirical relationship between the Strouhal number and the Reynolds number given by 1.0356 S(R) = 0.2684 − √ , R (14) which can be used to calculate the Strouhal number with only the Reynolds number. By using this empirical relationship to find what the Strouhal number should be for a given case, we can compare the calculated Strouhal number to our simulated Strouhal number to test the accuracy of the simulations. Since we know the diameter of the cylinder and the speed of the fluid for each case (see Table 1), we only need to find the shedding frequency to calculate the Strouhal number. We extract vorticity data from a single point located behind the cylinder in the Karman vortex sheet over the entire time domain. We then approximate the range over which the the fluid has reached a steady state, and we execute a fast-fourier-transform via numpy in Python to calculate the shedding frequency over the period. We use this frequency to calculate the Strouhal number. Likewise, while we can calculate the expected Strouhal number from the Reynolds number, we can also use this Strouhal number to calculate the expected frequency. The results are shown in Table 2. We – 16 – Ri Ri+1 (Si+1 /Si )OF (Si+1 /Si )emp (Si+1 /Si )OF (Si+1 /Si )emp 55 65 1.087 1.087 1.000 65 73 1.037 1.052 0.9857 73 102 1.019 1.127 0.9042 102 161 1.077 1.126 0.9565 Table 3:: Comparison of Strouhal numbers from OpenF oam and from the empirical formula from Fey et al. (1998) did not find the Strouhal number for the case of R = 32 because of the lack of shedding. At first look, the results in Table 2 are not promising. Our simulated frequencies and Strouhal numbers are larger than the expected values by approximately a factor of three. Due to these over-estimations, we also investigated how the Strouhal number changed with Reynolds number. Since the results of Fig. 5 indicate there is an expected amount of change in the shedding as Reynolds number changes, we can investigate if our shedding is changing properly even if the values are off. To do this, we take a ratio of the Strouhal number for a given case with the case directly below it. So, for instance, we find the ratio of R = 65 to R = 55 to test the change in Strouhal number as R changes. We do this for both the simulated and empirical Strouhal numbers, and compare them by taking a ratio of the resultant ratios. The closer this ratio is to one, the more accurate the simulation in terms of how the shedding changes as Reynolds number changes. The results are shown in Table 3, and are very promising. The worst case is the change between R = 73 to R = 102, but the three remaining cases are within 95% of the expected Strouhal ratio from the empirical relationship. At the change from R = 55 to R = 66 we even find the ratios are exactly the same, which is extremely promising. From this, we can see that while OpenFoam is not accurately predicting the Strouhal number, it is accurately modeling the change in the – 17 – Strouhal number. 5. Summary After simulating an incompressible, viscous flow around a cylinder using the icoF oam solver in OpenF oam, we compared our results to observations. Regarding a qualitative comparison, our results showed decent agreement with experimental results from Homann (1936), which were reproduced by Batchelor (1967). We also produced a qualitative comparison of simulations done for the same Reynolds number, but with different fluid parameters, and found good agreement. This showed OpenF oam obeyed the law of similarity. In order to further our analysis, we found the Strouhal number to better understand the shedding frequency in our model. We compared the Strouhal numbers from our simulations to Strouhal numbers calculated from an empirical relationship in Fey et al. (1998) which shows a relationship between the Strouhal and Reynolds numbers. We used this empirical relationship to find an expected value of the Strouhal number and did not find good agreement between our simulations and the expected results. To further our analysis, we explored the change in the Strouhal number as the Reynolds number changes. For this, we found very good agreement. We thus conclude that, using the conditions in this simulation, we are not able to accurately model the shedding in an incompressible, viscous flow around a cylinder, but we are able to accurately model how the shedding should change as Reynolds number changes. We acknowledge that the computational work reported in this paper was performed on the Shared Computing Cluster which is administered by Boston University’s Research Computing Services. We would also like to acknowledge the assistance of Prof. Sheryl – 18 – Grace, Connor Robinson, and Luisa Capannolo in providing assistance when it was needed. Additionally, we would like to acknowledge my penchant for using the royal “we”, even when I am the lone author; some habits are hard to kick. A. blockMeshDict