WOLFRAM

Wolfram Language & System Documentation Center
Fluid Dynamics Model Verification Tests

Fluid Dynamics Model Verification Tests

Stationary TestsTest Result Inspection
Transient TestsReferences

This notebook contains tests that verify that the fluid dynamics partial differential equations (PDE) model works as expected. To run all tests, SelectAll and press Shift+Enter. The results will then be in the section Test Result Inspection.

Note that these tests can also serve as a basis for developing your own fluid dynamics models. As such, the tests are grouped into stationary (time-independent) and transient (time-dependent) tests. In both categories, both Cartesian and axisymmetric test models can be found.

In each test case, the visualization section is there to provide post-processing results for inspection; however, it is not a necessary part of the test. In the interest of saving runtime and reducing memory consumption, the cells in the visualization section are set to not be evaluatable. To make these cells evaluatable, select the cells in question, choose Cell ▶ Cell Properties and make sure "Evaluatable" is ticked.

The fluid dynamics equations are used to solve for the velocity field and pressure in a fluid flow model. Please refer to the information provided in "Laminar Flow" for more general theoretical background for fluid flow analysis, as this monograph assumes basic knowledge of fluid dynamics.

To avoid keeping memory-intensive previous results, set the history length to 0:
Load the Finite Element package:

Stationary Tests

This section contains examples of stationary, which means non-time-dependent, fluid dynamics PDE models for the validation.

2D

This section contains examples of a 2D stationary fluid dynamics PDE model.

FluidDynamics-FEM-Stationary-2D-FluidDynamics-0001

The following test demonstrates a 2D stationary flow around a cylinder. The simulation domain is a rectangular region with a width of and a height of and a cylinder of radius located at a point . Notice that the -coordinate position of cylinder center is slightly below the center of the height of the domain. This is done on purpose to prevent symmetry of the domain.

The flow inside the region is driven by prescribed parabolic inflow on the left side of the domain. The other boundary conditions are no-slip boundary condition on the top and bottom and no traction boundary condition on the right hand side.

Test Reference

[1] https://wwwold.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html (last accessed: May 22th, 2025)

Equations
Set up the model variables and parameters:
Set up the PDE:
Region
Set up the model region:
Create and visualize the mesh:
Boundary Conditions
Set up the boundary condition on the inflow:
Set up the no-slip boundary condition on the wall:
Set up the no traction boundary condition on the outflow:
Set up the pressure boundary condition to make the solution unique:
Wrap all Dirichlet boundary conditions into one variable:
Solution
Solve the PDE model and monitor time/memory usage:
Test 1

This test compares the pressure in front of the circle at point with the pressure at the back of the circle at point .

Compute the pressure difference and define the reference pressure difference:
Verify the pressure difference:
Test 2

This test computes the drag acting on the circle. The drag can be computed as the first component of the traction

where is the Cauchy stress tensor, is the outward boundary unit normal and the integral is over the boundary of the circle. The resulting number is then scaled to obtain dimensionless drag.

Compute the gradient of the velocity:
Compute the strain rate:
Compute the shear rate:
Compute the stress tensor:
Set up the boundary unit normal:
Set up the traction:

Note the minus in this expression. The represents an outward unit boundary normal vector. However, this verification test requires an outward unit boundary normal vector to the cylinder in the domain; therefore, the vector now points in the opposite direction and the minus sign is necessary.

Compute the mean velocity and define the characteristic length:
Compute the dimensionless drag and define the reference drag:
Verify the drag:
Test 3

This test computes the lift acting on the circle. The lift can be computed as the second component of the traction

where is the Cauchy stress tensor, is the outward boundary unit normal and the integral is over the boundary of the circle. The resulting number is then scaled to obtain dimensionless lift.

Define the reference lift:
Compute the dimensionless lift:
Verify the lift:
Visualization

The following cells are marked as not evaluatable to save the runtime and memory consumed. To make these cells evaluatable, select the cells in question, choose Cell ▶ Cell Properties and make sure "Evaluatable" is ticked.

Plot the velocity field using the VectorPlot:
Plot the velocity field using the StreamPlot:
Plot the pressure using ContourPlot:

FluidDynamics-FEM-Stationary-2D-FluidDynamics-0002

The following test demonstrates a 2D stationary flow in a lid driven cavity. The domain is a unit rectangle filled with fluid. On the top there is a mechanism that drives the fluid with a fixed flow velocity in the positive direction. The remaining sides are walls and have a no-slip condition.

