Merge branch 'develop' - v7.1

krebeljk · Nov 9, 2020 · e688d60 · e688d60
2 parents 5914194 + a8cf8e7
commit e688d60
Show file tree

Hide file tree

Showing 48 changed files with 190 additions and 460 deletions.
diff --git a/README.md b/README.md
@@ -1,4 +1,4 @@
-# openInjMoldSim [`v7`](VERSION.md)
+# openInjMoldSim [`v7.1`](VERSION.md)
 
 This is an [OpenFOAM](https://en.wikipedia.org/wiki/OpenFOAM) solver for simulation of injection molding filling, packing and cooling stages.
 It is a modification of the compressibleInterFoam solver distributed with OpenFOAM.

diff --git a/VERSION.md b/VERSION.md
@@ -42,7 +42,6 @@ The naming of constant thermal conductivity as lambda was changed to kappa.
 - OpenFOAM 5
 
 ## Version v6
-
 - OpenFOAM 6
 
 ## Version v7
@@ -51,3 +50,10 @@ The naming of constant thermal conductivity as lambda was changed to kappa.
 - `pAux` auxiliary pressure field to balance the negative pressure during packing.
 - Calculation reordering to improve the restart behavior - work in progress.
 - Remove debug fields `p_rgh_resid` and `time`.
+
+## Version v7.1 (9. 11. 2020)
+- Resolve the boundary value viscosity issue that caused restart problems:
+  - Remove T.poly, T.air.
+  - Rename strig to shrRate.
+- Simplify the fiber code.
+- Improve demo case stability.
diff --git a/applications/solvers/multiphase/openInjMoldSim/TEqn.H b/applications/solvers/multiphase/openInjMoldSim/TEqn.H
@@ -4,7 +4,7 @@
         fvm::ddt(rho, T)
       + fvm::div(rhoPhi, T)
       - 1.0/mixture.Cp()*fvm::laplacian(mixture.kappa(), T)
-      - mixture.mu() * pow(strig,2)/mixture.Cp()//viscous heating
+      - mixture.mu() * pow(shrRate,2)/mixture.Cp()//viscous heating
       - T * mixture.cTExp()
         * (
             fvc::ddt(p)

diff --git a/applications/solvers/multiphase/openInjMoldSim/createFields.H b/applications/solvers/multiphase/openInjMoldSim/createFields.H
@@ -27,23 +27,11 @@ volVectorField U
     mesh
 );
 
-volScalarField shrRate
-(
-    IOobject
-    (
-        "shrRate",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    sqrt(2.0*symm(fvc::grad(U))&&symm(fvc::grad(U)))
-);
-
 #include "createPhi.H"
 
 Info<< "Constructing mojTwoPhaseMixtureThermo\n" << endl;
 mojTwoPhaseMixtureThermo mixture(mesh);
+mixture.correct();
 
 volScalarField& alpha1(mixture.alpha1());
 volScalarField& alpha2(mixture.alpha2());
@@ -52,7 +40,7 @@ Info<< "Reading thermophysical properties\n" << endl;
 
 volScalarField& p = mixture.p();
 volScalarField& T = mixture.T();
-volScalarField& strig = mixture.thermo1().strig();
+volScalarField& shrRate = mixture.shrRate();
 volScalarField& rho1 = mixture.thermo1().rho();
 const volScalarField& psi1 = mixture.thermo1().psi();
 volScalarField& rho2 = mixture.thermo2().rho();
@@ -105,28 +93,26 @@ volScalarField visc
         "visc",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    alpha1*mixture.thermo1().mu() + alpha2*mixture.thermo2().mu()
+    mixture.mu()
 );
 
 //kristjan - Calc Celsius
+dimensionedScalar twoSevenThree("twoSevenThree", dimensionSet(0,0,0,1,0,0,0), 273.15);
 volScalarField Tc
 (
     IOobject
     (
         "Tc",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
-    mesh,
-    dimensionedScalar("initTc", dimensionSet(0,0,0,1,0,0,0), 0.0)
+    T - twoSevenThree
 );
-dimensionedScalar twoSevenThree("twoSevenThree", dimensionSet(0,0,0,1,0,0,0), 273.15);
-Tc = T - twoSevenThree;
 
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransport.H b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransport.H
@@ -174,7 +174,7 @@ public:
         }
 
         //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar p, const scalar T, const scalar strig) const;
+        inline scalar mu(const scalar p, const scalar T, const scalar shrRate) const;
 
         //- Thermal conductivity [W/mK]
         inline scalar kappa(const scalar p, const scalar T) const;

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportI.H b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportI.H
@@ -163,15 +163,15 @@ inline Foam::scalar Foam::crossWLFTransport<Thermo>::mu
 (
     const scalar p,
     const scalar T,
-    const scalar strig
+    const scalar shrRate
 ) const
 {   
     if (T < TnoFlow_) { return etaMax_; } //leave immediately to avoid n0 = infinite
     const scalar n0 = D1_*::exp(-(A1_*(T-D2_-D3_*p))/(A2_+T-D2_));
-    scalar visk = n0/(1.0+pow(n0*strig/Tau_,1.0-n_));
+    scalar visk = n0/(1.0+pow(n0*shrRate/Tau_,1.0-n_));
     scalar Tstart = TnoFlow_ + deltaTempInterp_; //solidification onset temperature
 
-    //Info << "visk= " << visk  << " p=" << p << " T=" << T << " strig=" << strig << endl;   // debug Kristjan
+    //Info << "visk= " << visk  << " p=" << p << " T=" << T << " shrRate=" << shrRate << endl;   // debug Kristjan
 
     if (visk > etaMax_) { visk = etaMax_; }
     if (visk < etaMin_) { visk = etaMin_; }

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportTTC.H b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportTTC.H
@@ -177,7 +177,7 @@ public:
         }
 
         //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar p, const scalar T, const scalar strig) const;
+        inline scalar mu(const scalar p, const scalar T, const scalar shrRate) const;
 
         //- Thermal conductivity [W/mK]
         inline scalar kappa(const scalar p, const scalar T) const;

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportTTCI.H b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/crossWLFTransportTTCI.H
@@ -167,15 +167,15 @@ inline Foam::scalar Foam::crossWLFTransportTTC<Thermo>::mu
 (
     const scalar p,
     const scalar T,
-    const scalar strig
+    const scalar shrRate
 ) const
 {   
     if (T < TnoFlow_) { return etaMax_; } //leave immediately to avoid n0 = infinite
     const scalar n0 = D1_*::exp(-(A1_*(T-D2_-D3_*p))/(A2_+T-D2_));
-    scalar visk = n0/(1.0+pow(n0*strig/Tau_,1.0-n_));
+    scalar visk = n0/(1.0+pow(n0*shrRate/Tau_,1.0-n_));
     scalar Tstart = TnoFlow_ + deltaTempInterp_; //solidification onset temperature
 
-    //Info << "visk= " << visk  << " p=" << p << " T=" << T << " strig=" << strig << endl;   // debug Kristjan
+    //Info << "visk= " << visk  << " p=" << p << " T=" << T << " shrRate=" << shrRate << endl;   // debug Kristjan
 
     if (visk > etaMax_) { visk = etaMax_; }
     if (visk < etaMin_) { visk = etaMin_; }

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojBasicThermo.C b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojBasicThermo.C
@@ -178,20 +178,9 @@ Foam::mojBasicThermo::mojBasicThermo
 
     p_(lookupOrConstruct(mesh, "p")),
 
-    T_
-    (
-        IOobject
-        (
-            phasePropertyName("T"),
-            mesh.time().timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    ),
+    T_(lookupOrConstruct(mesh, "T")),
 
-    strig_(lookupOrConstruct(mesh, "strig")),
+    shrRate_(lookupOrConstruct(mesh, "shrRate")),
 
     alpha_
     (
@@ -235,20 +224,10 @@ Foam::mojBasicThermo::mojBasicThermo
 
     p_(lookupOrConstruct(mesh, "p")),
 
-    T_
-    (
-        IOobject
-        (
-            phasePropertyName("T"),
-            mesh.time().timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    ),
 
-    strig_(lookupOrConstruct(mesh, "strig")),
+    T_(lookupOrConstruct(mesh, "T")),
+
+    shrRate_(lookupOrConstruct(mesh, "shrRate")),
 
     alpha_
     (
@@ -490,14 +469,14 @@ const Foam::volScalarField& Foam::mojBasicThermo::p() const
 }
 
 
-Foam::volScalarField& Foam::mojBasicThermo::strig()
+Foam::volScalarField& Foam::mojBasicThermo::shrRate()
 {
-    return strig_;
+    return shrRate_;
 }
 
-const Foam::volScalarField& Foam::mojBasicThermo::strig() const
+const Foam::volScalarField& Foam::mojBasicThermo::shrRate() const
 {
-    return strig_;
+    return shrRate_;
 }
 
 const Foam::volScalarField& Foam::mojBasicThermo::T() const

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojBasicThermo.H b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojBasicThermo.H
@@ -69,10 +69,10 @@ protected:
             volScalarField& p_;
 
             //- Temperature [K]
-            volScalarField T_;
+            volScalarField& T_;
 
             //- Temperature [1/s]
-            volScalarField strig_;
+            volScalarField& shrRate_;
 
             //- Laminar thermal diffusuvity [kg/m/s]
             volScalarField alpha_;
@@ -279,10 +279,10 @@ public:
             virtual const volScalarField& p() const;
 
 	    //- Shear rate [1/s]
-	    virtual volScalarField& strig();
+	    virtual volScalarField& shrRate();
 
 	    //- Shear rate [1/s]
-	    virtual const volScalarField& strig() const;
+	    virtual const volScalarField& shrRate() const;
 
 
             //- Density [kg/m^3]

diff --git a/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojHeRhoThermo.C b/applications/solvers/multiphase/openInjMoldSim/crossWlfInSorodno/mojHeRhoThermo.C
@@ -32,7 +32,7 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
 {
     //const scalarField& hCells = this->he(); //Kristjan: governed by TEqn
     const scalarField& pCells = this->p_;
-    const scalarField& strigCells = this->strig_;
+    const scalarField& shrRateCells = this->shrRate_;
 
     scalarField& TCells = this->T_.primitiveFieldRef();
     scalarField& psiCells = this->psi_.primitiveFieldRef();
@@ -55,7 +55,7 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
         psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
         rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);
 
-        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli], strigCells[celli]);
+        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli], shrRateCells[celli]);
 
         alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
     }
