Simcenter STAR-CCM+ Estimation of Spacial and Temporal Discretization Errors Using Richardson Extrapolation and Simcenter STAR-CCM+

Author: Rahil Doshi

Motivation

Computational Field Dynamics (CFD) is a powerful tool for simulating complex fluid phenomena, such as multiphase flows, turbulence, heat transfer, and chemical reactions. However, CFD simulations are subject to various sources of errors and uncertainties, such as modelling assumptions, numerical discretization, boundary conditions, and initial conditions. Therefore, it is essential to verify and validate the CFD simulations to ensure they are reliable and accurate. Verification is the process of checking that the numerical solution is free of errors and converges to the exact solution of the mathematical model. Validation is comparing the numerical solution with experimental or analytical data to assess how well the mathematical model represents the physical reality.

Scope of the tutorial

In this tutorial, you will learn how to estimate spatial and temporal discretization errors using Richardson extrapolation. Richardson extrapolation estimates the exact solution and the discretization error from a series of numerical solutions obtained with different mesh sizes or time steps.

• You will apply this method to gravity-driven flow, a multiphase flow problem where water flows from one chamber to another through a connecting channel.
• You will use Simcenter STAR-CCM+ software to set up and run the simulations and analyze the results.
• You will also learn how to choose a sensitive quantity, which is a representative integral physical quantity of the case, how to create a geometry part and a derived part for the velocity section, how to calculate the spatial and temporal discretization errors and the grid and time convergence indices, and how to compare the errors with the acceptable limits.
• By the end of this tutorial, you will be able to assess and quantify the spatial and temporal errors using Richardson extrapolation and determine if your solution is mesh and time step independent.

Geometry of the validation case
 

Case description

Initial conditions: The chamber on the left is filled with water, and the one on the right and the connecting channel with air.

Boundary conditions: All boundaries are solid walls except for the horizontal top left surface, where a constant (atmospheric) static pressure is applied.

Material Properties:

1. Air: Ideal Gas
2. Water: Constant density liquid

The fluid flows through the connecting channel through the left chamber to the right chamber. Note that a transient case is required to quantify the temporal discretization error. Due to the boundary conditions, this case will eventually reach a quasi-steady state where temporal discretization won't play a role. However, in this example, the case will only be run up to a physical time where the solution is still time-dependent.

Results

For a mesh base size of 0.625cm and a time step of 1e-4s,

Spatial Discretization Error: 0.707% and Temporal Discretization Error: 0.029%
 

Verification and Validation Procedure:

1. To begin with, perform the VOF: Gravity-Driven Flow tutorial (freeSurface.sim) without making any changes to the original case, using the provided volume mesh file (gridfs.ccm) to understand the model for multiphase flow problems, its material properties, initial conditions, reference values and the required field functions for obtaining the results.
2. To perform the estimation of spatial and temporal discretization errors, it is necessary to choose an integral quantity sensitive to mesh and time-step change that can be studied to analyze the results. In this case, the maximum velocity is chosen. The purpose of validating the case is to quantify the discretization error (spatial and temporal) and check if the error is within the acceptable limits. The spatial and temporal discretization error can be quantified by performing a mesh and time study to approach a mesh and time-step independent solution. To start with a mesh study, recreate the geometry of the problem (freeSurface1.sim) in the 3D CAD modeller to have a geometry available for creating a mesh.
   • To recreate the geometry, open a new file in Simcenter STAR-CCM+ (File > New > Ok). Under Geometry > 3D-CAD Models, right-click and select New. 
   • Right click on Features > XY and select Create sketch to create a new sketch on the XY plane.
   • Using the create line feature under Create sketch entities, sketch the approximate geometry, assign horizontal and vertical constraints wherever necessary and later specify the dimensions below. 
   • Click OK when done.  
   • Right-click on Sketch 1 and choose 'extrude'. Under distance, enter 1.0 cm as the extrude distance and click OK.
   • To create a new part, right-click Geometry > 3D-CAD Models > 3D-CAD Model 1 and select the new geometry part. Click OK to confirm. A new part is now created under the Parts section.
   • Additionally, create a velocity section (derived part) to check the velocity profile at the channel between the two chambers. Refer to the documentation and the sim file to set up the case.
