View Issue Details

IDProjectCategoryView StatusLast Update
0002453OpenFOAMBugpublic2017-02-13 18:31
ReporterSahas Assigned Tohenry  
PrioritynormalSeveritytweakReproducibilityalways
Status resolvedResolutionfixed 
Fixed in Versiondev 
Summary0002453: Solver postProcess with scalarTransport seems not very useful
DescriptionFirstly, I try to calculate scalar transport after main computation finish:
simpleFoam -postProcess -func scalarTransport -time 1000
I see an output:
scalarTransport write:
DILUPBiCG: Solving for s, Initial residual = 0.999999, Final residual = 0.0728566, No Iterations 1
But there is no output of field "s" with actual values.

Secondly, for scalar transport in turbulent flow one should specify both diffusion coefficient D and turbulent Schmidt number Sct, because now turbulent scalar transport is incorrect: only nuEff() is used. I think, one should rewrite function Foam::functionObjects::scalarTransport::D in src/functionObjects/solvers/scalarTransport/scalarTransport.C in such a similar way:
    if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
    {
        const icoModel& model = mesh_.lookupObject<icoModel>
        (
            turbulenceModel::propertiesName
        );
        return model.nut()/Sct + D;
    }
Steps To ReproduceRun tutorial in tutorials/incompressible/pisoFoam/les/pitzDaily. After outputing in time 0.001 run:
pisoFoam -postProcess -func scalarTransport -time 0.001
There is no updated field s anywhere.
TagsNo tags attached.

Activities

henry

2017-02-06 15:32

manager   ~0007710

Running the "scalarTransport" functionObject as a post-processing function doesn't make a lot of sense because it evolves the scalar field over time not as a single-shot solution to convergence. Try the tutorials/incompressible/pisoFoam/LES/pitzDaily case and you will see the "scalarTransport" functionObject in operation. The same should be possible with simpleFoam.

In principle it would be possible to write a steady-state version of the "scalarTransport" functionObject to use on the results of a simpleFoam simulation.

For a turbulent simulation the diffusivity used in the "scalarTransport" functionObject could be set to an effective turbulent diffusivity but often users want to see the transport of scalars through the domain without dispersion and set D to 0 so the application of the effective turbulent diffusivity would have to be a run-time selectable option.

Sahas

2017-02-06 16:05

reporter   ~0007712

Thanks for reply!
Yes, you are right: I can run scalarTransport together with main solver.
But my passive scalar (that is the temperature in fact) should be transported including effects of both turbulent and molecular transport with different Prandtl numbers (0.9 and 0.7 accordingly). And unfortunately now there is no option to achieve the goal: I can use neither scalarTransportFoam (because case with constant D is only supported) nor objectFunction "scalarTransport".
Of course I can write my own solver so the problem is actually not a problem :)

henry

2017-02-06 16:24

manager   ~0007713

The simplest option would be to add support for the laminar and turbulent Schmidt numbers to the "scalarTransport" functionObject. If you choose this option please contribute this development for inclusion in OpenFOAM-dev.

henry

2017-02-08 11:00

manager   ~0007719

Please reopen if you would like to provide a patch to the scalarTransport functionObject.

Sahas

2017-02-08 14:31

reporter   ~0007720

I am working on the patch to the scalarTransport functionObject and I have found that FatalErrorInFunction do not work properly when functionObject is called during solver run. Is it correct behavior?

henry

2017-02-08 14:52

manager   ~0007721

No I have not had any problems with FatalErrorInFunction.

Sahas

2017-02-08 15:10

reporter   ~0007722

I've just added two lines in src/functionObjects/solvers/scalarTransport/scalarTransport.C:

