From 2d8537759f22e8ed76cf12e7d189d64c42895c4a Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 24 Nov 2022 08:43:09 +0000
Subject: [PATCH] driftFluxFoam: Added drift velocity to Courant number test

to aid stability in cases where the drift velocity is comparable to the mean
velocity.
---
 .../multiphase/driftFluxFoam/CourantNo.H      | 53 +++++++++++++++++++
 1 file changed, 53 insertions(+)
 create mode 100644 applications/solvers/multiphase/driftFluxFoam/CourantNo.H

diff --git a/applications/solvers/multiphase/driftFluxFoam/CourantNo.H b/applications/solvers/multiphase/driftFluxFoam/CourantNo.H
new file mode 100644
index 0000000000..2f18166804
--- /dev/null
+++ b/applications/solvers/multiphase/driftFluxFoam/CourantNo.H
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    CourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar CoNum = 0.0;
+scalar meanCoNum = 0.0;
+
+{
+    const scalarField sumPhi
+    (
+        fvc::surfaceSum
+        (
+            mag(phi) + mag(fvc::flux(mixture.Udm()))
+        )().primitiveField()
+    );
+
+    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Courant Number mean: " << meanCoNum
+    << " max: " << CoNum << endl;
+
+// ************************************************************************* //