3. For the first mesh, use a base size of 10cm (freeSurface1_m10.sim). Reduce the implicit unsteady time-step from 0.005s to 0.001s and increase the maximum inner iterations from 5 to 50 to make sure that the residuals tend to converge. In the first run a maximum velocity of 6.774m/s is obtained.
4. For the second run, reduce the mesh base size by half, 5cm (freeSurface1_m5.sim), keeping all the other parameters the same. For this case, the maximum velocity is 6.819m/s.
5. With the mesh base size of 2.5cm (freeSurface1_m2.5.sim) for the third run, the maximum velocity is now 7.019m/s. Looking back at the three values for maximum velocity, we observe a significant jump between runs 2 and 3, which indicates that the mesh base size of 2.5 needs to be sufficiently bigger and a finer mesh is required.
6. Performed another run with the same parameters as above but with a mesh base of 1.25cm (freeSurface1_m1.25.sim). There is a jump again in the maximum velocity to 7.272m/s. The CFL number also crosses 0.5 with a maximum of 0.736, which is unacceptable.
7. To limit the CFL number, the time step is further reduced to half (0.0005s) (freeSurface1_m1.25_ts0.0005.sim), and the maximum velocity obtained is 7.403m/s with a CFL number of 0.368.
8. The same time step of 0.0005s is also applied for the mesh base size of 2.5cm (freeSurface1_m2.5_ts0.0005.sim) and 5cm (freeSurface1_m5_ts0.0005.sim). Again, a big jump in the maximum velocities is required, indicating an even finer mesh base size.
9. The methodology of the Richardson Extrapolation requires, strictly speaking, a systematic mesh refinement, which means a refinement about a constant factor in the entire domain. In this example, however, the velocity differs between the differently sized meshes only significantly in a small portion of the domain (see proof below). If refinement in the rest of the domain doesn't change the results (not only because of the cell size but also because of the jump in cell size in case of local refinement), it won't affect the value of the representative scalar physical quantity used as the basis of the Richardson Extrapolation as well. Hence, where not needed, the "missing" refinement will be "invisible" in the computation of the discretization errors. Whether or not a non-systematic refinement will still give a reasonable estimate of the extrapolated result is strongly case-dependent and should be verified with a run using the "fully" systematically refined mesh in case of doubt.
   In the present case, we first need to verify that the velocity field "away" from the location where the maximum velocity is supposed to occur is not affected significantly by the mesh size. That can be done by firstly exporting the velocity field of all the above runs in steps 8 and 9 to tables (Velocity_1.25.csv, Velocity_2.5.csv, Velocity_5.csv) and secondly importing all the tables into the same simulation for comparison. This will help to refine the mesh in the relevant volume and keep the coarse mesh in the volume with little changes to reduce computational time. (Check freeSurface1_m5_0.625_ATS.sim for details from step 9 till step 14.)
   • To analyze the velocity data, import all tables from the above three .csv files in a single simulation under Tools > Tables. To do this, right-click Tables > New Table > File Table to select the .csv file from the dialogue box. 
   • 
   • Next, create three separate field functions, one for each .csv file, to interpolate the data from the table. Use the "interpolatePositionTable" function to define the field function and choose a vector under the value type of the field function. 
   • As seen from the above snippet, three separate vector field functions are created, i.e. Velocity_1.25cm, Velocity_2.5cm and Velocity_5cm. These field functions are used to create two other field functions, as shown in the pictures (User Field Function 1 and User Field Function 2).
   • Using these two scalar field functions, compute the magnitude of the velocity difference between the velocity data of Velocity_5cm and Velocity_2.5cm and Velocity_2.5cm and Velocity_1.25cm, respectively. 
   • The difference in velocities can now be visualized with a scalar scene. 