function bool Foam::functionObjects::scalarTransport::execute():
{
    Info<< type() << " write:" << endl;
    Info<< "Test Info" << endl;
    FatalErrorInFunction << "Test FatalErrorInFunction" << endl;

I compile it with "wmake libso" and go to tutorial tutorials/incompressible/pisoFoam/LES/pitzDaily.

If I run just "postProcess -func scalarTransport" I get output:
scalarTransport write:
Test info
--> FOAM FATAL ERROR:
Test FatalErrorInFunction

But if I run "pisoFoam -postProcess -func scalarTransport" I get only
scalarTransport write:
Test info

and there is no "FATAL ERROR", the program continues to run. What am I doing wrong?

henry

2017-02-08 15:20

manager   ~0007723

You are not using FatalErrorInFunction correctly, you will need to add

    << exit(FatalError);

Take a look at any other use of FatalErrorInFunction in OpenFOAM.

Sahas

2017-02-08 16:07

reporter   ~0007724

Thank you, the problem is solved. In the original file src/functionObjects/solvers/scalarTransport/scalarTransport.C FatalErrorInFunction (line 254) was called without exit(FatalError). That why I was confused.

henry

2017-02-08 17:03

manager   ~0007725

Thanks for spotting this, I have now corrected it in OpenFOAM-4.x and OpenFoam-dev.

Sahas

2017-02-08 17:30

reporter   ~0007726

I attach patch to the scalarTransport functionObject with Sc and Sct. I kept the old behavior and just added new functionality. I performed small tests and it seems to be working (I want to make some test for my needs but it can take a long time). But the code looks like spaghetti. Please check whether it is appropriate for including in OF.
scalarTransport.patch (5,463 bytes)   
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.C b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
index 954a5d6..51fbbb1 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.C
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
@@ -82,6 +82,58 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             )
         );
     }
+    else if (laminarSc_)
+    {
+        if (phi.dimensions() == dimMass/dimTime)
+        {
+
+            const volScalarField& mu =
+                mesh_.lookupObject<volScalarField>(muName_);
+
+            tmp<volScalarField> D(new volScalarField(Dname, mu/Sc_));
+
+            if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
+            {
+                if (!turbulentSc_)
+                {
+                    FatalErrorInFunction
+                        << "Turbulent Schmidt number Sct is not specified"
+                        << exit(FatalError);
+                }
+                const cmpModel& model = mesh_.lookupObject<cmpModel>
+                (
+                    turbulenceModel::propertiesName
+                );
+                D = D + model.mut()/Sct_;
+            }
+
+            return D;
+        }
+        else
+        {
+            const volScalarField& nu =
+                mesh_.lookupObject<volScalarField>(nuName_);
+
+            tmp<volScalarField> D(new volScalarField(Dname, nu/Sc_));
+
+            if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
+            {
+                if (!turbulentSc_)
+                {
+                    FatalErrorInFunction
+                        << "Turbulent Schmidt number Sct is not specified"
+                        << exit(FatalError);
+                }
+                const icoModel& model = mesh_.lookupObject<icoModel>
+                (
+                    turbulenceModel::propertiesName
+                );
+                D = D + model.nut()/Sct_;
+            }
+
+            return D;
+        }
+    }
     else if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
     {
         const icoModel& model = mesh_.lookupObject<icoModel>
@@ -167,6 +219,8 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
 
     phiName_ = dict.lookupOrDefault<word>("phi", "phi");
     rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
+    muName_  = dict.lookupOrDefault<word>("mu", "thermo:mu");
+    nuName_  = dict.lookupOrDefault<word>("nu", "nu");
     schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
 
     constantD_ = false;
@@ -175,6 +229,30 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
         constantD_ = true;
     }
 
+    laminarSc_ = false;
+    if (dict.readIfPresent("Sc", Sc_))
+    {
+        if (constantD_)
+        {
+            FatalErrorInFunction
+                << "Both D and Sc are specified in dictionary "
+                << dict.name() << abort(FatalError);
+        }
+        laminarSc_ = true;
+    }
+
+    turbulentSc_ = false;
+    if (dict.readIfPresent("Sct", Sct_))
+    {
+        if (!laminarSc_)
+        {
+            FatalErrorInFunction
+                << "Laminar Sc must be specified in dictionary "
+                << dict.name() << exit(FatalError);
+        }
+        turbulentSc_ = true;
+    }
+
     dict.readIfPresent("nCorr", nCorr_);
 
     if (dict.found("fvOptions"))
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.H b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
index 86c92bd..8bb9441 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.H
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,8 @@ Description
       the 'schemesField' entry,
     - The diffusivity can be set manually using the 'D' entry, or retrieved
       from the turbulence model (if applicable).
