IMPLEMENTATION OF A MODIFIED LID DRIVEN CAVITY IN OpenFOAM ﬁ

. The standard lid-driven cavity test case is one of the most used validation cases in CFD. Whilst comparisons with experimental and particularly DNS simulations are possible, there is no an- alytical solution, and the case is ill-posed when considering the boundary conditions. A modified lid driven cavity (MLDC) case exists in the literature [1] in which the lid velocity is non-uniform and which introduces a spatially varying body force, and for which there is a closed-form analytical solution to the Navier-Stokes equations which is a function of the Reynolds number. In this paper I present an im- plementation of the MLDC as a modification of the standard OpenFOAM case, using run time coding for the boundary conditions and fvOptions , and show how convergence to the solution is affected by numerical parameters of simpleFoam such as choice of matrix inversion. The existance of an analytical solution also allows the investigation of the relation between the solver residual and the true solution error.


Introduction
Canonical test cases such as the backward facing step [2] or flow around a square prism [3] occupy an important position in CFD. They provide flow conditions which are reasonably simple to set up (particularly in terms of the geometry), which illustrate particular types of flow, such as shear or recirculation, and for which there is known data, typically experimental data, to compare against. One of the most classic canonical cases is the lid driven cavity [4,5,6,7], which generates a large primary vortex in the middle of the domain. In the simplest case, which is of course available in the OpenFOAM [8] tutorials, this is a 2d problem solved on a square domain (typically of unit side), with one moving side, usually the top side, and three stationary walls. Depending on the Reynolds number (based on lid speed and box dimension), this exhibits a range of laminar and turbulent flow, for which there is experimental data available. Given the relative simplicity of the domain, this has also been used extensively for DNS simulation, and this also provides significant validation resources [9].
More complex versions and extensions of the case have also been explored, often to address the conceptual challenges of the simple setup. Extending the case into 3d introduces a variety of additional behaviour [10] including potential periodicity in the third dimension. The behaviour of the flow in the corners of the domain can also be investigated [11]. Under the correct circumstances, the flow exhibits small secondary vortices in the corners which are not necessarily picked up by lower fidelity simulations, and some authors have explored this for example using triangular domains [12,13].
The challenge of the corner flow is made more complex by a fundamental contradiction in the simple case, namely what is the motion of the corner vertices of the domain? For the corners between the moving and stationary walls, this is ambiguous; if they count as part of the moving wall then they share that velocity in the direction of movement of the wall, but if they are thought of as part of the stationary walls then the velocity is zero. This can be resolved by using a varying lid velocity which drops to zero at the corners. This was introduced by Shih et al [1], with a formulation which also provides for a closed form analytical solution. Whilst complex, this has obvious benefits as a test case, as recognised by Marchi et al [14]. Implementing this involves creating a new boundary (Dirichlet condition but spatially varying) and introducing a complex body force into the momentum equation. This can be done by extensively modifying the existing solvers [15], but this of course creates a new solver that has to be maintained. OpenFOAM's capacity for run time coding in input files, function objects and fvOptions provides an alternative approach of implementing the necessary coding in the case directory, and this is the subject of this paper. The existence of an actual analytical solution to the Navier-Stokes equations makes this an interesting validation test case in its own right, but also allows us to investigate the relationship between the solution residual typically used in CFD (which measures the numerical imbalance of the solved equations) and the actual solution error, for example defined as the integrated error over the entire domain.
1.1. Modified Equations. Shih et al [1] work in dimensionless variables, for which the Navier Stokes equations take the form where x * = x/L 0 and L 0 is the box dimension. This includes an additional body force B oriented in the j direction. On the top surface the wall velocity is The body force B has the form where G(y * ) = g(y * )g ′′′ (y * ) − g ′ (y * )g ′′ (y * ) = −24y 5 * + 8y 3 * − 4y * Here the primes represent derivatives with respect to x * and y * . The exact solution to this problem is given by Shi et al provide the information that B(0.5, 0.5, 1) = −3.356250 as an additional check on the solution.

