The validation of forces on boundaries (section 3.6) is done with a problem to which the analytic solution is known – the Stokes drag force on a sphere in laminar flow. Figure19illustrates what is being simulated in this section.
Figure 19: Laminar flow past a sphere of radius R = 16 in a cubic simulation box of side length L = 384 visualized with streamlines. The flow direction is left-to-right, although due to the symmetry atRe= 0.01 one could not tell from the image. The coloring indicates velocity. Visualization is done with the real-time OpenCL graphics engine ofFluidX3D.
9.2.1 Analytic Solution
For a sphere of radiusRresting at the coordinate origin, when it is in laminar flow with densityρ0and velocity~u0
at infinite distance from the origin, the analytic density and velocity fields [88, p.230-235][89][90, p.168-171][91, p.36-38] are known as the following:
~
x:= (x, y, z)T r:=p
x2+y2+z2 (181)
ρ(~x) =ρ0−9ρ0ν R(~u0◦~x)
2r3 (182)
~
u(~x) =~u0−3 4
R r + R3
3r3
~ u0+
R r3 −R3
r5
(~u0◦~x)~x
(183) The analytic solution for the drag forceF~ on the sphere is given by
F~ = 6π ρ ν R ~u0 (184)
9.2.2 Strategy
The goal of this test is to validate the accuracy of boundary forces in LBM (section3.6), meaning that the error must explicitly be the error of the force on the sphere and not the error of the velocity field. This in turn means that a volume force (section3.5) cannot be used to drive the flow, because the sphere is the only object in the simulation box that slows down the flow. When a volume-force-driven flow becomes stationary, any boundary geometry in the simulation box will always absorb exactly the amount of force that was put into the fluid via
9 Error Validation High Performance Free Surface LBM on GPUs
volume force, making the error always (almost) zero. Of course, one could then measure the stationary velocity profile and calculate the error from velocity, but that would only validate the bounce-back boundaries and not fulfill the purpose of validating the boundary forces.
Instead, the flow needs to be driven by moving bounce-back boundaries (section3.4.3) at the outer edge of the simulation box, which in turn have been validated with Poiseuille flow. Furthermore, it is not sufficient to set the velocity at these boundaries to the velocity~u0 infinitely far away from the sphere, because the simulation box can never be infinitely large – the error would then converge to about 10 % (force would be too large) even for the largest possible simulation box that fits into memory. Even in an infinitely large simulation box, the streamlines are never perpendicular to~u0; they always curve around the sphere, although curvature gets less further out. It would not make sense to enforce straight streamlines at finite distance to the sphere, thereby artificially constricting the flow. To avoid this, the velocity at the boundaries is set to the analytic solution~u(~x) from equation (183), allowing streamlines to curve out and back into the simulation box at the boundaries.
9.2.3 Error Definition and Convergence Criteria
The relative errorEof a measured propertyxsimof the simulation to the theoretically expected propertyxtheo
here is calculated as theL1 norm:
E(x) :=|xsim−xtheo| xtheo
(185) Since it is not clear that the error will only decrease during simulation, it is not sufficient to run the simula-tion until the error reaches a local minimum. Instead, the following definisimula-tion is used for determining error convergence: Both the absolute slope and the absolute curvature
of the error E must be smaller than E with 1 being a small number. The derivatives are calculated as the first- and second-order backward difference of the last three error values, which are computed everyN simulation time steps.
9.2.4 Simulation Parameters
All simulations are performed with D3Q19 TRT at Reynolds numberRe = 0.01. This leaves the freedom of choosing either the kinematic shear viscosityν or velocityu0:=|~u0|. Simulations have shown (figure20) that the errorE does not depend onν as long as 0< ν ≈100 and 0.0003≤ |~u0| ≤0.5 remain in reasonable ranges (figure16). For section9.2.5,ν = 1 (except forR= 32 whereν = 10) and for9.2.6, ν= 10 is used in order to not have too small velocity. The error is calculated as E(|F~|) every N = 100 LBM time steps. For achieving sufficient convergence,= 10−4 is chosen. The sphere radius is denoted as Rand the simulation box is cubic with side lengthL. The flow is created along thex-axis with positive velocity.
0 accordingly. The observation is that as long as no physical property is changed (Re=const), the error does not significantly change as well.
9 Error Validation High Performance Free Surface LBM on GPUs
9.2.5 Results –R=const,L is varied
Figure21shows the error behavior when the sphere radius is kept constant and the simulation box size is varied.
As expected, for a larger simulation box, the error goes down to a plateau induced by the staircase effect of sphere voxelation.
0 1 2 3 4 5 6 7 8 9 10
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64
E(F) / %
L / R
R = 4 R = 8 R = 16 R = 32
Figure 21: For a fixed sphere radius R ={4, 8,16,32} the box size is varied. As expected, for a larger box size L the error decreases down to a plateau. Moreover, the error is smaller for larger spheres, indicating the plateau being caused by the staircase-effect of the sphere in limited voxel resolution.
9.2.6 Results –L/R=const, R is varied
Figure22shows various simulations for exactly the same physics (L/R=const) with varying voxel resolution of the sphere. As expected, the error overall decreases with increasingR(voxel resolution), with the addition that the error is especially small whenRis a multiple of 12 – for which there surely is some geometrical explanation.
0 1 2 3 4 5 6 7 8 9 10
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64
E(F) / %
R / Δx
L = 4 R L = 8 R L = 16 R L = 32 R
Figure 22: Here the ratio of the simulation box size and the sphere radiusL/Ris fixed andRis varied, meaning that along a data line the simulated physics are identical. Confirming the indication of figure21that the plateau of the error is caused by the staircase-effect; here the observation is that for identical physics the error decreases when the sphere is better resolved (Ris larger). The data points are not monotonic however and certain integer radii (multiple of 12) will make the error especially small. Also in agreement to figure21, larger L/Rresults in a smaller error.
9 Error Validation High Performance Free Surface LBM on GPUs
9.2.7 Results – Velocity Field Errors
Figure23shows the velocity field around the sphere as well as its error for various sphere radii. The simulated velocity fields are indistinguishable from the theoretical solution with the bare eye; only separate plots of the error show that the error in velocity is largest in close vicinity to the sphere surface.
-128
Figure 23: A slice through the velocity field of laminar flow past a sphere of radiusR={4,8,16,32} (top to bottom) at Re = 0.01. The flow direction is from left to right. The columns are the normalized theoretical velocity magnitude |~utheo|/u0 (left), the normalized simulated velocity magnitude |~usim|/u0 (middle) and the error of the velocity magnitude E(|~u|) (right). For the first three rows, ν = 1 and for the last row ν = 10 in order to avoid the floating-point errors of too small fluid velocity. The error is largest in close vicinity to the sphere surface, where the velocity magnitude is smallest, and its distribution is not exactly symmetrical due to Re >0. Next to the voxelated sphere surface, single voxels with unusually high error are observed due to the staircase-effect.
9 Error Validation High Performance Free Surface LBM on GPUs
Figure 24: A slice through the velocity field of laminar flow past a sphere of radius R= 16 at Re= 0.01. The flow direction is left to right, indicated with arrows. The theoretical velocity field is plotted on the top and the simulated velocity field is plotted in the middle. On the bottom, the relative difference of the theoretical and simulated velocity fields is plotted as a vector field. Errors are largest in close vicinity of the sphere surface.
9 Error Validation High Performance Free Surface LBM on GPUs