From 0c8ca47199e0f67b9698e3fdc8b674e63990784a Mon Sep 17 00:00:00 2001
From: henry <Henry Weller h.weller@opencfd.co.uk>
Date: Tue, 18 Aug 2009 17:13:22 +0100
Subject: [PATCH] Added support for variable density incompressible solvers.

---
 .../functionObjects/forces/forces/forces.C    | 75 ++++++++++++++++---
 .../functionObjects/forces/forces/forces.H    |  6 ++
 2 files changed, 71 insertions(+), 10 deletions(-)

diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index 1ffc36cab0..6baace3863 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -61,7 +61,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
         const incompressible::RASModel& ras
             = obr_.lookupObject<incompressible::RASModel>("RASProperties");
 
-        return rhoRef_*ras.devReff();
+        return rho()*ras.devReff();
     }
     else if (obr_.foundObject<compressible::LESModel>("LESProperties"))
     {
@@ -75,7 +75,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
         const incompressible::LESModel& les
             = obr_.lookupObject<incompressible::LESModel>("LESProperties");
 
-        return rhoRef_*les.devBeff();
+        return rho()*les.devBeff();
     }
     else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
     {
@@ -97,7 +97,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 
         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 
-        return -rhoRef_*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
+        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
     }
     else if (obr_.foundObject<dictionary>("transportProperties"))
     {
@@ -108,7 +108,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 
         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 
-        return -rhoRef_*nu*dev(twoSymm(fvc::grad(U)));
+        return -rho()*nu*dev(twoSymm(fvc::grad(U)));
     }
     else
     {
@@ -121,6 +121,34 @@ Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::forces::rho() const
+{
+    if (rhoName_ == "rhoInf")
+    {
+        const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+        return tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    mesh.time().timeName(),
+                    mesh
+                ),
+                mesh,
+                dimensionedScalar("rho", dimDensity, rhoRef_)
+            )
+        );
+    }
+    else
+    {
+        return(obr_.lookupObject<volScalarField>(rhoName_));
+    }
+}
+
+
 Foam::scalar Foam::forces::rho(const volScalarField& p) const
 {
     if (p.dimensions() == dimPressure)
@@ -129,6 +157,13 @@ Foam::scalar Foam::forces::rho(const volScalarField& p) const
     }
     else
     {
+        if (rhoName_ != "rhoInf")
+        {
+            FatalErrorIn("forces::rho(const volScalarField& p)")
+                << "Dynamic pressure is expected but kinematic is provided."
+                << exit(FatalError);
+        }
+
         return rhoRef_;
     }
 }
@@ -149,11 +184,12 @@ Foam::forces::forces
     active_(true),
     log_(false),
     patchSet_(),
-    pName_(""),
-    UName_(""),
+    pName_(word::null),
+    UName_(word::null),
+    rhoName_(word::null),
     directForceDensity_(false),
     fDName_(""),
-    rhoRef_(0),
+    rhoRef_(VGREAT),
     CofR_(vector::zero),
     forcesFilePtr_(NULL)
 {
@@ -175,6 +211,12 @@ Foam::forces::forces
     }
 
     read(dict);
+
+    if (active_)
+    {
+        // Create the forces file if not already created
+        makeFile();
+    }
 }
 
 
@@ -222,18 +264,31 @@ void Foam::forces::read(const dictionary& dict)
             // Optional entries U and p
             pName_ = dict.lookupOrDefault<word>("pName", "p");
             UName_ = dict.lookupOrDefault<word>("UName", "U");
+            rhoName_ = dict.lookupOrDefault<word>("rhoName", "rho");
 
-            // Check whether UName and pName exists, if not deactivate forces
+            // Check whether UName, pName and rhoName exists,
+            // if not deactivate forces
             if
             (
                 !obr_.foundObject<volVectorField>(UName_)
              || !obr_.foundObject<volScalarField>(pName_)
+             || (
+                    rhoName_ != "rhoInf"
+                 && !obr_.foundObject<volScalarField>(rhoName_)
+                )
             )
             {
                 active_ = false;
+
                 WarningIn("void forces::read(const dictionary& dict)")
-                << "Could not find " << UName_ << " or "
-                    << pName_ << " in database." << nl
+                    << "Could not find " << UName_ << ", " << pName_;
+
+                if (rhoName_ != "rhoInf")
+                {
+                    Info<< " or " << rhoName_;
+                }
+
+                Info<< " in database." << nl
                     << "    De-activating forces."
                     << endl;
             }
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.H b/src/postProcessing/functionObjects/forces/forces/forces.H
index 96651af97e..bb82aac412 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.H
+++ b/src/postProcessing/functionObjects/forces/forces/forces.H
@@ -132,6 +132,9 @@ protected:
             //- Name of velocity field
             word UName_;
 
+            //- Name of density field (optional)
+            word rhoName_;
+
             //- Is the force density being supplied directly?
             Switch directForceDensity_;
 
@@ -157,6 +160,9 @@ protected:
         //- Return the effective viscous stress (laminar + turbulent).
         tmp<volSymmTensorField> devRhoReff() const;
 
+        //- Return rho if rhoName is specified otherwise rhoRef
+        tmp<volScalarField> rho() const;
+
         //- Return rhoRef if the pressure field is dynamic, i.e. p/rho
         //  otherwise return 1
         scalar rho(const volScalarField& p) const;