ENH: Final stages for mass transfer models and laser modification for interface reflection

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/TEqn.H

@@ -10,8 +10,6 @@
         fluid.kappa() + fluid.Cp()*turbulence->mut()/fluid.Prt()
     );
 
-    //dimensionedScalar S("S", dimEnergy/dimVolume/dimTime, 1.225e8);
-
     fvScalarMatrix TEqn
     (
         fvm::ddt(rhoCp, T)
@@ -24,6 +22,11 @@
       + fvOptions(rhoCp, T)
     );
 
+    if (T.mesh().time().outputTime())
+    {
+        kappaEff.write();
+    }
+
     TEqn.relax();
 
     fvOptions.constrain(TEqn);

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/icoReactingMultiphaseInterFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,8 +106,6 @@ int main(int argc, char *argv[])
             }
         }
 
-
-
         rho = fluid.rho();
 
         runTime.write();

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C

@@ -38,8 +38,7 @@ Foam::DTRMParticle::DTRMParticle
     const scalar I,
     const label cellI,
     const scalar dA,
-    const label reflectedId,
-    const scalar Imin,
+    const label transmissiveId,
     bool doCellFacePt
 )
 :
@@ -49,8 +48,7 @@ Foam::DTRMParticle::DTRMParticle
     I0_(I),
     I_(I),
     dA_(dA),
-    reflectedId_(reflectedId),
-    Imin_(Imin)
+    transmissiveId_(transmissiveId)
 {}
 
 
@@ -62,8 +60,7 @@ Foam::DTRMParticle::DTRMParticle(const DTRMParticle& p)
     I0_(p.I0_),
     I_(p.I_),
     dA_(p.dA_),
-    reflectedId_(p.reflectedId_),
-    Imin_(p.Imin_)
+    transmissiveId_(p.transmissiveId_)
 {}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -94,43 +91,28 @@ bool Foam::DTRMParticle::move
 
         // Boltzman constant
         const scalar sigma = physicoChemical::sigma.value();
+
+        label reflectedZoneId = td.relfectedCells()[celli];
+
         if
         (
-            (!td.relfectedCells()[celli] > 0 && reflectedId_ == 0)
-         || reflectedId_ > 0
+            (reflectedZoneId > -1)
+         && (
+                (transmissiveId_ == -1)
+             || (transmissiveId_ != reflectedZoneId)
+            )
         )
-        {
-            scalar a = td.aInterp().interpolate(position(), cell0);
-            scalar e = td.eInterp().interpolate(position(), cell0);
-            scalar E = td.EInterp().interpolate(position(), cell0);
-            scalar T = td.TInterp().interpolate(position(), cell0);
-
-            const scalar I1 =
-            (
-                I_
-                + ds*(e*sigma*pow4(T)/mathematical::pi + E)
-            ) / (1 + ds*a);
-
-            td.Q(cell0) += (I_ - I1)*dA_;
-
-            I_ = I1;
-
-            if ((I_ <= 0.01*I0_) || (I_ < Imin_))
-            {
-                break;
-            }
-        }
-        else
         {
             scalar rho(0);
+
             // Create a new reflected particle when the particles is not
             // transmissive and larger than an absolute I
-            if (reflectedId_ == 0 && I_ > Imin_)
+            if (I_ > 0.01*I0_)
             {
                 vector pDir = dsv/ds;
 
                 cellPointWeight cpw(mesh_, position(), celli, face());
-                //vector nHat = td.nHatCells()[celli];
+
                 vector nHat = td.nHatInterp().interpolate(cpw);
 
                 nHat /= mag(nHat);
@@ -138,12 +120,26 @@ bool Foam::DTRMParticle::move
                 // Only new incoming rays
                 if (cosTheta > SMALL)
                 {
-                    vector newDir = td.reflection().R(pDir, nHat);
-
-                    //scalar theta = acos(-pDir & nHat);
+                    vector newDir =
+                        td.reflection()
+                        [
+                            td.relfectedCells()[celli]
+                        ].R(pDir, nHat);
 
                     // reflectivity
-                    rho = min(max(td.reflection().rho(cosTheta), 0.0), 0.98);
+                    rho =
+                        min
+                        (
+                            max
+                            (
+                                td.reflection()
+                                [
+                                    td.relfectedCells()[celli]
+                                ].rho(cosTheta)
+                                , 0.0
+                            )
+                            , 0.98
+                        );
 
                     scalar delaM = sqrt(mesh_.cellVolumes()[cell0]);
 
@@ -155,8 +151,7 @@ bool Foam::DTRMParticle::move
                         I_*rho,
                         cell0,
                         dA_,
-                        reflectedId_,
-                        Imin_,
+                        transmissiveId_,
                         true
                     );
                     // Add to cloud
@@ -164,7 +159,8 @@ bool Foam::DTRMParticle::move
                 }
             }
 