+    - The laminar Sc and turbulent Sct Schmidt numbers can be specified
+      instead of diffusivity D
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
@@ -77,12 +79,30 @@ class scalarTransport
         //- Name of density field (optional)
         word rhoName_;
 
+        //- Name of dynamic viscosity field (optional)
+        word muName_;
+
+        //- Name of kinematic viscosity field (optional)
+        word nuName_;
+
         //- Diffusion coefficient (optional)
         scalar D_;
 
         //- Flag to indicate whether a constant, uniform D_ is specified
         bool constantD_;
 
+        //- Laminar Schmidt number (optional)
+        scalar Sc_;
+
+        //- Flag to indicate whether laminar Schmidt number is specified
+        bool laminarSc_;
+
+        //- Turbulent Schmidt number (optional)
+        scalar Sct_;
+
+        //- Flag to indicate whether turbulent Schmidt number is specified
+        bool turbulentSc_;
+
         //- Number of corrector iterations (optional)
         label nCorr_;
 
scalarTransport.patch (5,463 bytes)   

henry

2017-02-08 17:39

manager   ~0007727

I would rather the code did not look like spaghetti, is there any way you can improve it?

It might be preferable to use "alpha" diffusion coefficients which multiply the viscosities rather than Schmidt numbers as it is then easier to switch-off the terms by setting them to 0.

Sahas

2017-02-09 09:39

reporter   ~0007730

> I would rather the code did not look like spaghetti, is there any way you can improve it?

Is there a need to keep old behavior? With Schmidt numbers (or alphas), old behavior can be easily reproduced. I think that one should keep only two possible options: constant D or imposed alpha numbers.

henry

2017-02-09 09:49

manager   ~0007731

Agreed

Sahas

2017-02-10 13:58

reporter   ~0007735

I attached the new patch with alphaD and alphaDt introduced.
Apropos, the fatal error in functionObject does not lead to program interruption, is it correct behavior?
scalarTransport_ver2.patch (6,103 bytes)   
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.C b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
index 954a5d6..4a1d996 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.C
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.C
@@ -82,42 +82,44 @@ Foam::tmp<Foam::volScalarField> Foam::functionObjects::scalarTransport::D
             )
         );
     }
-    else if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
+    else
     {
-        const icoModel& model = mesh_.lookupObject<icoModel>
-        (
-            turbulenceModel::propertiesName
-        );
+        if (phi.dimensions() == dimMass/dimTime)
+        {
+            const volScalarField& mu =
+                mesh_.lookupObject<volScalarField>(muName_);
 
-        return model.nuEff();
-    }
-    else if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
-    {
-        const cmpModel& model = mesh_.lookupObject<cmpModel>
-        (
-            turbulenceModel::propertiesName
-        );
+            tmp<volScalarField> D(new volScalarField(Dname, mu*alphaD_));
 
-        return model.muEff();
-    }
-    else
-    {
-        return tmp<volScalarField>
-        (
-            new volScalarField
-            (
-                IOobject
+            if (mesh_.foundObject<cmpModel>(turbulenceModel::propertiesName))
+            {
+                const cmpModel& model = mesh_.lookupObject<cmpModel>
                 (
-                    Dname,
-                    mesh_.time().timeName(),
-                    mesh_.time(),
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                mesh_,
-                dimensionedScalar(Dname, phi.dimensions()/dimLength, 0.0)
-            )
-        );
+                    turbulenceModel::propertiesName
+                );
+                D = D + model.mut()*alphaDt_;
+            }
+
+            return D;
+        }
+        else
+        {
+            const volScalarField& nu =
+                mesh_.lookupObject<volScalarField>(nuName_);
+
+            tmp<volScalarField> D(new volScalarField(Dname, nu*alphaD_));
+
+            if (mesh_.foundObject<icoModel>(turbulenceModel::propertiesName))
+            {
+                const icoModel& model = mesh_.lookupObject<icoModel>
+                (
+                    turbulenceModel::propertiesName
+                );
+                D = D + model.nut()*alphaDt_;
+            }
+
+            return D;
+        }
     }
 }
 
