Practical EG501V Computational Fluid Dynamics - 2023/2024
Practical 2. CFD of turbulent pipe flow and heat transfer
1 Introduction
Building on the knowledge of our previous practical of modelling laminar flow through a circular pipe using ANSYS Fluent, this tutorial will show how to solve problems involving:
1. Turbulent fluid flow through a circular pipe and;
2. Heat transfer in a circular pipe.
In order to solve the momentum equation for the turbulent flow, we will use the k − ε closure model, which is covered in Lecture Notes 8 (LN08) on turbulence and turbulence modelling. Heat transfer is modelled by including the energy equation in the Fluent simulations.
Learning outcomes of this tutorial includes four key elements of CFD modelling, simulation and analysis:
1. Set up model Geometry;
2. Mesh generation;
3. Perform the Fluent analysis and;
4. Post-processing of your results (in CFD-Post).
2 Turbulent Pipe Flow
2.1 Problem description
We will consider a fluid flowing through the circular pipe of diameter D = 0.12936m and length L = 5m as sketched above. The inlet velocity is constant over the cross-sectional area and equal to U = 3.8755m/s. The working fluid has density ρ = 1.1644kg/m3 and fluid viscosity µ = 1.8487 × 10−5 kg/ms. These properties result in a Reynolds number of Re = ρUD/µ = 31577, which means the flow is turbulent. The outlet pressure po is equal to 1atmosphere (standard pressure). The conditions given above are from a classic experiment obtained in the well-known Princeton SuperPipe research facility at Princeton University, more details of this facility can be found here. We will use experimental data from this facility to validate our numerical results.
In the following will we will go through the workflow in ANSYS Fluent, consisting of the pre- processing, running the simulation and Post-processing. We assume that you have gained initial knowledge in setting up a model from the first Fluent practical and therefore this guide will not be as detailed as the previous Practical.
2.2 Pre-processing
Start a new project in the ANSYS workbench and save the project under a new *.wbpj file on your Homedrive.
Geometry
Create the geometry in DesignModeler or SpaceClaim with the dimensions given in the Problem Description. Note that, similar to the laminar pipe flow of practical 1, we can use an axisym- metric model domain. If you use DesingModeler Don’t forget to change the surface body to Fluid.
Mesh
We will start off by making a mesh of 450 divisions in the horizontal and 40 divisions in the vertical. Compared to the mesh of the laminar model in the previous practical we are not only increasing the number of divisions at our inlet and outlet boundaries, we also want a finer mesh (smaller cells) near the pipe wall, because for turbulent flows the velocity gradient (∂u/∂r) near the wall is much large compared to laminar flow. Closer towards the centreline the velocity gradient decreases so we can use a coarser mesh (larger cells). In order to to control the mesh size along an edge we will apply a mesh bias.
First of all, we have to create a structured mesh by applying Face Meshing to the mesh. Now, we still need to specify the number of cells on the boundaries of our domain. Instead of apply- ing the edge sizing to the inlet and outlet simultaneously, we will need to apply it to each edge independently because of the biasing we are going to apply.
Click Mesh Control in the mesh toolbar and select Sizing. Select the edge tool from the selection toolbar and click on the pipe inlet in the geometry, and then click on Apply in the Details View area. Set type to Number of Divisions and change the Number of Divisions to 40, set behaviour to Hard and, under Bias Type, choose the second option. This decreases the cell size in the direction of the wall (note that the bias direction depends on how your edge was sketched, ie, either from top-to-bottom or from bottom-to-top, so you might need to select the first option depending on how your geometry was created).
Finally, set the Bias Factor to 15. The bias factor indicates that the smallest division is 15times smaller than the largest division. Repeat the exercise for the outlet boundary conditions but this time apply the first option under Bias type. As an example, the Details View and the outlet should now look as follows:
If your resulting mesh is finer near the centreline than you may need to change the Bias Type. Now also apply the 450 regular divisions (ie, no biasing) to the pipewall and centreline.
Create named boundaries for the inlet, outlet, centreline, pipewall and the flow domain.
Model Setup
- General: set the solver to a steady, axisymmetric model, leave other settings at their default settings. You can investigate the mesh pressing Check, the output of the mesh check is displayed in the Console, which should display Done (so no error messages). If you select Report quality the console will give quality details, here the orthogonal quality is reported (value between 0 and 1, the higher is better) and the aspect ratio is the length/height of a cell. In this case the largest aspect ratio is about 36, which occurs for the cells closest to the pipewall.
- Models: set the Viscous model to the k-epsilon (2-eqn) model (k − ε turbulence model), and choose the Standard model. The near wall area is the area affected by viscosity and covers the region y+ < 30 where y+ is the non-dimensional distance from the wall. Recall from boundary layer flow (LN08) that y+ < 5 is the viscous sublayer while y+ > 30 is the log-layer, and the region inbetween is called the buffer layer. The k − ε model was not designed for the Near-Wall area (only for y+ > 30 where the flow is fully turbulent) and therefore we apply an additional “model” for the near wall region, the various model options are listed under Near-Wall treatment.
Standard wall functions are semi-empirical expressions that cover near-wall areas and essentially impose boundary conditions at some distance away from the wall in the “log- arithmic layer” (30 < y+ < 100), so that turbulence-model equations do not need to be solved very close to the wall. When you use the standard wall function the y+ value of the first grid point therefore (ideally) needs to fall in the range 30 < y+ < 100. Now, do not take these rules too strictly, for complex geometries and flow problems it will be very difficult to ensure that every grid point in your domain satisfies these y+ conditions! In practice, what you need to ensure as a user is that the majority of your cells satisfy the y+ conditions of your selected near-wall model - further verifcation/validation tests can demonstrate that the errors associated with the near-wall modelling are small and that your results are sensible.
Another approach is the two-layer approach which applies a separate model to the near- wall area, which allows the flow field to be resolved for grid points in the near-wall area, which then smoothly “blends” into the turbulence model for y+ > 30. The requirement for this approach is that the first point of the grid lies in the viscous sublayer region y+ < 5.
Here, we will select Enhanced Wall Treatment which applies either the standard wall function or the two-layer approach depending on the y+ of the first grid cell, therefore we don’t need to worry about the y+ value. The requirement is that the first grid cell is either in the viscous sublayer y+ < 5 or in the logarithmic layer y+ > 30, we will check later if we have satisfied this condition. Note that the Enhanced Wall Treatment can only be applied to hydraulically smooth walls, for hydraulically rough walls the flow has no viscous sublayer anymore therefore only the standard wall function can be applied. When applying standard wall functions the y+ of the first grid cell needs to be greater than 30.
Leave the model constants and other settings at their default values and press OK. En- sure that all the other models (Multiphase, Energy, etc.) in the Models task page are switched off.
- Materials: create a new fluid with the density and viscosity as outlined in the problem description and apply the new fluid to your fluid domain under cell zone conditions.
- Boundary conditions: We now need to specify the boundary conditions at the four borders of the domain. Note that the mathematical model now consists of two additional equations (for k and ε), therefore we also have to set up boundary conditions for these equations. Select Boundary conditions in the Navigation Pane and then:
1. centreline: set the centreline boundary Type to axis.
2. inlet: select inlet in the task page and click Edit, which will display the velocity inlet dialog box. Set the Velocity Specification Method to Components and set the Axial Velocity to 3.8755m/s. In addition we now also need to prescribe boundary conditions for the turbulence model. Therefore select Intensity and Hydraulic Di- ameter under the Specification Method and set a turbulence intensity of 2% and set the hydraulic Diameter equal to the pipe diameter, i.e. 0.12936m, see below. click OK to close the velocity inlet dialog box
3. outlet: Select outlet in the boundary conditions Task Page and check if the Type is set to pressure-outlet (this should be the case, but if it’s not, change the type to pressure-outlet). Click Edit and verify that the gauge pressure by default is set to 0Pa. You will also see Turbulence boundary conditions, which specify the turbulence present in case there is any backflow entering the domain through the outlet bound- ary. This will not be the case in the model so you can leave the Turbulence boundary conditions at their default values.
4. pipewall : select the pipewall in the boundary conditions Task Page and verify that the Type is set to wall. Verify that the settings are stationary wall with no slip shear condition.
- Reference values: set the reference density, length (diameter in this case), velocity and viscosity. The other values variables will not be used as part of this exercise.
2.3 Numerical solution
- Solution Methods: set the discretization for Momentum, Turbulent Kinetic Energy and Turbulent Dissipation Rate to Second Order upwind. Set the other settings as indicated below.
- Monitors: set the absolute criteria for all 5 equations to 10 −6 .
- If you are using 2021 R1 or above: Select Preferences from the File tab. In the window that appears, under General, in the box beside Default Format for I/O, select Legacy. This will allow Data File Quantities to be selected from File (after initialization), and for CFD-Post to load any additional quantities selected.
- Initialization: We now need to initialise all the variables in each grid cell. Select Solution Initialization from the Navigation Panel, select Standard Initialization and choose Inlet in the compute from drop-down menu. The values from for the turbulent kinetic energy and turbulent dissipation rate are automatically computed from the inlet turbulence boundary conditions. Press Initialize.
- Before we start the calculation, we want to save the skin friction coefficient again like we did in the first practical, but also the y+ value of the first cell above the wall. Go to File in the menu bar, select Data File Quantities and in the right column select Skin Friction Coefficient and Wall Yplus (i.e. y+ of first cell above the wall) and click OK.
- Now select Run Calculation from the Project Tree. Set the Number of Iterations to 1000, then press Calculate and wait till the solution has converged. Save your project and close the Fluent window.
2.4 Post-processing
You should now be able to plot in CFD-Post velocity vectors, velocity contours, velocity pro- files along the in- and outlet, and the skin friction coefficient along the pipe wall similar to the previous tutorial on laminar pipeflow. To improve the resolution of your line plots you can increase the Sample number under the Line Type setting of the location lines that you make (e.g. pipe inlet, pipe outlet etc).
In the set-up of the model k − ε model we applied the enhanced wall treatment approach to model the near wall flow. This approach was only valid when the first cell away from the wall was either in the region y+ < 5 or y+ > 30. Here we are going to plot y+ of the first grid cell to ensure this condition is satisfied. FLUENT conveniently stores they+ value of the first grid cell as the variable Yplus. Select the icon, give your chart a suitable name (which will be used as the title of the chart) and press OK. In the Details View, select pipe wall under location, in the X-axis tab select X under variable and in the Y-axis tab select Yplus, then click Apply. Your plot should look as follows:
Since y+ never exceeds 5 in the domain, the near-wall resolution is appropriate. Note that y+ is largest at the pipe inlet because the velocity gradient (∂ux/∂r), and hence the wallshear-stress, is largest here. See also p. 6-7 of LN08.
2.5 Tasks (Model Verification and Validation)
1. Change the mesh size to 60 × 800 and verify that your solution is 1) mesh-independent (if themeshis not converged refine it further), 2) converged and 3) that mass is conserved.
For assessing mesh-independence of your solution, you can use two methodologies; 1) you can apply a visual inspection (qualitative assessment) of, for example, a ux(y) ve- locity profile profile at the outlet of several simulations with continuous mesh refinement (similar to what you have done in the previous practical), or; 2) from two continuous mesh refinements, i and i + 1, we can quantitatively assess the difference between the velocity profiles ux(i)(y) and ux(i)+1 (y) (both are arrays of dimensions N), as follows:
(1)
A solution is assumed to be mesh-independent if ε is smaller than a prescribed value. There is unfortunately no universally accepted value for ε, it depends highly on the complexity of the fluid problem, which variable is of main interest, the required precision of your solution and experience. Here we are going to set ε ≤ 10−3, i.e. a mean absolute difference of less than 0.1%. Note that data from your plots in CFD-Post can be exported to comma- separated files (*.csv) by clicking in the Details of your plots.
2. Plot:
a) Contour plot of the axial velocity;
b) Profiles of the axial velocity at x = 0m, 1m, 3m and 5m.
Make decent looking graphs by increasing the default fontsize, line width and axes labels in the figures (found under the graph Details). For a clear looking contour plot you change the aspect ratio of the plot by apply a different Y-scale under the View tab. Export figures using the save picture 面 button. Give an interpretation of the results in the graphs.
3. Compare the skin friction coefficient at the outlet to the value you would obtain from the Moody diagram. Recall that the skin friction coefficient (Fanning friction factor) obtained from Fluent is equal to one fourth of the (Darcy-Weisbach) friction factor shown in most Moody diagrams.
4. There is no analytical solution to describe the turbulent velocity profile, but in order to validate the numerical solution, we can compare it with the SuperPipe experimental data for fully developed pipe flow. The experimental data can be found on MyAberdeen in an Excel file (containing only the data for half of the pipeline similar to our numerical results) or here. Compare your outlet velocity profile to the experimental velocity profile. Note that in the experimental data the vertical coordinate is given in non-dimensional form, and that r is the radial coordinate with r = 0 at the pipe centreline and r = R at the pipe wall (R is the radius of the pipe), and that the y-coordinate is the wall-normal coordinate with y = 0 at the pipe wall and y = R at the pipe centreline. Note that in CFD-Post the variable Y is the radial coordinate.
3 Pipe Flow with Heat Transfer
3.1 Problem description
In this problem we will heat a section of the pipewall to investigate the effect of the temperature on the fluid. Direct localised heating of pipelines is often used to keep the fluid above a prescribed temperature, such as the freeze point. The pipeline has a diameter D = 0.06m and length of L = 7m. The velocity at the inlet is constant over the cross-section with a value of U = 10m/s. The inlet (room) temperature of the fluid is 294.15K. The outlet pressure po is equal to 1atmosphere (standard pressure).
The wall is heated between x = 2m and x = 4m with a constant heat flux of 5kW/m2 . As- sume that the pipe wall is smooth. The heat capacity of the fluid is cp = 1006.43J/(kg K), the thermal conductivity is kp = 0.0242W/(m K), and molecular weight MW = 28.966g/mol, the fluid viscosity will be specified later.
Since much of this set-up is similar to the previous turbulent pipe flow exercise, except for the heat transfer part, only the main differences will be highlighted in the following description of the set-up.
3.2 Pre-processing
Geometry
To save some effort you can copy the geometry from the previous exercise (right-click the geometry and select Duplicate) and adjust the dimensions to match the dimensions of the new pipeline. Now we need to modify the sketch to indicate the heated wall section in the middle of the pipe. Click the Modify tab and select Split, select two points at the top of the rectangle and two at the bottom of the rectangle. Note that their location does not matter at the moment since we will dimensionalise them later.
Before we set the dimension we will first constrain the lower edges of the rectangle with the top edges of the rectangle. Click Constraints tab, and select Equal Length. Click the appropriate top and bottom edge sections and set them to be of equal length as shown below. Do the same for the middle sections.
Now set the two Horizontal dimensions of the heated section and click Generate to generate the surface.
Mesh
Drag a new mesh component to your Workbench and connect the Geometry component to the Mesh to start a new mesh. When you’ve opened up the Mesher the first thing to do is to create the Named Selections as shown below. Because the top and bottom borders now consist of multiple sections you have to select two sections for the pipewall (keep Ctrl pressed while clicking the line sections with your mouse), and three section for the centreline.
Create a mesh with 350 equally-spaced divisions applied to the top and bottom boundaries and 10 equally-spaced divisions on the inlet and outlet boundaries. In order to apply an equally- spaced division to the top boundary, we will select the three edges and then change the Type to Element Size and prescribe an element size of 0.02m as indicated below, which will give us 350 divisions over our 7m long pipeline. Repeat this exercise to apply 350divisions to the bottom boundary.
Model Setup
- General: set the solver to a steady, axisymmetric model, leave the other options at their default settings.
- Models: Because heat transfer is now part of our problem we need to switch on the energy model. Select Models, then select the Energy model, click edit, tick the selection box and click OK. The choice of the appropriate Viscous model will be part of the exercises below.
- Materials: Since the variations in absolute pressure are small in the pipeline, we will use a form. of the ideal gas law in which the pressure remains constant and therefore the density will only be affected by changes in the temperature (i.e., ρ = ρ(T)). Fluent calls this the “Incompressible ideal gas” law. By keeping the pressure constant we save some precious computational time. Under Density select incompressible-ideal-gas. Change any of the other properties to those given in the problem description, the viscosity will determined in the first exercise below. Click Change/Create and then Close the Create/Edit Materials windows. Apply the new fluid to your fluid domain under cell zone conditions.
- Boundary conditions: We now need to specify our boundary conditions at the four boundaries of our domain. Specifying the type of boundaries and prescribing the Momen- tum boundary conditions (velocity, pressure, turbulence - if applicable) is similar to what we did in Section 2 and therefore will not be repeated here. However, since we are also solving the energy equation, we have additional thermal boundary conditions. For the inlet, wall and outlet boundaries these can be specified in the Thermal tab:
1. inlet: Specify the temperature of the inlet fluid to 294.15K.
2. outlet : Here you need to specify again a backflow boundary condition in case flow enters the domain, you can set it to the inlet temperature 294.15K.
3. heated section: Select Heat Flux and set the value to that given in the problem description. Leave all other values at their default values.
4. pipewall : There is not heat transfer across the pipe wall (insulated wall), so set the Heat flux (W/m2 ) to 0.
3.2.1 Numerical solution
Before you proceed to the tasks in the next section we have given some reminders on the
numerical solution settings below. At this stage you should have sufficient information to
proceed - any missing information will be clarified in the tasks below.
- Solution Methods: set the discretization for the Momentum, Energy, Turbulent Kinetic Energy (if applicable), Turbulent Dissipation Rate (if applicable) to Second Order upwind.
- Monitors: set the absolute criteria for all equations (initially) to 10 −6 .
- Initialization: Select Solution Initialization from the Navigation Panel, select Standard Initialization and choose inlet in the Compute from drop-down menu.
- Similar to the previous exercise you want to save the skin friction coefficient and the Y+ values (recall that in Fluent the y+ value of the first cell from the wall is our main interest). Go to File in the menu bar, select Data File Quantities and in the right column select Skin Friction Coefficient and (if applicable) Wall Yplus and click OK (Note that after initialization you can also right-click on Run simulation and select the Data File Quantities).
3.3 Tasks
1. Assuming the density is 1.2kg/m3 at the inlet temperature, specify a value for the vis- cosity of the fluid so that the Re = 400. Perform a simulation with the specified mesh. Investigate the mesh sensitivity (keep horizontal divisions the same). Report the number of iterations needed to reach convergence for each simulation. Which mesh size will you use for any remaining simulations?
2. Based on your chosen mesh make the following figures:
a) Contour plot of the ux velocity;
b) Contour plot of the temperature;
3. Perform a simulation for Re = 20,000, and:
a) Report the main changes you have made to the Fluent set-up, investigate mesh sensitivity and convergence.
b) Make contour plots of the ux velocity and T;
c) Compare the velocity and temperature results to those obtained under item 2. Dis- cuss your results.
4. The bulk (or mean) temperature across the pipe cross-sectional area, Tm, characterises the average thermal energy state of the fluid. It is characterised as follows:
(2)
where m˙ is the mass flow rate, Ac is the pipe cross-sectional area and R is the pipe radius. Note that cp is taken outside the integral because (in our simulation) it is a constant. The temperature, T, velocity, u, and density, ρ (because ρ depends on T), are all a function of the pipe radius. The mass flow rate can be obtained from:
(3)
We can therefore obtain the following expression for Tm:
(4)
This expression can be solved in CFD-Post at any location along the pipeline. To do so select the expressions tab , right-click anywhere in the expressions window and select New. Give your expression a suitable name, e.g. “Tm0m” if you want to calculate the mean temperature at the inlet. In the Details window we can now enter Equation 4 above. The integral can be entered by right clicking in the empty area and select Functions→ CFD-Post→ lengthInt. By right-clicking within the integral brackets enter the required variables and specify the Location where you want to evaluate the in- tegral after the @ symbol, i.e. here you need to specify the line which you made earlier. Complete the expression to obtain the r.h.s. of Equation 4, note that π can be entered as “pi” and recall that the radial coordinate in CFD-Post is “Y” . When the expression is complete press , the value will appear in the “Value” box.
Evaluate the bulk temperature at every meter along the pipeline and make a plot of Tm against x and comment on your results.