-            reflectedId_++;
+            // Change transmissiveId of the particle
+            transmissiveId_ = reflectedZoneId;
 
             const point p0 = position();
 
@@ -186,11 +182,33 @@ bool Foam::DTRMParticle::move
               + ds*(e*sigma*pow4(T)/mathematical::pi + E)
             ) / (1 + ds*a);
 
-            td.Q(celli) += (Itran - I1)*dA_;
+            td.Q(celli) += (Itran - max(I1, 0.0))*dA_;
+
+            I_ = I1;
+
+            if (I_ <= 0.01*I0_)
+            {
+                break;
+            }
+        }
+        else
+        {
+            scalar a = td.aInterp().interpolate(position(), cell0);
+            scalar e = td.eInterp().interpolate(position(), cell0);
+            scalar E = td.EInterp().interpolate(position(), cell0);
+            scalar T = td.TInterp().interpolate(position(), cell0);
+
+            const scalar I1 =
+            (
+                I_
+                + ds*(e*sigma*pow4(T)/mathematical::pi + E)
+            ) / (1 + ds*a);
+
+            td.Q(cell0) += (I_ -  max(I1, 0.0))*dA_;
 
             I_ = I1;
 
-            if (I_ <= 0.01*I0_ || I_ < Imin_)
+            if ((I_ <= 0.01*I0_))
             {
                 break;
             }

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H

@@ -43,7 +43,6 @@ SourceFiles
 #include "reflectionModel.H"
 #include "interpolationCellPoint.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -65,6 +64,7 @@ class DTRMParticle
 {
 private:
 
+
     // Private data
 
         //- Size in bytes of the fields
@@ -85,11 +85,8 @@ private:
         //- Area of radiation
         scalar dA_;
 
-        //- Reflected index
-        label reflectedId_;
-
-        //- Minimum radiation intensity to which the particle is tracked [W/m2]
-        scalar Imin_;
+        //- Trasnmissive index
+        label transmissiveId_;
 
 
 public:
@@ -101,6 +98,8 @@ public:
     :
         public particle::TrackingData<Cloud<DTRMParticle>>
     {
+
+
         // Interpolators for continuous phase fields
 
             const interpolationCell<scalar>& aInterp_;
@@ -111,7 +110,8 @@ public:
             const interpolationCellPoint<vector>& nHatInterp_;
 
             const labelField& relfectedCells_;
-            const reflectionModel& reflection_;
+            //const reflectionModelTable& reflection_;
+            UPtrList<reflectionModel> reflection_;
 
             //- Heat source term
             volScalarField& Q_;
@@ -130,7 +130,7 @@ public:
                 const interpolationCell<scalar>& TInterp,
                 const interpolationCellPoint<vector>& nHatInterp,
                 const labelField&,
-                const reflectionModel&,
+                const UPtrList<reflectionModel>&,
                 volScalarField& Q
             );
 
@@ -147,7 +147,7 @@ public:
             inline const interpolationCell<scalar>& TInterp() const;
             inline const interpolationCellPoint<vector>& nHatInterp() const;
             inline const labelField& relfectedCells() const;
-            inline const reflectionModel& reflection() const;
+            inline const UPtrList<reflectionModel>& reflection() const;
 
             inline scalar& Q(label celli);
 
@@ -190,8 +190,7 @@ public:
             const scalar I,
             const label cellI,
             const scalar dA,
-            const label reflectedId,
-            const scalar Imin,
+            const label transmissiveId,
             bool doCellFacePt = true
         );
 
@@ -274,15 +273,6 @@ public:
         //- Return access to reflectedId
         inline label& reflectedId();
 
-        //- Return access to initial tet face
-        //inline label& tetFace0();
-
-        //- Return access to initial tet point
-        //inline label& tetPt0();
-
-        //- Return access to initial proc Id
-        //inline label& origProc0();
-
 
    // Tracking
 

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H

@@ -34,7 +34,7 @@ inline Foam::DTRMParticle::trackingData::trackingData
     const interpolationCell<scalar>& TInterp,
     const interpolationCellPoint<vector>& nHatInterp,
     const labelField& relfectedCell,
-    const reflectionModel& reflection,
+    const UPtrList<reflectionModel>& reflection,
     volScalarField& Q
 )
 :
@@ -92,7 +92,7 @@ Foam::DTRMParticle::trackingData::relfectedCells() const
 }
 
 