10. As observed, the majority of the velocity change occurs inside the channel. Hence, refining it outside this volume will have a negligible influence on the velocity magnitude used for the Richardson Extrapolation. To get a refined mesh only in the channel, a Volumetric Control can be used. To create that, right-click Geometry > Parts and choose New Shape Part > Cylinder. Enter the following parameter values and click Apply when done. 
11. To create a refined mesh in the cylindrical section, right-click on Geometry > Operations > Automated Mesh > Custom Controls and select New > Volumetric Control. Under Volumetric Control Parts, choose the Cylinder. Under Volumetric Control > Values > Custom Size, select the following properties for a base size of 0.00625 m. Note that a base size of 5cm is used for all subsequent cases outside of the refined volume. Right-click on Automated Mesh and select Execute to update the mesh.
12. The Adaptive Time-Step model can adjust the time-step size suitable for a given mesh. Note that a constant time step size is needed later for the runs used to extrapolate the time step independent solution. The purpose of using the Adaptive Time-Step model is only to conveniently determine the required time-step size for the subsequent run with a constant one. Right-click Models and click Edit. Add the Adaptive Time-Step model to the simulation. Under Adaptive Time-Step > Time-Step Providers, right-click to select Free Surface CFL Condition. Under properties, set the Max Condition limit to 0.4 as the criteria. 
13. To plot the time step size over the simulation time, first create a new report (named "dt" in the sim files) by right-clicking Report > New report > User > Minimum. Choose Time Step as the field function and part as Fluid, Pressure outlet. Then, right-click on the report and choose Create Monitor and Plot from Report. A plot is created under the Plots section. Run the simulation. The minimum value of the time-step for this mesh is 0.000322s.
14. Also, simulate an even finer mesh with a custom size of 0.003125m. (freeSurface1_m5_0.3125_ATS.sim). The smallest time-step for this mesh is 0.00016s. Since the maximum velocity value jumped from 7.609 m/s to 7.678 m/s from a minimum cell size of 0.00625m to 0.003125m, a minimum cell size of 0.003125m is chosen for the following series of runs.
15.  Simulate with a time-step of 1e-4s (freeSurface1_m5_1.25.sim). Remember to remove the Adaptive time step model now, as a constant time step size is used. The maximum velocity of 7.48m/s occurs at 0.3314s with a smaller peak later of 7.405m/s at 0.333s.
16. Similarly, run the simulation with a fine mesh size of 0.00625m and 0.003125m for the same time step (freeSurface1_m5_0.625.sim and freeSurface1_m5_0.3125.sim, respectively).
17. Surprisingly, the maximum velocity plots for the later runs show an unexpected peak in the maximum velocity magnitude. For the mesh size of 0.00625m, a peak of 9.255m/s at 0.3588s is seen. However, there is a local peak of 7.617m/s at 0.3306s, which is more consistent with the previous results. In the case of the mesh size of 0.003125m, the peak velocity shoots abruptly to 9.876m/s at 0.3556s. However, there is a similar local peak of 7.671m/s at 0.3298s. The Convective Courant Number has yet to cross 0.443. Ideally, the curve will follow a gradual slope without any additional peaks. A more detailed investigation is needed to investigate which cells are causing this sudden change in velocity.
18. First, create a scalar field function (named "Velocity Magnitude Liquid" in the freeSurface1_m5_0.3125_liquidcells.sim file) which neglects the velocity of those cells whose volume fraction of water is lower than 0.5 and assigns a velocity of those cells to 0m/s to those cells. This prevents the velocity values caused by air from influencing the results. Remember to change the field function to Velocity Magnitude Liquid under the Maximum Velocity report properties. This run shows a gradual increase in velocity, as expected, with a peak of 7.684m/s at 0.3311s. 
19. Similarly, run the other two simulations (freeSurface1_m5_0.625_liquidcells.sim and freeSurface1_m5_1.25_liquidcells.sim). The results for the mesh size of 0.0125m give the expected results, but this differs from the mesh size of 0.00625m. An unexpected hump is observed at 0.055s. The root cause must be identified, as this was never observed earlier. 
20. To check the hump in the maximum velocity plot, run the simulation again (freeSurface1_m5_0.625_liquidcells_hump_check.sim) but till only 0.055s (time till the maximum velocity at the hump region is seen). Here, we see the cell responsible for the peak in velocity. 
21. To analyze the root cause for the high velocity in this cell, start by creating a new vector scene to see the velocity vectors. It becomes evident that the velocity is perpendicular to the symmetric plane, which is unexpected and requires inspecting all quantities which influence velocity. Check the pressure field for any sudden pressure difference around the cell causing this behaviour. In this case, the pressure values in the neighbouring cells are comparable to that of the bad cell. Next, check for the pressure reconstruction gradient in that cell and its neighbours. By default, this field function is not registered as a field function. Enable Temporary Storage Retained in Solvers > Segregated Flow to have it registered. Now run the simulation for one more step to get updated values and select Pressure Recon by right-clicking on the colour bar in the vector scene. The pressure reconstruction gradient vector has a significant value with the direction pointing perpendicular to the symmetric plane. This is unexpected as this is a 2D case (velocity and pressure gradients are expected to be in the symmetry plane), and the fact that both the velocity and pressure reconstruction gradient at that cell point outwards the symmetry plane suggests that the gradient computation is associated with a significant error in that cell.   Velocity vector scene               Pressure reconstruction gradient
22. As a remedy, the Cell Quality Remediation model (refer to Remedying Cell Quality) can reduce the influence of potentially wrong gradients computed for cells with lousy quality. A scalar scene showing the Bad Cell indicator field function in the scalar scene helps identify any quality-related issues for that cell. The bad cell is orange in color, and its value for bad cell indicator is 2.0, indicating that the cell is suitable but is surrounded by bad neighbours (Cell Quality Remediation Field Functions Reference). This confirms the earlier assumption, and we now rerun the simulation with the cell quality remediation model (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun.sim). The Temporary Storage Retained option under Solvers > Segregated Flow can again be deactivated before starting the run. The hump in the velocity plot doesn't show up anymore, and the maximum velocity variation over time is now as expected.

