A curiously useful old formula

A curiously useful old formula

David Standingford, Lead Technologist, CFMS, talks about the usefulness of mesh convergence studies – or more strictly mesh independence studies, commonly used in 1990’s when numerical fluid dynamics simulation based on the available computing power was slow.

David CropMy wife and I have almost managed to convince our kids (10 and under) that the world prior to 1990 was all in black and white.  The books that they see from before this ancient time are largely colourless, and the palettes used in films look as though they were touched up after modern light became available.  Despite this, they refuse to believe that it was possible to exist without wireless networking and personal tablet computers – preferably expensive ones.

In the sepia-toned world of 1990, numerical fluid dynamics simulation based on the available computing power was slow and the challenges tackled were often straight extensions of the fundamental analytical work of the great exponents: Taylor, Batchelor, Lamb, Michell, Prandtl and the ilk.  It was possible and expected that deviation from well established results could be corrected with hard thought and judicious application of calculation – even then sometimes via specialist departments who could turn algorithms into computer code.  

When computing resource is scarce and simulations are limited to relatively benign linear or quasi-linear regimes, rigor is expected on all fronts.  Mesh convergence studies – or more strictly mesh independence studies – demonstrating that the solution obtained is not merely some artefact resulting from error cancellation, were seen as a basic requirement for any credible contribution to the literature.  Most graduate-level computational fluid dynamics (CFD) courses included this as a tenet and woe betide any study omitting a convergence report.   With the exception of some international workshops (e.g. the NASA aerodynamics conferences) it has been a long time since I saw this kind of requirement.  Instead, simulation data is piled on top of other sources of confidence in design work – somehow creating a warm feeling that everything is OK because the experimental trends broadly align with the computational ones.  

Perhaps this is because the simulation tasks are now so critically dependent upon the specific meshes that support them that mesh independence is an abandoned objective.  It might also be that mesh refinement using the straightforward operation of doubling the mesh density in each of the three dimensions results in meshes that are just too big to be used to obtain a solution.  This is especially true if the original mesh was sized to use the available resource to start with.

Many years ago I inherited a piece of FORTRAN code designed to calculate the linearised lift coefficient of a simple planar wing geometry using a variant of the vortex lattice method (developed by Lan).  Within the code, a simple formula has been used to combine the solutions from three separate calculations, using 17x17, 18x18 and 19x19 cells to represent the geometry.  By any modern standard, this is a trivial calculation but at the time the inversion of a 361x361 matrix (which is what the largest mesh results in) was hard work for a desktop computer.  The results from the three calculations (none of which was particularly close to the correct analytical answer) were labelled F1, F2 and F3, respectively.  A fourth result was determined as:

F4 = (F1*F3  – F2**2) / (F1+F3 - 2.0D0*F2) (1)

Note that in FORTRAN, (**) means exponent and the ‘D0’ after the 2 is often added for good measure to ensure that calculations are executed with a particular precision.  I was struck by the simplicity of the expression, and the astonishing accuracy of the resulting solution F4 (correct to 5 significant figures).

The extrapolation formula above works because the answer produced by the code is equal to the correct answer plus an error that varies systematically with mesh size.  As a formula, the code output can be written as: 

F (n) = A + B*(n**C)(2)

Where A is the real answer, n is an input parameter (in the above example 17, 18 or 19 corresponding to the number of mesh points in each dimension) and B, C are fixed parameters.  Providing (C < 0) then the expression (2) will converge with increasing n to the correct answer A for any value of B.  In order for formula (1) to work, it is only necessary that you choose successive values of n1, n2 and n3 (for F1, F2 and F3 respectively) so that (n3-n2) = (n2-n1).  You do not need to know the values for A, B or C. This is why the values of n1, n2 and n3 in the piece of code were chosen as 17, 18 and 19 respectively.

For the mathematically minded, it is fun to prove that formula (1) works.  Start with expression (2) for each of n1, n2 and n3 and observe that for the natural logarithm:

Log (1+x) ~ x, for small values of x.

For the results-oriented, observe that if you apply formula (1) to your own work and get a result F4 that is more accurate than the input values F1, F2 and F3 then you can have some confidence that your results would pass the mesh independence test.  If F4 is not an improvement, then there is reason to suspect that your answer is some artefact resulting from …

Well – you know the rest!

Screen Shot 2016 06 30 At 101057

Sign Up To Our Newsletter and Events