-inline const Foam::reflectionModel&
+inline const Foam::UPtrList<Foam::radiation::reflectionModel>&
 Foam::DTRMParticle::trackingData::reflection() const
 {
     return reflection_;
@@ -147,12 +147,6 @@ inline Foam::point& Foam::DTRMParticle::p0()
 }
 
 
-inline Foam::label& Foam::DTRMParticle::reflectedId()
-{
-    return reflectedId_;
-}
-
-
 inline Foam::point& Foam::DTRMParticle::p1()
 {
     return p1_;

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/Make/options

 EXE_INC = \
-    -DFULLDEBUG -g -O0 \
+    -I../phasesSystem/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \

applications/solvers/multiphase/icoReactingMultiPhaseInterFoam/laserDTRM/laserDTRM.C

@@ -124,19 +124,13 @@ Foam::tmp<Foam::volVectorField> Foam::radiation::laserDTRM::nHatfv
     const dimensionedScalar deltaN
     (
         "deltaN",
-        1e-8/pow(average(mesh_.V()), 1.0/3.0)
-    );
-
-    const dimensionedScalar minAlpha
-    (
-        "minAlpha", dimless, 1e-3
+        1e-7/pow(average(mesh_.V()), 1.0/3.0)
     );
 
     volVectorField gradAlphaf
     (
-        "gradAlphaf",
-         (alpha2 + minAlpha)*fvc::grad(alpha1)
-       - (alpha1 + minAlpha)*fvc::grad(alpha2)
+         alpha2*fvc::grad(alpha1)
+       - alpha1*fvc::grad(alpha2)
     );
 
    // Face unit interface normal
@@ -156,6 +150,48 @@ Foam::tmp<Foam::volScalarField>Foam::radiation::laserDTRM::nearInterface
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::radiation::laserDTRM::sPhase
+(
+    const volScalarField& alpha1,
+    const volScalarField& alpha2
+) const
+{
+   // positive : the phases have similar orientation other
+   // negative: the phases oppose each
+   return  fvc::grad(alpha1) & fvc::grad(alpha2);
+}
+
+
+
+void Foam::radiation::laserDTRM::initialiseReflection()
+{
+    if (found("reflectionModel"))
+    {
+        dictTable modelDicts(lookup("reflectionModel"));
+
+        forAllConstIter(dictTable, modelDicts, iter)
+        {
+            const phasePairKey& key = iter.key();
+
+            reflections_.insert
+            (
+                key,
+                reflectionModel::New
+                (
+                    *iter,
+                    mesh_
+                )
+            );
+        }
+
+        if (reflections_.size() > 0)
+        {
+            reflectionSwitch_ = true;
+        }
+    }
+}
+
+
 void Foam::radiation::laserDTRM::initialise()
 {
     // Initialise the DTRM particles
@@ -203,19 +239,6 @@ void Foam::radiation::laserDTRM::initialise()
                 new interpolation2DTable<scalar>(*this)
             );
 
-            // Check dimensions ndr and ndTheta
-//             if
-//             (
-//                (powerDistribution_->size() != ndTheta_)
-//             || (powerDistribution_().first().second().size() != ndr_)
-//             )
-//             {
-//                  FatalErrorInFunction
-//                     << " The table dimensions should correspond with ndTheta "
-//                     << " and ndr "
-//                     << exit(FatalError);
-//             }
-
             break;
         }
         case pdUniform:
