Validation
Validation of the LUMA solver, including its FSI capabilities has been performed for a number of cases and compared to available data in the literature. Available validation results are summarised below.
The flow past a circular cylinder at is simulated using LUMA v1.4.0. Drag and lift values over the cylinder are computed by LUMA with coefficients,
and
as well as the pressure coefficient over the cylinder surface calculated in post-processing. Results are compared to the set of experiments presented by Park and Kwon [1].
Consider a two dimensional cylinder of diameter located at the centre of a square
computational domain. Freestream velocity boundary conditions are used on all domain edges with a value set to
. The cylinder boundary uses the bounce-back boundary condition. The computational domain is shown in Figure 1.
Figure 1: Computational domain of a cylinder
The table below presents the primary setup configuration in non-dimensional units as used in the definitions file.
| Parameter | Value |
|---|---|
| L_DIMS | 2 |
| L_RESOLUTION | 8 |
| L_TIMESTEP | 0.0125 |
| L_BX | 50 |
| L_BY | 50 |
| L_RE | 2-160 |
| L_GEOMETRY_FILE | defined |
| L_INLET_ON | defined |
| L_FREESTREAM_TUNNEL | defined |
| L_NUM_LEVELS | 4 |
The validation of pressure coefficient around a cylinder is based on [1] who present pressure coefficient for . We also compare to the results of Dennis and Chang [2] for
and Norberg [3] at
. Surface pressure coefficient is calculated at cylinder surface between 0 and 180 degrees as
Figure 2a and 2b present the surface pressure coefficient computed with LUMA compared with Park [1]. Figure 2a considers while Figure 2b shows
. At low
it is possible to see good qualitative results, but as soon as
increase quantitative results improve considerably obtaining identical values to the reference.
Figure 2a: Pressure coefficient for
Figure 2b: Pressure coefficient for
The pressure coefficient at the base and stagnation
points. Results were compare with Park[1], Dennis [2], Williamson and Roshko [3], Norberg [4] and Henderson[5]. In all cases LUMA achieves excellent results as indicated in Figure 3.
Figure 3: Pressure coefficient at base and stagnation
points. * LUMA, . Park [1], + Dennis [2], Fonberg [6], x Williamson[3], ^ Norberg[4], v Henderson[5].
Figure 4 shows a comparison of drag coefficient for and Strouhal number. In both cases, good agreement with the literature is obtained.
Figure 4: Drag coefficient and Strouhal number
Figure 5b compares the length of time average separation bubble (depicted in Figure 5a). Again, the results from LUMA are in excellent agreement with the reference results.
Figure 5a: Separation bubble definition
Figure 5b: Length of time averaged separation bubble vs Re
[1] J. PARK, K. KWON, AND H. CHOI, Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160, KSME International Journal, 12 (1998), pp. 1200–1205.
[2] S. C. R. DENNIS AND G.-Z. CHANG, Numerical solutions for steady flow past a circular Cylinder at Reynolds numbers up to 100, Journal of Fluid Mechanics, 42 (2006), p. 471.
[3] WILLIAMSON AND ROSHKO, Vortex formation in the wake of an oscillating cylinder, Journal of Fluids and Structures, 2 (1988), pp. 355–381.
[4] C. NORBERG, An experimental investigation of the flow around a circular cylinder: influence of aspect ratio, Journal of Fluid Mechanics, 258 (1994), pp. 287–316.
[5] R. D. HENDERSON, Details of the drag curve near the onset of vortex shedding, Physics of Fluids, 7 (1995), pp. 2102–2104.
[6] B. FORNBERG, A numerical study of steady viscous flow past a circular cylinder, Journal of Fluid Mechanics, 98 (1980), p. 819
A similar investigation to the previous case can be conducted for an inclined flat plate in 2D. The problem presented by Tiara and Colonius [1] is simulated using LUMA v1.4.0.
We consider a computational domain with the same boundary conditions as the previous case. This time the freestream velocity is set to
. The flat plate is also treated as bounce back boundary condition and angles of incident of
. The computational domain is shown in Figure 1.
Figure 1: Computational domain of a flat plate.
The table below presents the primary setup configuration in non-dimensional units as used in the definitions file.
| Parameter | Value |
|---|---|
| L_DIMS | 2 |
| L_RESOLUTION | 8 |
| L_TIMESTEP | 0.0625 |
| L_BX | 20 |
| L_BY | 20 |
| L_RE | 300 |
| L_GEOMETRY_FILE | defined |
| L_INLET_ON | defined |
| L_FREESTREAM_TUNNEL | defined |
| L_NUM_LEVELS | 4 |
Figure 2 presents the drag and lift coefficient for the flat plate at different angles of incidence and results are compared with [1] showing excellent agreement.
Figure 2: Drag and Lift coefficient flat plate, Re=300.
[1] TAIRA, K., & COLONIUS, T. (2009). Three-dimensional flows around low-aspect-ratio flat-plate wings at low Reynolds numbers. Journal of Fluid Mechanics, 623, 187-207. doi:10.1017/S0022112008005314
This test case is a very popular benchmark for testing fluid-structure interaction solvers. The case, initially presented by Turek & Hron [1], consists of a rigid cylinder with an attached flexible plate embedded within a channel flow. A parabolic velocity profile is specified at the inlet which is ramped up over the first two seconds. No-slip conditions are used at the channel walls, and the outlet boundary condition is left up to the user. Figure 1 shows a schematic of the computational domain. Note that the cylinder is placed off-centre from the channel to allow the instabilities to develop naturally. The original benchmark provides data for three different Reynolds numbers (based on the diameter of the cylinder). Here we present tip displacements and forces for the intermediate Reynolds number (Re = 100) case, which exhibits steady periodic flapping of the plate, and compare against the benchmark data. As of v1.6.1 LUMA only supports fluid-structure interaction solutions in 2D, therefore only the 2D results are presented. The flow conditions for this case are given by:
- Mean Velocity = 1.0 m/s
- Density = 1000 kg/
- Kinematic viscosity = 0.001 m2/s
- Re = 100
Figure 1: Schematic of case with dimensions given in physical units.
To properly capture the profile of the tip dynamics and force histories a reasonably high sample rate for the output writers is required (approximately every 0.01s is enough). To allow the tip dynamics to reach a steady periodic state the user should specify enough time steps to cover at least 20 seconds of physical time. It has been found that even at relatively low grid resolutions (20 lattice sites across the diameter of the cylinder) LUMA gives reasonably good results. Either regularised or bounceback can be used at the walls to impose a no-slip velocity however care must be taken when setting the domain height so as to ensure the distance between the points where the no-slip condition is applied is 0.41m. The macro L_WALL_THICKNESS_BOTTOM is useful for this. Since the Reynolds number is set assuming a non-dimensional length scale of unity, to preserve the dimensions given in the original benchmark (where the Reynolds number is based on D=0.1m) the input Reynolds number must be scaled accordingly.
Since this case requires the use of an immersed body a geometry input file is required. This file should contain two lines – one describing the cylinder dimensions and one describing the plate dimensions and material properties. All flexible bodies must be of type IBM. Care must be taken not to have IBM markers very close to each other in the region where the plate is attached to the cylinder (the L_UNIVERSAL_EPSILON_CALC flag helps to mitigate against instabilities which may arise because of this). A relaxation parameter of approximately 0.5 works well for this intermediate case.
The table below presents the relevant options and parameters set in the definitions.h file. Note that LUMA accepts dimensionless quantities as inputs however the reference values have been carefully chosen here to preserve the dimensions of the original benchmark.
| Parameter | Value |
|---|---|
| L_OUT_EVERY | 20 |
| L_OUT_EVERY_FORCES | 20 |
| L_HDF5_OUTPUT | Defined |
| L_LD_OUT | Defined |
| L_NO_FLOW | Defined |
| L_TOTAL_TIMESTEPS | 60000 |
| L_DIMS | 2 |
| L_RESOLUTION | 200 |
| L_TIMESTEP | 0.0005 |
| L_BX | (2.5+0.5*(L_WALL_THICKNESS_LEFT+L_WALL_THICKNESS_RIGHT)) |
| L_BY | (0.41+(L_WALL_THICKNESS_BOTTOM+L_WALL_THICKNESS_TOP)) |
| L_PHYSICAL_U | 1.0 |
| L_PHYSICAL_RHO | 1000.0 |
| L_PARABOLIC_INLET | Defined |
| L_UX0 | 1.0 |
| L_UY0 | 0.0 |
| L_RHOIN | 1.0 |
| L_RE | 1000 |
| L_GEOMETRY_FILE | Defined |
| L_VTK_BODY_WRITE | Defined |
| L_IBM_ON | Defined |
| L_UNIVERSAL_EPSILON_CALC | Defined |
| L_NB_ALPHA | 0.25 |
| L_NB_DELTA | 0.5 |
| L_RELAX | 0.5 |
| L_WRITE_TIP_POSITIONS | Defined |
| L_WALL_LEFT | eVelocity |
| L_WALL_RIGHT | ePressure |
| L_WALL_BOTTOM | eSolid |
| L_WALL_TOP | eSolid |
| L_REGULARISED_BOUNDARIES | Defined |
| L_VELOCITY_RAMP | 2 |
| L_PRESSURE_DELTA | 0.0 |
| L_NUM_LEVELS | 0 |
The geometry file should be configured with the following information:
Line 1: FILAMENT_ARRAY IBM 0 0 1 0.2525 0.2 0.0 0.0 0.0 0.0 0.35 0.02 0.0 0.0 0.0 FLEXIBLE 10 CLAMPED 10000.0 1.4e6
Line 2: CIRCLE_SPHERE IBM 0 0 0.2025 0.2 0.0 0.05 RIGID
Figures 2 and 3 show the velocity and vorticity contours obtained with LUMA for the intermediate Reynolds number case. Figure 4 shows the tip dynamics compared against the original benchmark data. The results show very good agreement even for a relatively coarse grid.
Figure 2: Velocity counters for 2D FSI Benchmark.
Figure 3: Vorticity contours for 2D FSI benchmark.
Figure 4: Tip displacement histories for LUMA and benchmark data.
[1] S. Turek and J. Hron, “Proposal for Numerical Benchmarking of Fluid-Structure Interaction between an Elastic Object and Laminar Incompressible Flow” Fluid-Structure Interaction. Lecture Notes in Computational Science and Engineering, vol 53. Springer, Berlin, Heidelberg , 2006.
This test case is used to validate the solver performance by simulating the wind tunnel experiment by Meinders and Hanjalic [1]. The experiment consists of a grid of cubes, in the streamwise and spanwise direction respectively, in a turbulent channel flow. Figure 1 shows a schematic view of the central region of the experimental setup which can be used to configure the simulation. Meinder and Hanjalic provide reliable experimental data that has been used in previous simulation validation studies [2,3]. Vertical and horizontal velocity profiles together with vector fields and streak lines at different 2D slices around the cube are presented and used for comparison here. The representative dimensions of the experiment in physical units are:
- Channel height: h=0.051m
- Height of the cubes: H=0.015m
- Separation between cubes: S = 4H
- Bulk velocity:
m/s
- Bulk Re number:
Figure 1a: Top view of a section of the wind tunnel experiment by [1]. The dot-dash blue rectangle encloses the numerical domain used by [2] together with its boundary conditions.
Figure 1b: Side view of a section of the wind tunnel experiment by [1]. The dot-dash blue rectangle encloses the numerical domain used by [2] together with its boundary conditions.
This test case simulated the flow around the central cube of the cube array following the setup in [2] and using a Smagorinsky turbulence model configuration according to [4]. LUMA v1.5 was used to run this case. Figure 1 shows the size of the domain and the boundary conditions. The boundary condition applied to the walls is no-slip.
The flow is driven by a pressure gradient between the inlet and outlet of the domain. This pressure gradient is modelled as a body force and applied to all the cells of the domain using LUMA’s gravity force feature. As the value of the pressure gradient force is not given explicitly in the literature, this value must be iteratively determined through trial and error. A series of test runs are conducted until the volumetric flow rate in the simulation is close to the volume flow rate in the experiment
The table below presents the relevant LUMA parameters for this case in dimensionless units, the basic configuration units used in the LUMA definitions file. Using these parameters, the volume flow rate is approximately 8% lower than the experiment but as only qualitative agreement is sought in this test case this is deemed sufficient. The length, velocity and time scales used to convert between dimensionless and physical units depend on the choice of reference values. These are selected from the problem definition as:
- Reference length: Dimension of the cube H = 0.015m
- Reference velocity: Bulk velocity upstream of the cubes
- Reference time
s
| Parameter | Value |
|---|---|
| L_GRAVITY_FORCE | 0.00589 |
| L_GRAVITY_DIRECTION | eXDirection |
| L_SMAG | 0.3 |
| L_RESOLUTION | 35 |
| L_TIME_STEP | 0.0008 |
| L_BX | 4 |
| L_BY | 3.4 |
| L_BZ | 4 |
| L_RHO_IN | 1.0 |
| L_RE | 1.0 / 0.000261 |
| L_WALLS_ON | Defined |
| L_WALLS_ON_2D | Defined |
| L_WALL_THICKNESS | L_COARSE_SITE_THICKNESS |
The parameters for the LBM simulation obtained with the input data are:
- Lattice viscosity: 0.000256
- Relaxation time: 0.500767
This section presents the vertical and horizontal velocity profiles from LUMA compared to the experimental data in [1] and [3]. These results show a good qualitative agreement despite the settings used achieving a volume flow rate 8% lower than the target volumetric flow rate.
Figure 2: Mean streamwise velocity ![U/Ub] on the z-plane z/H = 0 (left) and y-plane y/H = 0
Figure 3: Streamwise Reynolds stress on the z-plane z/H = 0 (left) and y-plane y/H = 0
Figure 4: Spanwise Reynolds stress on the z-plane z/H = 0 (left) and y-plane y/H = 0
Tests will be repeated in the future with a more precise value of the forcing term and this test case updated.
[1] E. R. H. K. Meinders, “Vortex structure and heat transfer in turbulent flow over a wall-mounted matrix of cubes,” International Journal of Heat and Fluid Flow , 1999.
[2] S. T. F. Schmidt, “Comparison of numerical methods applied to the flow over wall-mounted cubes,” International Journal of Heat and Fluid Flow, 2002.
[3] R. P. Hellsten A., “8th ERCOFTAC/IAHR/COSTWorkshop on Refined TurbulenceModelling,” 8th ERCOFTAC/IAHR/COST - Workshop on Refined Turbulence Modelling, 1999.
[4] H. G. S. S. L. L. S. Yu, “DNS and LES of decaying isotropic turbulence with and without frame rotation using lattice Boltzmann method,” Journal of Computational Physics, 2005.