@@ -167,11 +169,45 @@ bool Foam::functionObjects::scalarTransport::read(const dictionary& dict)
 
     phiName_ = dict.lookupOrDefault<word>("phi", "phi");
     rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
+    muName_  = dict.lookupOrDefault<word>("mu", "thermo:mu");
+    nuName_  = dict.lookupOrDefault<word>("nu", "nu");
     schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
 
+    bool alphaDPresent = false;
+    if (dict.readIfPresent("alphaD", alphaD_))
+    {
+        alphaDPresent = true;
+    }
+
+    if (dict.readIfPresent("alphaDt", alphaDt_))
+    {
+        if (!alphaDPresent)
+        {
+            FatalErrorInFunction
+                << "Laminar alphaD must be specified in dictionary "
+                << dict.name() << exit(FatalError);
+        }
+    }
+    // Set turbulent alphaDt to 0 if it is not specified
+    else
+    {
+        alphaDt_ = 0.0;
+    }
+
     constantD_ = false;
     if (dict.readIfPresent("D", D_))
     {
+        if (alphaDPresent)
+        {
+            FatalErrorInFunction
+                << "Both D and alphaD are specified in dictionary "
+                << dict.name() << exit(FatalError);
+        }
+        constantD_ = true;
+    }
+    // If D and alphaD are not specified then D=0 is applied
+    else if (!alphaDPresent)
+    {
         constantD_ = true;
     }
 
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.H b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
index 86c92bd..7992419 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.H
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,8 +33,11 @@ Description
     - To specify the field name set the 'field' entry
     - To employ the same numerical schemes as another field set
       the 'schemesField' entry,
-    - The diffusivity can be set manually using the 'D' entry, or retrieved
-      from the turbulence model (if applicable).
+    - The constant diffusivity can be set manually using the 'D' entry
+    - The 'alphaD' and 'alphaDt' entries can be specified instead of 'D';
+      these are inverse Schmidt numbers (laminar and turbulent) defined as:
+      alphaD = 1/Sc, alphaDt = 1/Sct.
+      So diffusivity is evaluated as: D = nu*alphaD + nut*alphaDt
 
 See also
     Foam::functionObjects::fvMeshFunctionObject
@@ -77,12 +80,24 @@ class scalarTransport
         //- Name of density field (optional)
         word rhoName_;
 
+        //- Name of dynamic viscosity field (optional)
+        word muName_;
+
+        //- Name of kinematic viscosity field (optional)
+        word nuName_;
+
         //- Diffusion coefficient (optional)
         scalar D_;
 
         //- Flag to indicate whether a constant, uniform D_ is specified
         bool constantD_;
 
+        //- Inverse laminar Schmidt number alphaD = 1/Sc (optional)
+        scalar alphaD_;
+
+        //- Inverse turbulent Schmidt number alphaDt = 1/Sct (optional)
+        scalar alphaDt_;
+
         //- Number of corrector iterations (optional)
         label nCorr_;
 
scalarTransport_ver2.patch (6,103 bytes)   

henry

2017-02-10 14:07

manager   ~0007736

Thanks for the new patch, I will work on it over the weekend.

FatalError is functionObject does not terminate the program it terminates the functionObject construction which is then not added to the active functionObject list.

henry

2017-02-12 17:27

manager   ~0007751