Spatial Discretization Error using Richardson Extrapolation

1. Richardson's extrapolation (Wikipedia ) can be used to quantify the spatial and temporal discretisation error. The variation of a representative global scalar quantity is computed for systematically refined meshes/time step sizes, and a "mesh independent" value is derived. Comparing the extrapolated value with the computed ones for different mesh/time step sizes allows us to quantify the respective discretization error. Here, the global maximum of the velocity magnitude computed for mesh sizes of 0.003125m, 0.00625m and 0.0125m are used to quantify the spatial discretization error. The Richardson extrapolation is using approximations of the order O(h³)  or O(h⁶) using the value of the global quantity for any two distinct mesh/time step sizes (h and h/t) using the formulae                                                    N_j(h) = (t^j-1 * N_j-1(h/2) - N_j-1(h)) / (t^j-1 ​​​​-1) or                                                 N_j(h) = (t^2(j-1) * N_j-1(h/2) - N_j-1(h)) / (t^2(j-1) -1) respectively.
2. In this example, the mesh size (in the relevant spatial region) is refined systematically by reducing the cell edge length by a factor of 0.5 for the subsequent level. The spatial discretization error is approximated with O(h³). ​​​​​​
3. Table for generating Richardson's extrapolation

   	O(h)	O(h2)	O(h3)
   h	N1(h)=N(h)
   h/2	N1(h/2) =N(h/2)	N2(h)=2*N1(h/2)-N1(h)
   h/4	N1(h/4) =N(h/4)	N2(h/2) =2*N1(h/4)-N1(h/2)	N3(h)=(4*N2(h/2)-N2(h))/3