Implementation in OpenFOAM
2.1. Dimensional equations. Implementation of these equations in OpenFOAM requires that they first be converted back to dimensional quantities, which can be achieved by multiplying through by U 2 0 /L 0 (dimension [LT −2 ]), where U 0 is the maximum lid velocity. In particular this gives B(x, y, together with the solution For notational convenience, the non-dimensional coordinates x * , y * have been kept here, however the actual coordinates x, y must first be converted to non-dimensional form when evaluating these expressions. After this, the equations (14) -(18) can be implemented as follows : • Equation (14) can be implemented as a coded boundary condition in the velocity field. The coding for this is reproduced in Table 1. • Equation (15) can be implemented as a coded fvOption. The coding for this is reproduced in Table 3. • The solution (equations (16) -(18)) can be implemented and compared with the computed solution through a user coded functionObject. This is discussed in section 2.5. Note that the evaluation of the B term and of the overall analytical solution involves evaluation of a number of derivative terms which are common to both sets of equations, so the coding of these functions is placed in a separate file shiEqn.H which can be included into the fvOptions and functionObjects definitions as required. Evaluation of these quantities like this involves repeated creation and destruction of the appropriate geometricField objects, which is hardly efficient, but this is not an important consideration here.

2.2.
Force on the lid. In fact we can take the analysis further and evaluate the force on the moving boundary. Taking the derivative of equation (16) ∂U evaluated at y * = 1 which is the moving boundary. Since for an elemental area on the moving lid A = dx × δz, where δz is the cell thickness in the mesh, this gives which can be compared with the force evaluated from the functionObject forces; with the rho entry used to set the fluid density. Alternatively we can rearrange this to formulate a dimensionless force coefficient by dividing through by 1 2 ρU 2 0 , giving For numerical convenience, U 0 was chosen as 3 m/s, giving theoretical forces as decimals of 8, and the viscosity varied to give Re = 3, 30, 300 as the test cases.
2.3. Coded boundary. The code section in Listing 1 calculates the position-dependent speed of the moving lid. Parameters of the MLDC (boxSide ≡ L 0 and boxVel ≡ U 0 ) are provided in the transportProperties dictionary. Since dimension checking is provided at the geometricField level not at the level of boundary fields, the values but not the dimensions of these quantities are needed. Note that the boxSide parameter is linked to the geometry of the domain, so if this parameter is changed (in the separate file shiEqn.H) then blockMesh needs to be rerun on the case. §  (3) provides the additional source term Equation (15) in the momentum equation. The majority of the code is actually coded in a separate file shiCalc.H as it needs to be shared with the functionObject trueSoln (section 2.5); codeOptions being set appropriately so that the file is available at compilation. Listing (5) presents the implementation of equations (5) -(10), with the Reynolds number being calculated and output. (1) to calculate the true field solutions for p and U, and (2) to calculate the error in the calculated fields, defined as the cell weighted average of the error term normalised by the cell weighted average of the true solution, i.e.
As will be seen later, the one problem with this normalisation is that the pressure field ranges from negative to positive, and so the denominator in equation 23 and 25 may be quite small (and sensitive to other aspects of the solution), giving an unreasonably large p err normalisation.
2.6. Numerical Parameters. OpenFOAM, being an implicit CFD code, uses matrix inversion techniques to solve the individual equations. Depending on the exact equation being solved, and significantly its mathematical properties (such as symmetry/asymmetry), different inversion/solver algorithms are used, particularly between the pressure equation and the other equations. In OpenFOAM, the pressure equation can be solved either using a variant of Algebraic Multigrid (GAMG) or the preconditioned Conjugate Gradient solver (PCG), whilst the velocity equation (and other transport equations) use either the direct solver smoothSolver or the preconditioned Biconjugate Gradient solver (PBiCGStab). In the installation tutorials, the pairing GAMG/smoothSolver is typically used (particularly for the existing Lid Driven Cavity cases) and is denoted set A in table 1, whilst the combination PCG/PBiCGStab is denoted set B. Table 1 also sets out the basic numerical parameters used, which are the common defaults (except in section 6 where the tolerances were tightened to ensure that the matrix continued to solve properly).
Other numerical parameters and differencing schemes are as standard from the tutorial files; the intention here is to illustrate the relationship between matrix residual and solution error for the commonly used solver settings, not to perform an in-depth analysis of the matrix methods themselves. Similar plots can be produced for the other Reynolds number cases, but have not been included to save space.
To identify actual discrepancies between the calculated and theoretical fields, error fields can be evaluated, normalised by the cell weighted values. This has been done and the calculated fields and error fields for pressure and velocity magnitude are displayed in figures 2 and 3 respectively. Note that since the range of p goes from negative to positive, the average value is considerably smaller than the extremum values, so this normalisation in fact exaggerates the level of the error.  . v a l u e ( ) ; 3 s c a l a r UWeight = magSqr ( UTheor ) ( ) 4 . weightedAverage ( mesh ( ) .V( ) ) . v a l u e ( ) ; 5 6 const   easy comparison, the forces are normalised by the appropriate value of F theor calculated from equation (2.2). Convergence for the cases Re = 3, 30 seems to follow very much the same pattern, however the calculated force for the Re = 300 case seems to overshoot before returning. The simulation terminates when the matrix residuals drop below 10 −4 ; it is fairly obvious from this graph that the forces are still evolving at this stage and a final solution has not been achieved. Figure 4.b. shows the correlation between the solution errors for the dependent variables and the matrix residuals for the same variables. Here the solid lines give the correlation between the pressure error p err (as defined by equation (23) and the initial pressure residual (p 0 extracted from running foamLog). The dashed and dash-dot lines demonstrate the correlation between U err (defined by equation (24) and the Ux 0 and Uy 0 residuals. As can be seen there are substantial differences between the a. b. c. d.
e. f.  pressure curves although some of this may be down to overall normalisation -the pressure field of course ranges from negative to positive so the denominator in equation 23 is possibly small. The other curves relating U err to the U residuals all more or less follow the same path. What is interesting here is how the solution errors change as the calculation proceeds, i.e. as the matrix residuals reduce in magnitude. This corresponds to tracing the curves from right to left as the SIMPLE loop executes. As can be seen, initially the curves are quite flat; improving the matrix residual 1 → 0.1 → 0.01 has comparatively little effect on the solution accuracy. The solid black line on the figure shows a 1:1 correlation between the Error and Residual, which would be ideal. The actual solution only really starts to substantially improve when the matrix residuals drop to < 10 −3 , and as observed above, the errors are still reducing when the matrix residuals reach 10 −4 and the algorithm terminates. Figure 5 shows the same graphs for the Conjugate Gradient solution. Largely this shows the same behaviour. Convergence in this case is slightly quicker (in terms of number of iterations of the SIMPLE loop) but not significantly so, and roughly the same comments may be made about the relationship between matrix residuals and the solution errors.
3.3. Richardson Extrapolation. Figure 6 (top) presents results for the different Reynolds numbers based on convergence of the force results. Instead of stopping the simulation when the matrix residuals drop below 10 −4 , the simulation was run until the calculated lid force reached a stable value. This was achieved for the coarse meshes to 6 significant figures but required significantly more effort for the finer meshes (N = 320 meshes required 20,000 iterations to converge to achieve this). Tolerances on the matrix solvers were tightened to 1 × 10 −9 to ensure that the matrices continued to solve throughout; otherwise the matrix combination A was used. Richardson extrapolation was then applied by fitting a quadratic curve to the data for each Reynolds number and then extrapolating to a zero reduced mesh spacing (the mesh spacing for N = 320 being 1.0) to find the infinite mesh resolution value for the force. Table 2 shows the zero value for the different Reynolds numbers, indicating an error in this parameter < 0.1%.
The results can also be analysed to calculate the actual numerical error in the algorithm. As is well known, the error term for any numerical scheme can be written as where h is the cell dimension and c the order of the scheme. Since the standard settings (used here) for simpleFoam comprise a blend of 1st and 2nd order numerics we expect the effective value of c to lie 1 < c < 2. In order to compare the values for all three Reynolds numbers we can divide this expression through by F theor and evaluate  Figure 7 shows the evolution of both matrix residuals and solution errors against iteration for the three different Reynolds number cases. Judging by the solution errors the simulations have converged at 2500 iterations and no improvement in solution accuracy is possible, although the matrix residuals continue to decrease until they reach the set level of 10 −9 . The error in U reaches a level of 10 −7 − 10 −8 ; however the error in force F and the pressure error remain orders of magnitude higher. Since it is the viscous force being evaluated (related to the velocity gradients) these two errors are not obviously related. Figure 8 shows the correlation between error and residual for different quantities; errors in force and pressure plotted against pressure residual (top graph) and error in U plotted against U x , U y residuals (bottom graph, A curves). Only the first 2500 iterations (up to convergence) are plotted. The straight lines give a 1:1 relationship for comparison. As before (section 3.2, figure 4), initial improvements in residual have little effect on the solution; residual convergence beyond 10 −3 is more effective, whilst the final residual convergence beyond 10 −6 is ineffectual again. The somewhat odd behaviour of the force error for the Re = 300 case is because the error goes negative for a time (the absolute error is being plotted here).
Since the force on the lid depends on the gradient of the velocity, figure 8 also shows the correlation between the error in F and the U y and U x residuals (solid and dashed lines, respectively, lower figure, B curves). As expected the correlations are somewhat smoother here; the lines for Re = 3 and Re = 30 are basically coincident, whilst those for Re = 300 diverge somewhat as the force error goes negativeagain, the absolute error is being plotted (in order to use logarithmic axes).

Conclusions
Validation and verification are important aspects of CFD. The distinction between these activities can be summed up as : verification is the act of checking that the mathematical equations have been validly implemented in the code, whilst validation is the act of checking that the solution matches mathematical or physical reality -that the mathematical model is itself correct. In both activities, canonical cases, which illustrate specific types of flow in simple geometries, hold a significant position. Because of the complex nature of the Navier Stokes equations, mathematical solutions are rarely available, so experimental or DNS numerical data has often been used. However where analytical results are possible they can provide a very high quality of comparison. Juretic [16] for example has demonstrated validation of OpenFOAM against algebraic solutions for a planar jet, as well as analytical solutions for simpler physics such as creeping flow and heat conduction.
This paper presents an implementation in OpenFOAM of a modified Lid Driven Cavity test case. The modifications, due to Shi et al [1,14] include a non-uniform lid velocity (thus avoiding the ambiguity in the velocity at the corners of the domain) and a body force, and allow a closed form analytical solution. This has been extended to calculate the force on the lid, which provides a convenient physical and quantitative comparison between the solution and the numerical results. In addition to validating the classic simpleFoam solver for both multigrid and Conjugate Gradient solvers, this provides an excellent test case to examine the relation between the solution residual from the matrix solver and the actual mathematical error in the results. This shows the following : • Monitoring residuals alone is not necessarily sufficient to identify convergence of a simulation. In particular even a relatively tight residual tolerance (of 10 −4 ) does not guarantee convergence of physical properties such as boundary forces. • In general the relationship between the matrix residuals and the physical solution error is monotonic but not linear. Initial reductions in residual may not be accompanied by equivalent reductions in error. Where there is a mathematical link between the quantity being solved and the physical parameter being monitored (so between the velocity field and the force which is a function of the gradient of that field) there may be a closer relationship between the residual and the error. • Fully converging the forces and applying Richardson extrapolation across a range of mesh resolutions, the error in the solution can be reduced to less than 0.1%, with a numerical error order of around 1.5, indicating a combination of 1st and 2nd order schemes. Further tuning of the numerical parameters might improve this, and this could be a valuable test case for comparing the actual effects of various numerical schemes (both differencing schemes and matrix solvers).