@@ -235,9 +258,6 @@ void Foam::radiation::laserDTRM::initialise()
     p1 = p0;
     if (mesh_.nGeometricD() == 3)
     {
-        //scalar r0 = dr/2.0;
-        //scalar r1Max0 = drMax/2.0;
-
         for (label ri = 0; ri < ndr_; ri++)
         {
             scalar r1 = SMALL + dr*ri;
@@ -275,8 +295,6 @@ void Foam::radiation::laserDTRM::initialise()
                 // calculate target point using new deviation rl
                 p1 = lPosition + finalPos + (0.5*maxTrackLength_*lDir);
 
-                //scalar p = magSqr(p0 - lPosition);
-
                 scalar Ip = calculateIp(rP, thetaP);
 
                 scalar dAi = (sqr(r2) - sqr(r1))*(theta2 - theta1)/2.0;
@@ -290,7 +308,7 @@ void Foam::radiation::laserDTRM::initialise()
                 {
                     // Create a new particle
                     DTRMParticle* pPtr = new DTRMParticle
-                        (mesh_, p0, p1, Ip, cellI, dAi, 0 , 0.01*Ip, true);
+                        (mesh_, p0, p1, Ip, cellI, dAi, -1, true);
 
                     // Add to cloud
                     DTRMCloud_.addParticle(pPtr);
@@ -338,12 +356,12 @@ Foam::radiation::laserDTRM::laserDTRM(const volScalarField& T)
     (
         Function1<point>::New("focalLaserPosition", *this)
     ),
+
     laserDirection_
     (
         Function1<vector>::New("laserDirection", *this)
     ),
 
-
     focalLaserRadius_(readScalar(lookup("focalLaserRadius"))),
     qualityBeamLaser_
     (
@@ -354,11 +372,10 @@ Foam::radiation::laserDTRM::laserDTRM(const volScalarField& T)
     laserPower_(Function1<scalar>::New("laserPower", *this)),
     powerDistribution_(),
 
-    reflection_(reflectionModel::New(*this, mesh_).ptr()),
-    reflectionSwitch_(lookupOrDefault("reflection", false)),
-    initialPhase_(lookupOrDefault("initialPhase", word::null)),
-    alpha1_(lookupOrDefault("alpha1", word::null)),
-    alpha2_(lookupOrDefault("alpha2", word::null)),
+    reflectionSwitch_(false),
+
+    alphaCut_( lookupOrDefault<scalar>("alphaCut", 0.5)),
+
     Qin_
     (
         IOobject
@@ -413,7 +430,7 @@ Foam::radiation::laserDTRM::laserDTRM(const volScalarField& T)
     ),
     Q_
     (
-     IOobject
+        IOobject
         (
             "Q",
             mesh_.time().timeName(),
@@ -425,6 +442,8 @@ Foam::radiation::laserDTRM::laserDTRM(const volScalarField& T)
         dimensionedScalar("Q", dimPower/dimVolume, 0.0)
     )
 {
+    initialiseReflection();
+
     initialise();
 }
 
@@ -462,11 +481,10 @@ Foam::radiation::laserDTRM::laserDTRM
     laserPower_(Function1<scalar>::New("laserPower", *this)),
     powerDistribution_(),
 
-    reflection_(reflectionModel::New(*this, mesh_).ptr()),
-    reflectionSwitch_(dict.lookupOrDefault("reflection", false)),
-    initialPhase_(lookupOrDefault("initialPhase", word::null)),
-    alpha1_(lookupOrDefault("alpha1", word::null)),
-    alpha2_(lookupOrDefault("alpha2", word::null)),
+    reflectionSwitch_(false),
+
+    alphaCut_( lookupOrDefault<scalar>("alphaCut", 0.5)),
+
     Qin_
     (
         IOobject
@@ -533,6 +551,7 @@ Foam::radiation::laserDTRM::laserDTRM
         dimensionedScalar("Q", dimPower/pow3(dimLength), 0.0)
     )
 {
+    initialiseReflection();
     initialise();
 }
 
@@ -573,25 +592,11 @@ void Foam::radiation::laserDTRM::calculate()
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedScalar("zero", dimless, 0)
+            dimensionedScalar("zero", dimless, -1)
         )
     );
     volScalarField& reflectingCellsVol = treflectingCells.ref();
 
-    // Reset the fields
-    Qin_ == dimensionedScalar("zero", Qin_.dimensions(), 0);
-    Q_ == dimensionedScalar("zero", Q_.dimensions(), 0);
-
-    a_ = absorptionEmission_->a();
-    e_ = absorptionEmission_->e();
-    E_ = absorptionEmission_->E();
-
-    const interpolationCell<scalar> aInterp(a_);
-    const interpolationCell<scalar> eInterp(e_);
-    const interpolationCell<scalar> EInterp(E_);
-    const interpolationCell<scalar> TInterp(T_);
-
-    labelField reflectingCells(mesh_.nCells(), 0);
 
     tmp<volVectorField> tnHat
     (
@@ -611,54 +616,91 @@ void Foam::radiation::laserDTRM::calculate()
     );
     volVectorField& nHat = tnHat.ref();
 
+    // Reset the fields
+    Qin_ == dimensionedScalar("zero", Qin_.dimensions(), 0);
+    Q_ == dimensionedScalar("zero", Q_.dimensions(), 0);
+
+    a_ = absorptionEmission_->a();
+    e_ = absorptionEmission_->e();
+    E_ = absorptionEmission_->E();
+
+    const interpolationCell<scalar> aInterp(a_);
+    const interpolationCell<scalar> eInterp(e_);
+    const interpolationCell<scalar> EInterp(E_);
+    const interpolationCell<scalar> TInterp(T_);
+
+    labelField reflectingCells(mesh_.nCells(), -1);
+
     autoPtr<interpolationCellPoint<vector>> nHatIntrPtr;
 
+    UPtrList<reflectionModel> reflactionUPtr;
+
     if (reflectionSwitch_)