@@ -66,8 +66,8 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
     volScalarField::Boundary& TBf =
         this->T_.boundaryFieldRef();
 
-    volScalarField::Boundary& strigBf =
-        this->strig_.boundaryFieldRef();
+    volScalarField::Boundary& shrRateBf =
+        this->shrRate_.boundaryFieldRef();
 
     volScalarField::Boundary& psiBf =
         this->psi_.boundaryFieldRef();
@@ -85,8 +85,8 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
     {
         fvPatchScalarField& pp = pBf[patchi];
         fvPatchScalarField& pT = TBf[patchi];
-	fvPatchScalarField& pstrig = strigBf[patchi];
-	//fvPatchScalarField& pU = strig.boundaryField()[patchi]; // tukaj boundary
+	fvPatchScalarField& pshrRate = shrRateBf[patchi];
+	//fvPatchScalarField& pU = shrRate.boundaryField()[patchi]; // tukaj boundary
 
         fvPatchScalarField& ppsi = psiBf[patchi];
         fvPatchScalarField& prho = rhoBf[patchi];
@@ -109,7 +109,7 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
 
                 ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                 prho[facei] = mixture_.rho(pp[facei], pT[facei]);
-                pmu[facei] = mixture_.mu(pp[facei], pT[facei], pstrig[facei]);
+                pmu[facei] = mixture_.mu(pp[facei], pT[facei], pshrRate[facei]);
                 palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);
             }
         }