This test considers a non-Newtonian power-law fluid with various values of parameters. The equations for the power-law fluid can be found in the Laminar Flow monograph. The power-law exponent and a dimensionless Reynolds-like number are used as parameters for the flow. These parameters are varied to obtain different behavior of the fluid. The Reynolds-like number is defined as

where is mass density, is characteristic length, is power-law viscosity, is reference velocity and is power-law exponent. In the case of this test, the characteristic length and reference velocity are both equal to . The power-law viscosity is also set to . This means that the value of the Reynolds-like number is equal to the value of the mass density.

Test Reference

[2] P. Neofytou. A 3rd order upwind finite volume method for generalised Newtonian fluid flows. Advances in Engineering Software, vol. 36, issue 10: pp. 664-680 (2005).

Equations
Set up the model variables and parameters:
Set up the PDE:
Region
Set up the domain:
Create and visualize the mesh:
Boundary Conditions
Set up the boundary condition on the top of the domain:
Set up the no-slip boundary condition on the rest of the domain:
Set up the boundary condition for the pressure to make the solution unique:
Wrap all boundary conditions into one variable:
Solution
Set up a parametric function which solves the equations for given parameters:

The ExtrapolationHandler is necessary due to the usage of NMaxValue in the tests.

Define a Reynolds number parameter lists and create a list of parameter pairs:
Solve the equations for given parameters:
Test 1

This test computes and verify the dimensionless maximum of the wall shear stress at the bottom part of the boundary. The wall shear stress is the tangent part of the traction; however, in this simple domain it is equal to the shear rate. Note that the wall shear stress is sensitive to mesh refinement and a fine mesh is used.

The precise definition of the wall shear for more general domains can be found in the Blood Flow Modeling in Cerebral Aneurysm application example.

Extract velocity components from the list of solutions:
Compute velocity gradients:
Compute strain rates:
Compute shear rates:
Find maximum wall shear stress at the bottom of the domain:
Compute the dimensionless wall shear stress maxima:
Define reference values for the dimensionless wall shear stress:
Verify the wall shear stress:
Visualization

The following cells are marked as not evaluatable to save the runtime and memory consumed. To make these cells evaluatable, select the cells in question, choose Cell ▶ Cell Properties and make sure "Evaluatable" is ticked.

Create StreamPlots for each flow and partition the resulting plots into groups of three:
Create an interactive plot which allows to change the Reynolds-like number and the power-law exponent:

2D Axisymmetric

This section contains examples of a 2D stationary axisymmetric fluid dynamics PDE model.

FluidDynamics-FEM-Stationary-2D-Axisymmetric-FluidDynamics-0001

The following test case considers an axisymmetric flow in a cylindrical pipe driven by a pressure gradient. The model domain is reduced from cylinder to rectangle thanks to the assumed symmetry. The resulting flow is then compared to the exact solution.

Since the exact solution for the problem is knows, all parameters in this test such can be chosen arbitrarily.

Test Reference

[3] M. O. Deville. An Introduction to the Mechanics of Incompressible Fluids, Chapter 3. Springer (2022).

Equations
Set up the model variables and parameters:
Set up the PDE:
Region
Set up the model region:
Boundary Conditions
Set up the boundary condition on the inflow:
Set up the symmetry boundary condition:
Set up the no-slip boundary condition on the wall:
Set up the do-nothing boundary condition on the outflow:
Set up the pressure boundary condition on the outlet:
Set up the boundary condition on the inlet:
Wrap all Dirichlet boundary conditions into one variable:
Solution
Solve the PDE model and monitor time/memory usage:

The "ExtrapolationHandler" is necessary due to the usage of NMaxValue in one of the tests. NMaxValue may evaluate the function outside of it's domain. To allow for that and to avoid possible messages the "ExtrapolationHandler" is specified.

Test 1

The derivative of the pressure with respect to must be equal to .

Compute the pressure derivative:
Verify the pressure derivative:
Test 2

The flow has to have parabolic profile all along the domain.

Set up and plot the profile function:
Verify the flow profile:
Test 3

The maximum velocity of the flow has to match the maximum velocity of the exact solution.

Compute the maximum velocity:
Verify the maximum velocity:
Test 4

This test verifies that the shear stress acting on the wall matches the shear stress of the exact solution.

Compute the exact shear stress:
Define a function which computes the shear stress:
Verify the shear stress:

Transient Tests

This section contains examples of transient, which means time-dependent, fluid dynamics PDE models for the validation.

2D Axisymmetric

This section contains examples of a transient axisymmetric fluid dynamics PDE model.

FluidDynamics-FEM-Transient-2D-Axisymmetric-FluidDynamics-0001