I played with your proposed patch and found it rather complex, particularly in the way the laminar viscosities are handled. I have created a simpler implementation which I think covers your requirements while maintaining backward-compatibility:

commit 0714ccecd60957242c3f3a1d83e4ef6eb89077e6
Author: Henry Weller <http://cfd.direct>
Date: Sun Feb 12 17:19:27 2017 +0000

    functionObjects::scalarTransport: Added support for optional laminar and turbulent diffusion coefficients
    
    Description
        Evolves a passive scalar transport equation.
    
        - To specify the field name set the \c field entry
        - To employ the same numerical schemes as another field set
          the \c schemesField entry,
        - A constant diffusivity may be specified with the \c D entry,
    
        - Alternatively if a turbulence model is available a turbulent diffusivity
          may be constructed from the laminar and turbulent viscosities using the
          optional diffusivity coefficients \c alphaD and \c alphaDt (which default
          to 1):
          \verbatim
              D = alphaD*nu + alphaDt*nut
          \endverbatim
    
    Resolves feature request https://bugs.openfoam.org/view.php?id=2453

Let me know if there are any problems.

Sahas

2017-02-13 16:36

reporter   ~0007752

I have tested new functionality and it's ok, thank you!
Note that your implementation does not catch laminar case with constant Sc number and varying viscosity. But I think that this case is not too frequent.

henry

2017-02-13 17:05

manager   ~0007754

If you choose the laminar transport model in constant/turbulenceProperties:

simulationType laminar;

and set alphaD then it should handle the constant Sc number and varying viscosity case but I have not tested this. If this isn't working let me know.

Sahas

2017-02-13 17:24

reporter   ~0007756

I've found a misprint in scalarTransport.C, line 173:
alphaD_ = dict.lookupOrDefault("alphat", 1.0);
"alphat" -> "alphaD"

And indeed, setting "simulationType laminar" did the trick.

henry

2017-02-13 18:29

manager   ~0007757

Thanks for spotting the typo, I will push the fix shortly.

Issue History

Date Modified Username Field Change
2017-02-06 15:17 Sahas New Issue
2017-02-06 15:32 henry Note Added: 0007710
2017-02-06 16:05 Sahas Note Added: 0007712
2017-02-06 16:24 henry Note Added: 0007713
2017-02-08 11:00 henry Assigned To => henry
2017-02-08 11:00 henry Status new => closed
2017-02-08 11:00 henry Resolution open => fixed
2017-02-08 11:00 henry Note Added: 0007719
2017-02-08 14:31 Sahas Status closed => feedback
2017-02-08 14:31 Sahas Resolution fixed => reopened
2017-02-08 14:31 Sahas Note Added: 0007720
2017-02-08 14:52 henry Note Added: 0007721
2017-02-08 15:10 Sahas Note Added: 0007722
2017-02-08 15:10 Sahas Status feedback => assigned
2017-02-08 15:20 henry Note Added: 0007723
2017-02-08 16:07 Sahas Note Added: 0007724
2017-02-08 17:03 henry Note Added: 0007725
2017-02-08 17:30 Sahas File Added: scalarTransport.patch
2017-02-08 17:30 Sahas Note Added: 0007726
2017-02-08 17:39 henry Note Added: 0007727
2017-02-09 09:39 Sahas Note Added: 0007730
2017-02-09 09:49 henry Note Added: 0007731
2017-02-10 13:58 Sahas File Added: scalarTransport_ver2.patch
2017-02-10 13:58 Sahas Note Added: 0007735
2017-02-10 14:07 henry Note Added: 0007736
2017-02-12 17:27 henry Note Added: 0007751
2017-02-13 16:36 Sahas Note Added: 0007752
2017-02-13 17:05 henry Note Added: 0007754
2017-02-13 17:24 Sahas Note Added: 0007756
2017-02-13 18:29 henry Note Added: 0007757
2017-02-13 18:31 henry Status assigned => resolved
2017-02-13 18:31 henry Resolution reopened => fixed
2017-02-13 18:31 henry Fixed in Version => dev