From b5142dca74d5d0bbb62b6a725806a793960c2f7b Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 9 May 2023 14:56:56 +0100 Subject: [PATCH] modules/solid: Improved diffusion number calculation The calculation of the diffusion number has been put into a form consistent with finite-volume, rather than the finite-difference form that was used previously. This difference in formulations is analogus to that of the Courant number in the fluid solvers. Whilst a textbook will typically define the courant number as equal to 'U*deltaT/deltaX', in a finite-volume context it is more appropriate to define it as 'Sum(phi)/V*deltaT' (where 'Sum' is a sum over the cell's faces). Similarly, the finite-difference Fourier number, 'kappa/rho/Cp*deltaT/deltaX^2', is more consistently expressed for finite-volume as 'Sum(Sf*kappa*deltaX)/(V*rho*Cp)*deltaT'. This makes the calculation of the diffusion number less sensitive to the presence of small, poor quality faces, and therefore makes time-step adjustment more robust on arbitrary polyhedral meshes. --- applications/solvers/modules/solid/solid.C | 23 ++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/applications/solvers/modules/solid/solid.C b/applications/solvers/modules/solid/solid.C index 3bf01eac02..e2755de4c3 100644 --- a/applications/solvers/modules/solid/solid.C +++ b/applications/solvers/modules/solid/solid.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "solid.H" +#include "fvcSurfaceIntegrate.H" #include "fvMeshMover.H" #include "localEulerDdtScheme.H" #include "addToRunTimeSelectionTable.H" @@ -61,19 +62,25 @@ void Foam::solvers::solid::correctDiNum() : mag(thermo.Kappa())() ); - const surfaceScalarField kapparhoCpbyDelta + const volScalarField::Internal DiNumvf ( - sqr(mesh.surfaceInterpolation::deltaCoeffs()) - *fvc::interpolate(kappa) - /fvc::interpolate(thermo.rho()*thermo.Cp()) + fvc::surfaceSum + ( + mesh.magSf() + *fvc::interpolate(kappa) + *mesh.surfaceInterpolation::deltaCoeffs() + )()() + /(mesh.V()*thermo.rho()()*thermo.Cp()()) + *runTime.deltaT() ); - DiNum = max(kapparhoCpbyDelta).value()*runTime.deltaTValue(); - const scalar meanDiNum = - average(kapparhoCpbyDelta).value()*runTime.deltaTValue(); + const scalar meanDiNum = gAverage(DiNumvf); + const scalar maxDiNum = gMax(DiNumvf); Info<< "Diffusion Number mean: " << meanDiNum - << " max: " << DiNum << endl; + << " max: " << maxDiNum << endl; + + DiNum = maxDiNum; }