Modeling the Backward Facing Step using ANSYS FLUENT

Modeling the Backward Facing Step using ANSYS FLUENT ENME 488 August 20, 2021 Professor: Dr. Meilin Yu Advisor: Dr. Kan Liu Prepared by Aryeh Laks 1 Tables ...................................................................................................................................... 2 Figures ..................................................................................................................................... 2 Objectives................................................................................................................................. 4 Introduction.............................................................................................................................. 4 Basic Theories ........................................................................................................................... 4 General................................................................................................................................. 4 Velocity profile ...................................................................................................................... 6 Pressure gradient ................................................................................................................... 7 Geometry of the simulation vs. the geometry of the article .......................................................... 9 Reynolds number of the simulation vs. the Reynolds number of the article....................................11 Discrete phase modeling ........................................................................................................11 Experimental Methods ..............................................................................................................12 Results & Discussion..................................................................................................................20 Conclusion ...............................................................................................................................27 Appendix A: MATLAB functions ...................................................................................................28 Code for inlet velocity profile ..................................................................................................28 Code for reading ASCII Data exported from ANSYS Fluent and graphing it......................................30 Appendix B: Air & Particles (Multiphase Model) ............................................................................32 Appendix C: Pictures .................................................................................................................37 Appendix D: References .............................................................................................................39 Tables Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT [1] ....................... 4 Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS [1] ................................................. 5 Table 3 – Flow properties for the BFS [1] ......................................................................................10 Table 4 – DPM Boundary Conditions [9] .......................................................................................35 Figures Figure 1 – Boundary conditions for the log law ............................................................................... 6 Figure 2 - Velocity profile of the inlet channel using the log law ........................................................ 7 Figure 3 – geometry of the mesh used for ANSYS FLUENT ................................................................ 8 Figure 4 – Origin of the BFS according to Fessler and Eaton [1] ........................................................10 Figure 5 – Geometry of the BFS according to Fessler and Eaton (1999) [1] .........................................10 Figure 6 – Close up view of the mesh used for ANSYS FLUENT .........................................................13 Figure 7 – General settings and gravitational acceleration...............................................................13 2 Figure 8 – Viscous Setting (SST k-omega, Mixture) .........................................................................14 Figure 9 – Model constants ........................................................................................................14 Figure 10 – Properties of air .......................................................................................................15 Figure 11 – Inlet boundary conditions ..........................................................................................15 Figure 12 – Outlet boundary conditions .......................................................................................16 Figure 13 – Solution method, scheme and momentum ...................................................................17 Figure 14 – Solution controls ......................................................................................................17 Figure 15 – Residual monitors.....................................................................................................18 Figure 16 – Initialization settings .................................................................................................19 Figure 17 – Convergence of residuals vs. number of iterations during runtime ...................................19 Figure 18 – Contour plot of x component of velocity ......................................................................20 Figure 19 – Quiver plot of the velocity at the positions 𝑥/𝐻 ≈ −1,2,… ,10........................................21 Figure 20 – Plot of the velocity profile created to match the style of the article..................................21 Figure 21 – Plot of the velocity profiles of the experiment and simulation provided by the article [1] ....22 Figure 22 – Contour plot of x component of velocity with particles ..................................................22 Figure 23 – Particle tracking .......................................................................................................22 Figure 24 – Particle tracking zoomed in ........................................................................................23 Figure 25 – Plot of the velocity profile of the fluid part of the particle phase in the style of the article ...23 Figure 26 – Quiver plot of the velocity profile of the particle part of the particle phase at the positions 𝑥/𝐻 ≈ −1,2,… ,10 ...................................................................................................................24 Figure 27 – Plot of the velocity profile of the particle part of the particle phase in the style of the article ..............................................................................................................................................24 Figure 28 – Particle tracking with no gravity..................................................................................25 Figure 29 – Particle tracking zoomed in with no gravity ..................................................................25 Figure 30 – Plot of the velocity profile of the fluid part of the particle phase with no gravity................25 Figure 31 – Comparison of the velocity profiles of the dispersed phase in the main direction 𝑢𝑃, 𝑥𝑦 behind the BFS from simulations in ANSYS FLUENT and OpenFOAM on mesh A and experiment according to Fessler and Eaton (1999) [1] ...................................................................................................26 Figure 33 – Multiphase model ....................................................................................................32 Figure 34 – Set injection properties .............................................................................................33 Figure 35 – Inert particle properties.............................................................................................34 Figure 36 – Discrete phase model conditions for the inlet boundary conditions .................................34 Figure 37 - Solution Methods......................................................................................................35 Figure 38 – Solution controls ......................................................................................................36 Figure 39 – Multiphase initialization settings.................................................................................36 3 Objectives The objective of this experiment was to simulate the backward facing step (BFS) using ANSYS FLUENT. The intention of this paper is to assist one who is new to ANSYS FLUENT set up a simulation of the backward facing step. Introduction The proper modeling of fluids and particles in fluids has a wide range of applications. An example of one application is the modeling of COVID 19 particles in the air. A simulation of the backward facing step was attempted using ANSYS FLUENT. The reason that the backward facing step was chosen is that it is a welldocumented experiment and can therefore serve as benchmark [1]. This paper is very reliant on an article titled “Assessment of particle-tracking models for dispersed particle-laden flows [1]” for both the assumed value of the numerous parameters and the comparison between the results of the experiment. It will be constantly referred to as “the article” throughout this paper. Basic Theories General The methods, schemes and boundary conditions that are relevant for the simulation are organized in Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT below. Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT [1] Numerical methods steady Solver Pressure-based Pressure–velocity coupling SIMPLE Multigrid Algebraic multigrid (AMG) Under-relaxation factors (correspond to the default settings in ANSYS FLUENT) Spatial discretization • • • 0.3 – pressure 0.7 – momentum 0.8 – turbulence Pressure Standard Convective terms Second-order upwind Gradient Green–Gauss node-based 4 interpolation in all model equations Miscellaneous Second-order upwind The values relevant for the simulation were taken from an article that models the backward facing step using ANSYS FLUENT. They are shown below in Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS. Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS [1] Nodes Cells (with equidistant intervals along the main flow direction x and a relative refinement along the step contour) Dynamic viscosity Properties Reynolds Diameter Density Mass flow rate Particle Injection velocity Maximum velocity (at 𝑦/ℎ = 0.5) Inlet patch location (𝑥/𝐻 ≈ −4 upstream of the step) Mass loading and volume fraction of solid material Boundary and initial conditions Velocity Maximum velocity Inlet Turbulence kinetic energy, k [2] Turbulence dissipation rate, ω [3] Outlet Pressure Wall Velocity 𝑚2 𝜈= 𝑠 𝑘𝑔 𝜇 = 1.838 ∙ 10−5 𝑚∙𝑠 18600 𝑑𝑃 = 7.0 ∙ 10−5 𝑚 𝑘𝑔 𝜌𝑃 = 8800 3 𝑚 𝑘𝑔 𝑚̇ 𝑃 = 1.58 ∙ 10−5 𝑠 𝑚 𝑢 ⃗ 𝑃 = 10.5𝑖 𝑠 𝑚 𝑈𝑥,0 = 10.5𝑖 𝑠 𝑥 𝑃 = −0.1068 𝑚 𝛼𝑃 = 5 ∙ 10−9 𝑚 𝑢 𝑎𝑣𝑔 = 9.39𝑖 𝑠 𝑚 𝑈𝑥,0 = 10.5𝑖 𝑠 𝑚2 𝑘 = 0.45 2 𝑠 1 𝜔 = 2800 𝑠 𝑝 = 0 𝑃𝑎 𝑚 𝑢 ⃗ = 〈0,0,0〉 𝑠 1.5 ∙ 10−5 Kinematic viscosity Air 32000 Although the kinematic viscosity was listed in the article [1] as 1.5 ∙ 10−5 𝑚2 /𝑠, the ANSYS FLUENT software requires the dynamic viscosity. The dynamic viscosity can be calculated using the following formula: 𝜈𝜌 = 𝜇 5 Equation 1 Where 𝜈 is the kinematic viscosity, 𝜌 is the density, and 𝜇 is the dynamic viscosity. Although the density of the air is not given in the article as well, it assumed to be 1.225 kg/m 3 , the preset of the ANSYS FLUENT software for the density of the air. 1.225 kg/m3 is also the density of air at sea level & 15 ⁰C according to the ISA [4]. Therefore, using Equation 1, 1.5 ∙ 10−5 × 1.225 = 𝜇 = 1.838∙ 10−5 𝑘𝑔/(𝑚 ∙ 𝑠). Velocity profile The article [1] also mentions that the simulation was performed with a turbulent profile at the inlet patch with a maximum velocity at the midpoint of the inlet patch and a zero velocity at the edges (walls) of the inlet patch in its simulation. Therefore, the log law was used to model the profile [5] [6]: 𝑦𝑣 ∗ 1 ) + 𝐵] 𝑢 = 𝑣 ∗ [ ln ( 𝜈 𝜅 Equation 2 Where 𝑣 ∗ is the friction velocity, 𝜅 is the Von Karman constant, 𝐵 is a constant, 𝑢 is the mean flow velocity at height 𝑦 above the boundary, and 𝜈 is the kinematic viscosity. In this case 𝜅 = 0.41 and 𝐵 = 5.0. The friction velocity can be solved using the knowledge that the maximum velocity is at the midpoint of the inlet patch. The figure below (see Figure 1 – Boundary conditions for the log law) illustrates. Therefore, 𝑢(𝑦 = ℎ̃) = 𝑢 𝑚𝑎𝑥. Substituting that into Equation 2, one obtains: 1 ℎ̃𝑣 ∗ ) + 𝐵] − 𝑢𝑚𝑎𝑥 = 0 𝑣 ∗ [ ln ( 𝜅 𝜈 Equation 3 0.02𝑣 ∗ 1 ) + 5.0] − 10.5 = 0 𝑣∗ [ ln ( 0.41 1.5 ∙ 10−5 Figure 1 – Boundary conditions for the log law The friction velocity can then be solved using MATLAB’s built-in root solving function (see Code for inlet velocity profile). The friction velocity obtained is 𝑣 ∗ = 0.503. The value for 𝑣 ∗ can be substituted into 6 Equation 2 to obtain the velocity as a function of the position in the inlet. If the values are substituted into Equation 2 and simplified, one obtains: 𝑢 = 1.227ln 𝑦 + 15.30, 0 < 𝑦 ≤ 0.02 Equation 4 This equation is valid for the bottom half of the inlet channel. The upper half of the inlet channel is the mirror image of the lower half. A picture of the velocity profile of the inlet channel was generated by the MATLAB in the appendix (see Code for inlet velocity profile) and is shown below. This MATLAB code also generates an Excel Spreadsheet containing the x and y components of the velocity profile as well as the x and y position coordinates. The information from the Excel Spreadsheet can then be copied and pasted into the Windows notepad application and the text file can be saved as a “.prof” file. At that point, the ANSYS FLUENT program should be able to read it. Figure 2 - Velocity profile of the inlet channel using the log law Pressure gradient The pressure gradient is approximated by using the assumption that the inlet is a laminar, fully developed profile. Therefore, the relation between the velocity and the pressure gradient can be modeled by the following equation: 𝑢= 1 𝜕𝑝 2 ( ) 𝑦 + 𝑐1 𝑦 + 𝑐2 2𝜇 𝜕𝑥 7 Equation 5 One can see that the equation (Equation 5) is a quadratic equation with a parabolic shape. As the article provides three boundary conditions, it is easiest to substitute 1 𝜕𝑝 ( ) in place of the constant 𝑐0 and use 2𝜇 𝜕𝑥 the following equations when determining the velocity profile of the inlet patch: 𝑢 = 𝑐0 𝑦2 + 𝑐1 𝑦 + 𝑐2 Equation 6 Where 𝑢 is the velocity in the x direction, 𝑦 is the position along the y axis, and 𝑐0, 𝑐1 , and 𝑐2 are constants. The value of the constants can be determined using the following boundary conditions: 𝑢(𝑦 = 0.0267 𝑚) = 0 𝑢(𝑦 = 0.0467 𝑚) = 10.5 𝑢(𝑦 = 0.0667 𝑚) = 0 𝑚 𝑠 Equation 7 Where 𝑦 = .0267 𝑚 is the location of the wall on the bottom of the inlet patch, 𝑦 = .0467 𝑚 is the location of the midpoint of the inlet patch, and 𝑦 = .0667 𝑚 is the location of the wall on the top of the inlet patch. The figure below (Figure 3 – geometry of the mesh used for ANSYS FLUENT) may be useful in understanding the relationship in the boundary conditions between the velocity in the x direction and the position along the y axis. The figure also describes the mesh arrangement used in the simulation. A larger picture of figure is shown in the Appendix as well (see Appendix C: Pictures). Figure 3 – geometry of the mesh used for ANSYS FLUENT Substituting the boundary conditions (Equation 7) into the equation for steady, laminar flow between fixed parallel plates (Equation 5), one obtains: 0 = (0.0267)2 𝑐0 + (0.0267)𝑐1 + 𝑐2 10.5 = (0.0467)2 𝑐0 + (0.0467)𝑐1 + 𝑐2 0 = (0.0667)2 𝑐0 + (0.0667)𝑐1 + 𝑐2 Equation 8 Rearranging into matrix form, one procures: (0.0267)2 [(0.0467)2 (0.0667)2 0.0267 1 𝑐0 0 0.0467 1] [𝑐1 ] = [10.5] 0 0.0667 1 𝑐2 8 Equation 9 𝐴𝐶 = 𝑈 Consequently, if one solves the matrix using inversion (𝐶 = 𝐴−1 𝑈 or in MATLAB: C = A\U), one can solve the constants in the equation for steady, laminar flow between fixed parallel plates (Equation 5), yielding: 𝑢(𝑦) = −26250𝑦2 + 2452𝑦 − 46.75 Equation 10 Using Equation 5 and Equation 6, one can see that: 𝑐0 = 1 ∆𝑝 1 𝜕𝑝 ( )≈ 2𝜇 ∆𝑥 2𝜇 𝜕𝑥 Equation 11 Rearranging, one obtains an approximation of the pressure gradient: ∆𝑝 = 𝑝1 − 𝑝2 = −2𝜇𝑐0 ∆𝑥 Equation 12 Where ∆𝑝 is the difference of pressure between two points (points 1 and 2 with point 1 being upstream), 𝑐0 is a constant (-26250), 𝜇 is the dynamic viscosity, and ∆𝑥 is the distance between two points (points 1 and 2 with point 1 being upstream). As this simulation uses a longer inlet with the distance from the inlet to the step being 1.4 m, ∆𝑥 is assumed to be 1.4 m. Therefore, if the values are substituted in, one obtains −2(1.838 ∙ 10−5 )(−26250)(1.4) = 𝑝1 − 𝑝2 = 1.35 − 0 (see Code for inlet velocity profile). Consequently, a pressure gradient of 1.35 Pa was used for the simulation. Geometry of the simulation vs. the geometry of the article Incidentally, it is important to note that while the origins of the x and y axis of the mesh used for ANSYS FLUENT match up to coordinates assumed by the article (see Figure 4 – Origin of the BFS according to Fessler and Eaton ), the inlet channel of the BFS is lengthened from 0.1 m (see Figure 5 – Geometry of the BFS according to Fessler and Eaton (1999) ) to 1.4 m (Figure 3 – geometry of the mesh used for ANSYS FLUENT) 9 Figure 4 – Origin of the BFS according to Fessler and Eaton [1] Figure 5 – Geometry of the BFS according to Fessler and Eaton (1999) [1] Moreover, the article says that the turbulent profile at the inlet patch is located at 𝑥/𝐻 ≈ −4 upstream of the step, where H is the step height at 0.0267 m, shown in the table below (Table 3 – Flow properties for the BFS [1]Table 3 ). Solving for x, one obtains 0.1068 m upstream of the step. This translates into 𝑥 = −0.1068 𝑚 for the location of the turbulent profile at the inlet patch on the x axis. However, as the inlet was lengthened, the inlet patch of the long inlet, 𝑥 = −1.4 𝑚 was used whenever the location of the turbulent profile at the inlet patch on the x axis was applicable. Table 3 – Flow properties for the BFS [1] Channel flow BFS flow Channel height, h 40 mm Step height, H 26.7 mm Channel width, w 457 mm w/H 17:1 h/H 5:3 Additionally, the article mentions that the number of nodes used in the simulation was 32000. The number of nodes for the lengthened inlet used in this simulation was 76,054. The reason the inlet was lengthened was to even out any irregularities arising from the rough modeling of the velocity profile using the log law before reaching the backwards facing step. 10 Reynolds number of the simulation vs. the Reynolds number of the article Another interesting detail to note was the way the article calculated the Reynolds number. The classic way to calculate the Reynolds number in a pipe or channel is: 𝑅𝑒 = 𝜌𝑈𝐷 𝑈𝐷 = 𝜇 𝜈 Equation 13 Where 𝑅𝑒 is the Reynolds number, 𝜌 is the density, 𝜇 is the dynamic viscosity, 𝜈 is the kinematic viscosity, 𝑈 is the velocity, and 𝐷 is the pipe or hydraulic diameter. The way the article calculates the Reynolds number is: 𝑅𝑒 = 𝑈𝑥,0 𝐻 𝜈 Equation 14 Where 𝜈 is the kinematic viscosity found in Table 2 – Values extracted from article (Sec 4.1-2) modeling 𝑚 BFS , 𝑈𝑥,0 = 10.5 is the maximum velocity also found in Table 2 – Values extracted from article (Sec 4.1𝑠 2) modeling BFS , and 𝐻 = 0.0267 𝑚 is the step height found in Table 3 – Flow properties for the BFS [1]. When those values are substituted into Equation 14, one obtains 18690, approximately equal to the value that the article gives for the Reynolds number (18600). Interestingly enough, the article chooses the step height despite fact that the step height does not have anything to do with the BFS’s inlet height or outlet height. [1] The more logical way to calculate the Reynolds number is: 𝑅𝑒𝑎𝑣 = 𝑢𝑎𝑣 𝐷ℎ 𝑢 𝑎𝑣(2ℎ) = 𝜈 𝜈 Equation 15 Where 𝑅𝑒𝑎𝑣 is the average Reynolds number, 𝜈 is the kinematic viscosity, 𝑢 𝑎𝑣 is the average velocity, 𝐷ℎ is the hydraulic diameter, and ℎ is the channel height. [5] The average velocity can be determined through the following equation: 𝑢 𝑎𝑣 = 0.5ℎ 1 ∫ 𝑢(𝑦) 𝑑𝑦 0.5ℎ 0 Equation 16 The average velocity was calculated numerically (see Code for inlet velocity profile) as 9.273 m/s. Incidentally, this is close to the average velocity of 9.39 m/s of the article. Discrete phase modeling 11 The discrete phase is modeled numerically by using a methodology called the Lagrangian–Eulerian approach. The fluid is assumed to be isothermal and incompressible while the particles are assumed to be spherical and point masses, allowing for the torque to be omitted. Furthermore, the particle volume fraction is usually assumed negligible compared to the carrier phase volume and particle-particle interactions are usually neglected as well. The particle motion is solved by summing the force balance, which is written in a Lagrangian frame on the particles. The Navier–Stokes equations are solved for the continuous carrier phase (using Reynolds-averaged Navier-Stokes – RANS) and the particles are treated as a dispersed phase and are tracked individually. Particle boundary layers, wake flow regions and similar flow processes on the particle scale or smaller are not directly solved and are instead modeled via functional correlations with empirical parameters. [1] [7] Experimental Methods The purpose of the Experimental Methods section is to detail the setup of the ANSYS FLUENT software when air without particles is modeled in the BFS geometry (one phase). It is also the intention that this section should be able to act as a tutorial for one who is new to ANSYS FLUENT. See Appendix B: Air & Particles for when air being modeled with particles in the discrete phase model (two phases) using ANSYS FLUENT. A picture of the close up of the mesh is shown below (see Figure 6 – Close up view of the mesh used for ANSYS FLUENT). The dimensions of the mesh were adapted from the article [1]. The figure taken from the article below (Figure 5 – Geometry of the BFS according to Fessler and Eaton (1999) ) provides the proportions. Earlier, it was shown that the channel height, h, is 40 mm (Table 3 – Flow properties for the BFS [1]). Therefore, 𝑥 𝐿,𝑈 = 2.5 ∙ 40 𝑚𝑚 = 0.1 𝑚 and 𝑥 𝐿,𝑈 = 35 ∙ 40 𝑚𝑚 = 1.4 𝑚. However, as mentioned earlier (see Figure 3 – geometry of the mesh used for ANSYS FLUENT), 𝑥 𝐿,𝑈 was changed to 1.4 m. 12 Figure 6 – Close up view of the mesh used for ANSYS FLUENT In the general parameters, the model’s gravitational acceleration was set to -9.8 m/s2 in the y direction (see Figure 8 – Viscous Setting (SST k-omega, Mixture)). Figure 7 – General settings and gravitational acceleration The model was set to Viscous (SST k-omega, Mixture) (see Figure 8 – Viscous Setting (SST k-omega, Mixture)) and the Model constants were left unchanged (see Figure 9 – Model constants). 13 Figure 8 – Viscous Setting (SST k-omega, Mixture) Figure 9 – Model constants The viscosity of the air was changed to 1.838 ∙ 10−5 kg/m∙s while the density of the air was left as 1.225 kg/m3 (see Figure 10 – Properties of air). 14 Figure 10 – Properties of air The velocity profile was imported. Using the MATLAB code (see Code for inlet velocity profile), there should be 4 components: x (a uniform position of -1.4 m), y (position in m), x-velocity (the velocity profile generated by the log law in m/s), and y-velocity (a uniform velocity of 0 m/s). For the boundary condition of the inlet, the type was set to Velocity Inlet. The pressure was set to 1.35 Pa, the x component of the velocity was set to x velocity from the imported velocity profile, and the y component of the velocity was set to y velocity from the imported velocity profile. (see Figure 11 – Inlet boundary condition). Figure 11 – Inlet boundary conditions 15 For the boundary condition of the outlet, the type was set to Pressure Inlet. The pressure was set to 0 (Figure 12 – Outlet boundary conditions below). The turbulence specification method was set to K and Omega for both the inlet and the outlet. The turbulence kinetic energy was set to k = 0.45 m2 /s2 and the turbulence dissipation rate was set to ω = 2800 1/s to match the article (see Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS earlier as well as Figure 11 – Inlet boundary conditions and Figure 12 – Outlet boundary conditions). Figure 12 – Outlet boundary conditions For the solution method, the scheme was set to SIMPLE, the gradient was set to Green–Gauss node-based, and the Pressure was set to Standard to match the article (see Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT earlier) and the other schemes were assumed to be set to SecondOrder Upwind (see Figure 13 – Solution method, scheme and momentum below). 16 Figure 13 – Solution method, scheme and momentum For the controls, the under-relaxation factors were left unchanged (see Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT). However, the advanced solution controls were modified. The Stabilization method was set to GMRES. For the Fixed Cycle Parameters in the Algebraic multigrid (AMG), the number of Pre-Sweeps and Post-Sweeps was set to 10 and 1, respectively. Figure 14 – Solution controls 17 The absolute criteria of the residual monitors in the monitors section were all set to 10−5 (see Figure 15 – Residual monitors). Figure 15 – Residual monitors For the initialization, the gauge pressure was initialized to 1.35 Pa. The x velocity was initialized to 9.39 m/s, the mean velocity of the fully developed velocity profile of the inlet patch, and the y velocity was set to 0 m/s as found in the article [1] (see Figure 16 – Initialization settings below). The turbulence kinetic energy, k [2], was set to 0.45 m2 /s2 and the turbulence dissipation rate, ω [3], was set to 2800 1/s as found in the article (see Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS earlier). The initialization was computed from the inlet. 18 Figure 16 – Initialization settings The number of iterations was set to 10000. However, the residuals converged before 2500 iterations (see Figure 17 – Convergence of residuals vs. number of iterations during runtime). Figure 17 – Convergence of residuals vs. number of iterations during runtime 19 Results & Discussion The velocity was first calculated without including the particle phase. The x component of the velocity is shown in the figure below (see Figure 18 – Contour plot of x component of velocity). Incidentally, one can see that the “.prof” file was properly read as the inlet channel has a higher velocity in the middle of the channel, as opposed to the walls of the inlet which have a lower velocity. One can also see backflow with fluid moving towards the negative x direction (the blue colored region) directly behind the BFS as one would expect. Another positive sign of the simulation conforming to reality is that the velocity drops to around zero m/s along the walls further behind the BFS (the blue colored region). Figure 18 – Contour plot of x component of velocity The x velocity solution data was then exported as an ASCII file (with commas included) and a program was designed in MATLAB to import the ASCII file and convert it into standard MATLAB data (see Code for reading ASCII Data exported from ANSYS Fluent and graphing it). The data was then processed to plot the velocity profiles at 𝑥/𝐻 ≈ 𝑛 = −1,2,… ,10 (see Figure 19 – Quiver plot of the velocity at the positions 𝑥/𝐻 ≈ −1,2,… ,10). One can see that the reattachment happens at 𝑥/𝐻 ≈ 7 as specified by the article [1]. Incidentally, one can see that there are a larger number of data points towards the bottom of the wider channel. This is due to the design of the BFS used in this simulation. As mentioned earlier (Figure 6 – Close up view of the mesh used for ANSYS FLUENT, see also see Appendix C: Pictures), the BSF was designed with a smaller mesh size towards the bottom of the BFS; hence a higher concentration of points in that location. 20 Figure 19 – Quiver plot of the velocity at the positions 𝑥/𝐻 ≈ −1,2,… ,10 Another figure was generated using the MATLAB code in order to compare the results of this simulation to the experiment and modeling provided by the article (see Figure 20 – Plot of the velocity profile created to match the style of the article below). For the y axis, the y position was scaled by dividing by the step height, 𝐻. Therefore, the y axis is unitless. The graph of the x axis was created by adding the unitless, scaled x velocity (2 ∙ 𝑢 𝑥/𝑈𝑚𝑎𝑥) to the unitless, scaled x position (𝑥/𝐻). (See Code for reading ASCII Data exported from ANSYS Fluent and graphing it) Figure 20 – Plot of the velocity profile created to match the style of the article The simulation of this experiment matches very well both in the profile shape of the velocity and in the reattachment point to the simulation of the article. The figure provided by the article is shown below (Figure 21 – Plot of the velocity profiles of the experiment and simulation provided by the article ). 21 Figure 21 – Plot of the velocity profiles of the experiment and simulation provided by the article [1] The velocity was then calculated with the particle phase. The x component of the velocity is shown in the figure below (Figure 22 – Contour plot of x component of velocity with particles). Figure 22 – Contour plot of x component of velocity with particles The velocity magnitude as well as the path of the particle streams through the fluid are shown in the figure below (Figure 23 – Particle tracking). A larger picture of figure is shown in the Appendix as well (see Appendix C: Pictures). A figure that magnifies on the particle streams by a slight amount is also shown below (Figure 24 – Particle tracking zoomed in). Figure 23 – Particle tracking 22 Figure 24 – Particle tracking zoomed in Some interesting observations arise from the figures above. One observation is that the gravity exerted on the particles affects the trajectory of the particles. One can see (Figure 23 – Particle tracking and Figure 24 – Particle tracking zoomed in) that the particles eventually sink towards the floor of the inlet before the reaching BFS as well as sinking towards the floor along the channel downstream of the BFS upon exiting the inlet channel. Another observation is that the trajectory of the particles appears to have no effect on the flow of the fluid as well. One cannot see a noticeable difference between the simulation with the Discrete Phase Model (Figure 22 – Contour plot of x component of velocity with particles) and the simulation of the fluid without particles (Figure 18 – Contour plot of x component of velocity) earlier. It appears that although the particles are affected by gravity, they do not affect the fluid motion of the DPM (Figure 22 – Contour plot of x component of velocity with particles) when compared to the simulation of the fluid without particles (Figure 18 – Contour plot of x component of velocity). Another observation is the particles appear to bounce after contacting with the floor of the channel downstream of the BFS (see Figure 23 – Particle tracking and Figure 24 – Particle tracking zoomed in). This is likely due to the discrete phase model conditions of the walls. Since the reflect condition was selected for the walls, the particles bounce off the floor with a reduced momentum and velocity in the x direction and a reduced velocity in the opposite of the previous direction along the y axis [8]. A further observation is that although the air still exhibits backflow, none of the particles themselves appear to have backflow. The solution data of the particle phase’s fluid velocity was then exported as an ASCII file to MATLAB. One can see that there is no visible difference between the MATLAB plot of the DPM (Figure 25 – Plot of the velocity profile of the fluid part of the particle phase in the style of the article) and the MATLAB plot of the fluid without particles (Figure 20 – Plot of the velocity profile created to match the style of the article). One can see that the reattachment happens at 𝑥/𝐻 ≈ 7 as expected [1]. Figure 25 – Plot of the velocity profile of the fluid part of the particle phase in the style of the article 23 The solution data of the particle phase’s particle velocity was also exported as an ASCII file to MATLAB. It appears from the two figures below (Figure 26 – Quiver plot of the velocity profile of the particle part of the particle phase at the positions 𝑥/𝐻 ≈ −1,2,… ,10 and Figure 27 – Plot of the velocity profile of the particle part of the particle phase in the style of the article) that the x velocity of the particle phase did not get exported properly by ANSYS. It could be because the particles did not completely fill the space of the BFS geometry. Figure 26 – Quiver plot of the velocity profile of the particle part of the particle phase at the positions 𝑥/𝐻 ≈ −1,2,… ,10 Figure 27 – Plot of the velocity profile of the particle part of the particle phase in the style of the article The experiment was also simulated without gravity. One can see below (Figure 27 – Plot of the velocity profile of the particle part of the particle phase in the style of the article and Figure 28 – Particle tracking with no gravity) that the particles don’t sink towards the bottom of the channels and instead spread out 24 downstream after the step. However, it is apparent that there is no backflow, and therefore the particles do not completely fill the space. Consequentially, it is logical that the solution data of the particle phase’s particle velocity did not graph properly and resembled the two figures of the particle phase when gravity was present (Figure 26 – Quiver plot of the velocity profile of the particle part of the particle phase at the positions 𝑥/𝐻 ≈ −1,2,… ,10 and Figure 27 – Plot of the velocity profile of the particle part of the particle phase in the style of the article earlier) Figure 28 – Particle tracking with no gravity Figure 29 – Particle tracking zoomed in with no gravity As expected, there is no difference between the fluid velocity when there is no gravity (see Figure 30 – Plot of the velocity profile of the fluid part of the particle phase with no gravity) when compared with the fluid velocity when there gravity is present, as the particles don’t affect the fluid flow. Figure 30 – Plot of the velocity profile of the fluid part of the particle phase with no gravity The simulation performed does well when compared to the experimental data of the article. One can see from the figure below (Figure 31 – Comparison of the velocity profiles of the dispersed phase in the main direction 𝑢 𝑃,𝑥 (𝑦) behind the BFS from simulations in ANSYS FLUENT and OpenFOAM on mesh A and 25 experiment according to Fessler and Eaton (1999)) that the fluid flow of the particle phase (see Figure 30 – Plot of the velocity profile of the fluid part of the particle phase with no gravity) shares the same geometry as the experimental data of Fessler and Eaton (1999) found in the article [1]. However, the simulation found in this report differs with Mesh A (A two-dimensional mesh that most resembles the mesh used in this experiment) found in the figure below. Figure 31 – Comparison of the velocity profiles of the dispersed phase in the main direction 𝑢𝑃,𝑥 (𝑦) behind the BFS from simulations in ANSYS FLUENT and OpenFOAM on mesh A and experiment according to Fessler and Eaton (1999) [1] However, the article [1] points out that the block profile distribution of the fluid velocity found in Mesh A is not realistic. The particle cloud should spread span-wise to the flow and disperse into the recirculation zone as indicated by the experiment performed by Fessler and Eaton. Therefore, the article concludes that when turbulent particle laden flow is simulated, ANSYS FLUENT simulations should performed using three-dimensions even if it is reasonable to assume a two-dimensional fluid flow structure. The article suggests that this can be due to the three-dimensional character of turbulent dispersion getting lost in a two-dimensional fluid flow structure. Although it does not appear necessary, it would be interesting to perform the simulation using a three-dimensional mesh to see if there are any minor differences between the simulation performed in this paper and a simulation using three dimensions. The article also suggests that there may be a bug in the two-dimensional dispersion model of ANSYS FLUENT. This can explain why there is no error in the simulation as this experiment was performed using ANSYS 2020 and ANSYS 2021 while the article was published in 2016. Even if the article is correct in its suspicion of a bug, it is possible that ANSYS fixed the bug between 2016 and 2020. 26 Conclusion A simulation of the backward facing step was attempted using ANSYS FLUENT and compared with an article titled “Assessment of particle-tracking models for dispersed particle-laden flows” [1]. The BFS was modeled both with particles and without particles. The results of this simulation agree with those from the literature. The outcome of Fessler and Eaton’s experiment and this simulation appears to have no visible difference in the geometry of the velocity profile regardless of whether there were particles injected into the airflow or not. The reattachment point of the flow matched the article as well. Although the article points possible concerns when assuming a two-dimensional fluid flow structure for a turbulent particle laden flow, using a two-dimensional mesh did not produce any noticeable differences between the simulation in this experiment and the experiment of Fessler and Eaton. This experiment was instrumental in allowing me to learn the basics of using ANSYS FLUENT software. Proper knowledge of how to use computational software is very important in many fields of engineering. 27 Appendix A: MATLAB functions Code for inlet velocity profile clc, close all, clear all, format compact %% Given Values h = 40/1000; H = 26.7/1000; uWall = 0; uMax = 10.5; % % % % channel height, m step height, m velocity at the walls of the inlet patch, m/s maximum velocity, m/s yLow = H; % location of the bottom wall of the inlet patch, m yMid = H + .5*h; % location of the maximum velocity of the inlet patch, m yUpp = H + h; % location of the top wall of the inlet patch, m xLoc = -1.4; dx = 1.4; % location of the velocity profile of the inlet patch along the x-axis, m % distance for approximating pressure gradient nu = 1.5e-5; kappa = 0.41; B = 5.0; n = 100; rho = 1.225; ReArt = 18600; % % % % % % kinematic viscosity for power law, m^2/s Von Karman constant for power law constant for power law number of points density, kg/m^3 Reynold's number from the article, unitless z = 0; w = 0; dP = 7e-5; T = 273 + 25; mPdot = 1.58e-5; % % % % % position in z direction velocity in z direction particle diameter, m fluid temperature, assumed to be standard mass flow rate, kg/s x = ones(n,1)*xLoc; % x-position, m y = linspace(H,H+h,n)'; % y-position, m %% Solving the Boundary Conditions and Creating Velocity Profile - Power Law vstar0 = @(vstar) vstar*(1/kappa*log(h/2*vstar/nu) + B) - uMax; % friction velocity equation vstar = fzero(vstar0,0.4); % friction velocity, m/s % x-velocity equation, m/s uxPowerLaw = @(y) vstar.*(log(y.*vstar./nu)/kappa + B); u = zeros(n,1); m/s u(1:n/2) = uxPowerLaw( y(1:n/2)-H ); u(n/2+1:end) = uxPowerLaw( h-y(n/2+1:end)+H ); u(1) = uWall; u(n/2) = uMax; u(end) = uWall; % preallocation for x-velocity, % % % % % lower half of inlet, m/s upper half of inlet, m/s BC no slip condition, m/s BC mid channel condition, m/s BC no slip condition, m/s %% Find Reynolds's Number uAvg = ((h/2)^-1)*integral(uxPowerLaw,0,h/2); ReAvg = uAvg*(2*h)/nu; ReArtTry = uMax*H/nu; nuTry = (2*h)*uAvg/ReArt; number, m^2/s mu = nu * rho; muTry = nuTry * rho; number, kg/m*s % average velocity, m/s % expected Reynold's number % how the article calculated the Reynold's number % kinematic viscosity using article's Reynold's % expected dynamic viscosity , kg/m*s % dynamic viscosity using article's Reynold's fprintf('The Reynolds number using the method described in the article is %i\n', ReArtTry) fprintf('The Reynolds number using the method described by the professor is %.3e\n', ReAvg) %% Solving the Boundary Conditions and Creating Velocity Profile - Parabolic, Calculating Pressure Gradient 28 A = [yLow^2,yLow,1;... yMid^2,yMid,1;... yUpp^2,yUpp,1]; U = [uWall;uMax;uWall]; c = A\U; c0 = c(1); c1 = c(2); c2 = c(3); % % % % % c0*yLow^2 + c1*yLow c0*yMid^2 + c1*yMid c0*yUpp^2 + c1*yUpp velocity BCs, m/s constants generated + c2 = uWall, m/s + c2 = uMax, m/s + c2 = uWall, m/s by BCs dp = -2*mu*c0*dx; fprintf('\nThe approximate pressure gradient for a distance of %.3f m is %.3f Pa\n', dx,dp) % uxParabolic = @(y) c0.*y.^2 + c1.*y + c2; % x-velocity equation, m/s % u = uxParabolic(y); % *3/2 % x-velocity points, m/s %% Plot Velocity Profile v = zeros(n,1); % y-velocity, m/s quiver(x,y,u,v,.5) xlabel('x location (m)') ylabel('y location (m)') title('Velocity Components') legend('Inlet patch velocity profile') %% Write Velocity Profile to Excel titleExcel = ['((inlet point ',num2str(n),')']; C1 = {titleExcel;'(x'}; C2 = num2cell(x); C3 = {')'}; C4 = {'';'(y'}; C5 = num2cell(y); C6 = {')'}; C7 = {'';'(x-velocity'}; C8 = num2cell(u); C9 = {')'}; C10 = {'';'(y-velocity'}; C11 = num2cell(v); C12 = {')';'';')'}; C = [C1;C2;C3; C4;C5;C6; C7;C8;C9; C10;C11;C12]; % The information on the Excel Spreadsheet should be copied & pasted to a % txt file and saved as a .prof file filename = 'xvelocitylonginlet.xlsx'; writecell(C,filename,'Sheet','xvelocitylonginlet'); fprintf('\nThe velocity profile for the long inlet is exported to an Excel Spreadsheet named %s\n', filename) %% Write Injection Profile to Excel % ((x y z u v w diameter temperature mass-flow-rate) streamID) name = 'xinjectionlonginlet'; filename = [name,'.xlsx']; % C1 = {titleExcel} for i = 1:length(x) Cinj(i,1) = cellstr(['((',num2str(x(i)),' ',num2str(y(i)),' ',num2str(z),' ',... num2str(u(i)),' ',num2str(v(i)),' ',num2str(w),' ',... num2str(dP),' ',num2str(T),' ',num2str(mPdot),') streamID:',num2str(i),')']); end writecell(Cinj,filename,'Sheet',name); fprintf('\nThe injection profile for the long inlet is exported to an Excel Spreadsheet named %s\n', filename) 29 Code for reading ASCII Data exported from ANSYS Fluent and graphing it clc, close all, clear all, format compact tic %% Given Values h = 40/1000; H = 26.7/1000; xin = -1; xout = 14; a = -1; b = 10; % % % % % % channel height, m step height, m location of the inlet, m location of the outlet, m start of quiver graphing interval, m/m end of quiver graphing interval, m/m %% Load and Preprocess ASCII File % filename = 'OutputXVelLongInletNoPartFinal'; filename = 'OutputXVelPhase1PartNoGinj'; premat = readcell(filename); premat(:,1) = []; % remove column of ANSYS cell numbers strCell = premat(1,:); % variable names mat = cell2mat(premat(2:end,:)); % converts cell structure to matrix xPos = mat(:,1); % x position, m yPos = mat(:,2); % y position, m xVel = mat(:,3); % x velocity, m/s yVel = zeros(size(mat,1),1); % y velocity, m/s xdH = xPos/H; xdHsimp = round(xdH*1,1); % x/H, x scaled by H to match the article % keep rounded points of x/H %% Channel Shape and Quiver Plot hold on xlabel('x/H (m/m)') ylabel('y (m)') % shape of long inlet x = [0 xout xout xin xin 0 0]; % m y = [0 0 h+H h+H H H 0]; % m patch(x,y,'w') % quiver plot for i = a:b Irange = find(xdHsimp == i*1); xPquiv = xdH(Irange); yPquiv = yPos(Irange); xVquiv = xVel(Irange); yVquiv = yVel(Irange); q = quiver(xPquiv,yPquiv,xVquiv,yVquiv,200); q.ShowArrowHead = 'off'; q.Marker = '.'; end % % % % % % index for quiver column x position for quiver, m/m y position for quiver, m x velocity for quiver, m/s y velocity for quiver, m/s quiver plot column % mark base of the quiver arrow %% Plot figure hold on axis equal xlabel('2 \cdot u_x/U_{max} + x/H') ylabel('y/H') patch(x,y/H,'w','LineWidth',2) % shape of long inlet scaled by H axis([xin xout 0 h/H+1]) % axis limit of graph scaled by H Uscale = 2/10.5; % 2/Umax for i = [2,5,7,9,12] Irange = find(xdHsimp == i*1); % index for column line([i i],[0 h/H+1]) % y position line yP = yPos(Irange)/H; % y position for x velocity xV = xVel(Irange)*Uscale + i; % x velocity according to x position 30 yPxV = sortrows([yP,xV]); plot(yPxV(:,2),yPxV(:,1),'k') % matrix for sorting velocity according to y position % plot(x velocity + position, y position) end toc 31 Appendix B: Air & Particles (Multiphase Model) The purpose of this section is to detail the setup of the ANSYS FLUENT software when air is modeled with particles in the BFS geometry (multiphase phase). It is also the intention that this section should be able to act as a tutorial for one who is new to ANSYS FLUENT. In the models section, the multiphase model was switched from off to mixture (see Figure 32 – Multiphase model below). Figure 32 – Multiphase model An injection file was created using MATLAB (see Code for inlet velocity profile). The velocity profile used to create the injection velocity profile was the same as the velocity profile of the inlet fluid (see Velocity profile). The diameter and flow rate of both points are set to 7e-5 m and 1.58e-5 kg/s respectively. The position and velocity in the z direction were assumed to be 0 and the temperature was assumed to be 25 ⁰C or 298 Kelvin (in retrospect, it would have been better to use 15 ⁰C or 278 Kelvin to stay consistent with the density according to the ISA standards [4]). The values for the diameter and flow rate were obtained directly from the article and collected in Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS. The MATLAB code generated an Excel Spreadsheet which was then copied and pasted to notepad and saved as an “.inj” file. The file was then imported using the Set injection properties found in the Discrete 32 Phase section under the models section (see Figure 33 – Set injection properties). The injection type was set to file and the file was selected using the file button. Incidentally, the material was set to copper. The name of material can be left unchanged as well so long as the material is set to the correct density. Figure 33 – Set injection properties The density of the inert particles in the materials section was changed to 8800 kg/m3 (see Figure 34 – Inert particle properties below) to match the density of the copper particles mentioned in the article and collected in Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS. 33 Figure 34 – Inert particle properties For the discrete phase model conditions in the boundary conditions section (DPM tab in Figure 35 – Discrete phase model conditions for the inlet boundary conditions below), the escape condition was selected for the inlet and outlet and the reflect condition was selected for the walls. The phase was also set to mixture. Figure 35 – Discrete phase model conditions for the inlet boundary conditions 34 For more elaboration of the discrete phase model conditions, a table has been included below (Table 4 – DPM Boundary Conditions ). More information can be found online at ANSYS FLUENT 12.0 User's Guide 23.4.1 Discrete Phase Boundary Condition Types [8]. Table 4 – DPM Boundary Conditions [9] Boundary DPM BC Physical Meaning inlet escape - outlet escape Particle leaves the domain – Tracking ends wall trap Particle is removed but its current mass and energy is imparted to the gas phase (if the particles will stick to this boundary). wall_pipe-[*] reflect Particle rebounds off of the pipe wall with the normal and tangential coefficients of restitution set in the reflect panel. For Methods under the Solutions section, as found in Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT, the scheme was set to Coupled, the gradient was set to Green-Gauss Node Based, and all Spatial discretization settings besides the volume fraction were set to Second Order Upwind (see Figure 36 - Solution Methods below) [9]. Figure 36 - Solution Methods For the Pseudo Transient Explicit Relaxation Factors in the controls, the Momentum, Volume Fraction, and Turbulent Kinetic Energy were changed to 0.7, 5 ∙ 10−9 , and 0.8 respectively (see Table 1 – Numerical methods, schemes and boundary conditions in ANSYS FLUENT). In the advanced solution controls, the Stabilization method were all set to GMRES. For the Fixed Cycle Parameters in both the Scalar Parameters 35 and the Coupled Parameters of the Algebraic multigrid (AMG), the number of Pre-Sweeps and PostSweeps were set to 10 and 1, respectively (see Figure 37 – Solution controls). Figure 37 – Solution controls For initialization under the Solutions section, the phase 2 volume fraction to was set to 5 ∙ 10−9 as found in the article and collected in Table 2 – Values extracted from article (Sec 4.1-2) modeling BFS. (see Figure 38 – Multiphase initialization settings below) Figure 38 – Multiphase initialization settings 36 Appendix C: Pictures 37 38 Appendix D: References [1] F. Greifzu, C. Kratzsch, T. Forgber, F. Lindner and R. Schwarze, "Assessment of particle-tracking models for dispersed particle-laden flows," ENGINEERING APPLICATIONS OF COMPUTATIONAL FLUID MECHANICS, 2016. [2] CFD Online, "Turbulence kinetic energy," CFD-Wiki, the free CFD reference, 06 2011. [Online]. Available: https://www.cfd-online.com/Wiki/Turbulence_kinetic_energy. [Accessed May 2021]. [3] CFD Online, "Specific turbulence dissipation rate," CFD-Wiki, the free CFD reference, 06 2011. [Online]. Available: https://www.cfd-online.com/Wiki/Specific_turbulence_dissipation_rate. [Accessed May 2021]. [4] A. M. Helmenstine, "What Is the Density of Air at STP?," ThoughtCo.com, 02 2020. [Online]. Available: https://www.thoughtco.com/density-of-air-at-stp-607546. [Accessed May 2021]. [5] D. M. Yu, Log_law_turbulence.pdf, Baltimore, 2021. [6] F. Stern, "Chapter 8 Flow in Conduits," Fall 2014. [Online]. Available: http://user.engineering.uiowa.edu/~fluids/posting/Lecture_Notes/Chapter8.pdf. [Accessed August 2021]. [7] P. S. SHANKARA, "CFD SIMULATION AND ANALYSIS OF PARTICULATE DEPOSITION ON GAS | OhioLINK Electronic Theses & Dissertations (ETD) Center," 2010. [Online]. Available: https://etd.ohiolink.edu/apexprod/rws_etd/send_file/send?accession=osu1262290700&dispositio n=attachment. [Accessed August 2021]. [8] ANSYS, "ANSYS FLUENT 12.0 User's Guide - 23.4.1 Discrete Phase Boundary Condition Types," ENEAGRID PROJECTS WEB PAGES, 01 2009. [Online]. Available: https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/node699.htm. [Accessed May 2021]. [9] ANSYS, Introduction to ANSYS FLUENT - Workshop 6 - Using the Discrete Phase Model (DPM), ANSYS, 2010. [10] P. M. Gerhart, A. L. Gerhart and J. I. Hochstein, "6.9.1 Steady, Laminar Flow between Fixed Parallel Plates," in Munson, Young, and Okiishi’s Fundamentals of Fluid Mechanics, United States of America, John Wiley & Sons, 2016, p. 323. 39