4. Results: 

   	O(h)	O(h2)	O(h3)
   h=1.25 cm	N1(h)=7.39912 m/s
   h/2=0.625 cm	N1(h/2) =7.65012 m/s	N2(h)=2*7.65012 -7.39912=7.90292 m/s
   h/4=0.3125 cm	N1(h/4) =7.68348 m/s	N2(h/2) =2*7.68348 -7.65012=7.71594 m/s	N3(h)= (4*7.71594-7.90292)/3=7.65361 m/s

5. Percentage error: (abs (measured value-extrapolated value)/extrapolated value) *100
   • for mesh size 1.25cm, error=3.325%
   • for mesh size 0.625cm, error=0.0338%
   • for mesh size 0.3125cm, error=0.390%
6. If cell sizes are in the correct order of magnitude, a monotonic decrease of the estimated spatial discretization error can be expected. That's not the case here. The Cell Quality Remediation model can be causing this as it affects how face values are computed from cell values. To check this hypothesis, max velocity magnitude values for different mesh sizes are used without the Cell Quality Remediation model (freeSurface1_m5_0.625_liquidcells.sim). In this case, the hump in the velocity variation is ignored, known to be caused by a bad cell.
7. Results: 

   	O(h)	O(h2)	O(h3)
   h=1.25 cm	N1(h)=7.39912 m/s
   h/2=0.625 cm	N1(h/2) =7.63376 m/s	N2(h)=2*7.63376 -7.39912=7.8684 m/s
   h/4=0.3125 cm	N1(h/4) =7.68348 m/s	N2(h/2) =2*7.68348 -7.63376=7.7332 m/s	N3(h)= (4*7.7332-7.8684)/3=7.68813 m/s

8. Percentage error: (abs (measured value-extrapolated value)/extrapolated value) *100
   • for mesh size 1.25cm, error=3.759%
   • for mesh size 0.625cm, error=0.707%
   • for mesh size 0.3125cm, error=0.060%
9. The errors are now decreasing monotonically, as expected.

Temporal Discretization Error using Richardson Extrapolation

1. The next step is to approximately quantify the temporal discretization error. The spatial discretization must remain unchanged for the runs with systematically refined time step sizes. Since the spatial discretization error of the mesh with a base size of 0.625cm is already as low as 0.707%, this mesh will be used for the subsequent time step size study.
2. Results are computed for time step sizes of 10e-5s (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun.sim), 5e-5s (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun_ts5e-5.sim) and 2.5e-5s (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun_ts2.5e-5.sim). The maximum time step size is chosen to meet the CFL constraint of the HRIC convection scheme (CFL number normal to the free surface < 0.5). The corresponding global maximum velocity magnitudes are 7.65102m/s, 7.64885m/s and 7.64856m/s, respectively. Note: We can use the simulation with the CQR model for the time step study as this model does not affect the temporal discretization error, only the spatial discretization error. 
3. The results for the Richardson's extrapolation are:

   	O(h)	O(h2)	O(h3)
   h=10e-5 s	N1(h)=7.65102 m/s
   h/2=5e-5 s	N1(h/2) =7.64885 m/s	N2(h)=2*7.64885-7.65102=7.6468 m/s
   h/4=2.5e-5 s	N1(h/4) =7.64856 m/s	N2(h/2) =2*7.64856-7.64885=7.64827 m/s	N3(h)= (4*7.64827-7.64668)/3=7.64876 m/s

4. Percentage error: (abs (measured value-extrapolated value)/extrapolated value) *100
   • for time step of 10e-5 s, error=29.0242e-3%
   • for time step of 5e-5 s, error=0.6537e-3%
   • for time step of 2.5e-5 s, error=3.1377e-3%