@@ -124,7 +124,7 @@ void Foam::mojHeRhoThermo<BasicPsiThermo, MixtureType>::calculate()
 
                 ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                 prho[facei] = mixture_.rho(pp[facei], pT[facei]);
-                pmu[facei] = mixture_.mu(pp[facei], pT[facei], pstrig[facei]);
+                pmu[facei] = mixture_.mu(pp[facei], pT[facei], pshrRate[facei]);
                 palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);
             }
         }

diff --git a/applications/solvers/multiphase/openInjMoldSim/mojConstTransport/mojConstTransport.H b/applications/solvers/multiphase/openInjMoldSim/mojConstTransport/mojConstTransport.H
@@ -144,7 +144,7 @@ public:
         }
 
         //- Dynamic viscosity [kg/ms]
-        inline scalar mu(const scalar p, const scalar T, const scalar strig) const;
+        inline scalar mu(const scalar p, const scalar T, const scalar shrRate) const;
 
         //- Thermal conductivity [W/mK]
         inline scalar kappa(const scalar p, const scalar T) const;

diff --git a/applications/solvers/multiphase/openInjMoldSim/mojConstTransport/mojConstTransportI.H b/applications/solvers/multiphase/openInjMoldSim/mojConstTransport/mojConstTransportI.H
@@ -98,7 +98,7 @@ inline Foam::scalar Foam::mojConstTransport<Thermo>::mu
 (
     const scalar p,
     const scalar T,
-    const scalar strig
+    const scalar shrRate
 ) const
 {
     return mu_;

diff --git a/...ons/solvers/multiphase/openInjMoldSim/mojTwoPhaseMixtureThermo/mojTwoPhaseMixtureThermo.C b/...ons/solvers/multiphase/openInjMoldSim/mojTwoPhaseMixtureThermo/mojTwoPhaseMixtureThermo.C
@@ -48,16 +48,6 @@ Foam::mojTwoPhaseMixtureThermo::mojTwoPhaseMixtureThermo
     thermo1_(NULL),
     thermo2_(NULL)
 {
-    {
-        volScalarField T1(IOobject::groupName("T", phase1Name()), T_);
-        T1.write();
-    }
-
-    {
-        volScalarField T2(IOobject::groupName("T", phase2Name()), T_);
-        T2.write();
-    }
-
     thermo1_ = mojRhoThermo::New(mesh, phase1Name());
     thermo2_ = mojRhoThermo::New(mesh, phase2Name());
 
@@ -78,10 +68,8 @@ Foam::mojTwoPhaseMixtureThermo::~mojTwoPhaseMixtureThermo()
 
 void Foam::mojTwoPhaseMixtureThermo::correct()
 {
-    thermo1_->T() = T_;
     thermo1_->correct();
 
-    thermo2_->T() = T_;
     thermo2_->correct();
 
     psi_ = alpha1()*thermo1_->psi() + alpha2()*thermo2_->psi();

diff --git a/applications/solvers/multiphase/openInjMoldSim/openInjMoldSim.C b/applications/solvers/multiphase/openInjMoldSim/openInjMoldSim.C
@@ -50,7 +50,7 @@ Description
 
 int main(int argc, char *argv[])
 {
-    Info << "openInjMoldSim v7" << endl;
+    Info << "openInjMoldSim v7.1" << endl;
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -91,10 +91,9 @@ int main(int argc, char *argv[])
             solve(fvm::ddt(rho) + fvc::div(rhoPhi));
 
             //update state
-            strig = sqrt(2.0*symm(fvc::grad(U))&&symm(fvc::grad(U)));
-            shrRate = strig;
+            shrRate = sqrt(2.0*symm(fvc::grad(U))&&symm(fvc::grad(U)));
             mixture.correct();
-            visc = alpha1*mixture.thermo1().mu() + alpha2*mixture.thermo2().mu();
+            visc = mixture.mu();
             mojKappaOut = mixture.kappa();
 
             //Kristjan: Elastic deviatoric stress equation