The following test case considers an axisymmetric flow in a cylindrical pipe driven by a pressure gradient. The model domain is reduced from cylinder to rectangle thanks to the assumed symmetry. The resulting flow is then compared to the exact solution.

Since the exact solution for the problem is knows, all parameters in this test such can be chosen arbitrarily.

Test Reference

[4] M. O. Deville. An Introduction to the Mechanics of Incompressible Fluids, Chapter 3. Springer (2022).

Equations
Set up the model variables and parameters:
Set up the PDE:
Region
Set up the model region:
Boundary Conditions
Set up the boundary condition on the inflow:
Set up the symmetry boundary condition:
Set up the no-slip boundary condition on the wall:
Set up the do-nothing boundary condition on the outflow:
Set up the pressure boundary condition on the outlet:
Set up the boundary condition on the inlet:
Wrap all Dirichlet boundary conditions into one variable:
Initial Condition
Set up the initial condition:
Solution
Solve the PDE model and monitor time/memory usage:
Test 1

This test verifies that the solution has a constant pressure gradient at the time .

Compute the derivative of the pressure in the direction:
Verify that the derivative of the pressure is equal to :
Test 2

This test compares the computed flow profile with the profile of the exact solution. Unlike in the stationary case, in the transient case the exact profile involves an infinite series. This means that an approximation using a truncated sum has to be used. This is most problematic for small times as the profile converges to the parabolic profile of the stationary case for time going to infinity. For example, the approximated exact profile does not satisfy the initial condition. As a result, this test is less accurate than other tests.

Set up the profile function which approximates the exact profile:
Plot the approximation of the exact flow profile:
Set up a function which computes the norm of the difference of the computed solution and exact profile:
Compute the norm at each time step:
Verify the flow profile:

3D

This section contains examples of a 3D transient fluid dynamics PDE model.

FluidDynamics-FEM-Transient-3D-FluidDynamics-0001

This test computes a flow inside a cube. An exact solution is used as the initial and boundary conditions. The resulting flow is then compared to the exact solution.

Since the exact solution for the problem is known, most parameters in this test can be chosen arbitrarily. The only exception is the Reynolds number. This is due to the fact that the exact solution derived in the reference article uses a different dimensionless form of the Navier-Stokes equations. This form is equivalent to one implemented in the Wolfram Language only for Reynolds number equal to 1.

Test Reference

[5] C. R. Ethier and D. A. Steinman. Exact fully 3D Navier-Stokes solutions for benchmarking. International Journal for Numerical Methods in Fluids, vol. 19: pp. 369-375 (1994).

Equations
Set up the model variables and parameters:
Set up the PDE:
Region
Set up the model region and create and visualize the mesh:
Exact Solution
Define the exact solution:
Verify that the defined functions solve the Navier-Stokes equations:

The solution defined in the previous cell contains two parameters which have to be specified for numerical computations.

Define parameters found in the exact solution:
Boundary Conditions

The boundary conditions are based on the exact solution.

Set up the boundary condition:
Initial Condition

The exact solution is also used as the initial condition.

Set up the initial condition:
Solution
Solve the PDE model and monitor time/memory usage:
Test 2

This test verifies that the solution is identical to the exact solution.

Set up a function which computes the norm of the difference of the computed solution and the exact solution:
Compute the norm of the difference at time steps:
Verify that the norm of the difference is small:
Visualization

The following cells are marked as not evaluatable to save the runtime and memory consumed. To make these cells evaluatable, select the cells in question, choose Cell ▶ Cell Properties and make sure "Evaluatable" is ticked.

Plot the velocity at each timestep using VectorPlot3D:
Animate the resulting plots:
Plot the velocity at each timestep using StreamPlot3D:
Animate the resulting plots:

Test Result Inspection

This section contains the evaluation of the test runs. It works by collecting all instances of TestResultObject and generating a TestReport.

Extract TestResultObject from the notebook and generate a TestReport:
Inspect the failed tests run:

If the preceding table is empty, all tests succeeded.

References

1.  https://wwwold.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html (last accessed: May 22th, 2025)

2.  P. Neofytou. A 3rd order upwind finite volume method for generalised Newtonian fluid flows. Advances in Engineering Software, vol. 36, issue 10: pp. 664-680 (2005).

3.  M. O. Deville. An Introduction to the Mechanics of Incompressible Fluids, Chapter 3. Springer (2022).

4.  C. R. Ethier and D. A. Steinman. Exact fully 3D Navier-Stokes solutions for benchmarking. International Journal for Numerical Methods in Fluids, vol. 19: pp. 369-375 (1994).