5. The CFL constraint of the HRIC convection scheme leads to a very small temporal discretization error, even for the largest time step size of the series.
6. To check what a reasonable time step size to discretize the temporal variation of the solution variables of this case is, multi-stepping for volume fraction transport can be used to eliminate the time step size constraint due to the HRIC convection scheme.
7. Three runs are conducted using implicit multi-stepping for VOF with the following settings: 
   • Solvers->Implicit Unsteady->Time-Step: 10E-5 s and Solvers->Segregated VOF->Implicit Multi-Step->Number of Steps: 1 (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun_multistepping10e-5.sim)
   • Solvers->Implicit Unsteady->Time-Step: 20E-5 s and Solvers->Segregated VOF->Implicit Multi-Step->Number of Steps: 2 (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun_multistepping20e-5.sim)
   • Solvers->Implicit Unsteady->Time-Step: 40E-5 s and Solvers->Segregated VOF->Implicit Multi-Step->Number of Steps: 4 (freeSurface1_m5_0.625_liquidcells_hump_checked_rerun_multistepping40e-5.sim)
8. The results for the Richardson's extrapolation are:

   	O(h)	O(h2)
   h=40E-5 s	N1(h)=7.59351 m/s
   h/2=20E-5 s	N1(h/2) =7.6571 m/s	N2(h)=2*7.6571-7.59351=7.72069 m/s
   h/4=10E-5 s	N1(h/4) =7.65102 m/s	N2(h/2) =2*7.65102-7.6571=7.64494 m/s

9. Percentage error: (abs (measured value-extrapolated value)/extrapolated value) *100
   • for time step of 40e-5 s, error=0.34358%
   • for time step of 20e-5 s, error=0.49096%
   • for time step of 10e-5 s, error=0.41117%
10. The larger global time step only adds a little to the error since it doesn't affect velocity magnitude significantly. The non-nonlinearity in time is already discretized appropriately, even for the largest time step size. The slight change observed is due to the noise during the simulation runs. All the error values are already very low and within acceptable limits.

 

Conclusion

This tutorial demonstrated how to use the Richardson Extrapolation to estimate the spatial and temporal discretization error in CFD simulations on a practical example, including how to still get value out of the exercise even though the methodology wasn't (couldn't be) applied rigorously. The method was applied to gravity-driven flow, a multiphase flow problem. Simcenter STAR-CCM+ software was used for setting up, running, and analyzing the simulations. A sensitive quantity was chosen to validate the case, and a geometry part and a derived part were created for the velocity section. The spatial and temporal discretization errors and the grid and time convergence indices were calculated and compared with the acceptable limits.

The spatial discretization error was 0.707%, and the temporal discretization error was 0.029%, within the recommended range of less than 1%. This indicated that the solution was mesh and time step independent and had a high degree of accuracy. 

This method can verify and validate other CFD problems involving complex fluid phenomena. By doing so, the reliability and accuracy of CFD simulations can be ensured.

Note: 

The method used in this tutorial is a simplified version of Richardson extrapolation, which assumes that the discretization error is proportional to the mesh size or the time step raised to some power. This may only be true for some cases, especially for highly nonlinear or unsteady problems. Therefore, it is recommended to check the validity of this assumption by plotting the discretization error versus the mesh size or the time step on a log-log scale and verifying that the slope is approximately constant. If the slope is not constant, then the method may not be applicable or accurate. More advanced verification and validation methods, such as generalized Richardson extrapolation, grid refinement studies, and uncertainty quantification, may be required. For more information on these methods, please refer to Guide: Guide for the Verification and Validation of Computational Fluid Dynamics Simulations.

All the simulation files in this tutorial are in the compressed folder attached to this article. The Simcenter STAR-CCM+ version used in this tutorial is 2022.1 (17.02).

See also: 

• Richardson extrapolation
• Guide: Guide for the Verification and Validation of Computational Fluid Dynamics Simulations
• What is the effect of changing the temporal discretization in unsteady simulations?
• Swiss Simcenter STAR-CCM+ Knife: [5] Pixel comparison: Integrate image processing tools in your CFD workflow.