diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C index f0feddc97a698826958f3568a1a96bbc7906c9a9..06e6826fc852566ec3c122383003b034e8582c47 100644 --- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C @@ -160,10 +160,27 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL; // Momentum source due to particle forces - const vector FCoupled = - mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U); - const vector FNonCoupled = - mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U); + const vector FCoupled = mass*td.cloud().forces().calcCoupled + ( + cellI, + dt, + rhoc_, + rho, + Uc_, + U, + d + ); + + const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled + ( + cellI, + dt, + rhoc_, + rho, + Uc_, + U, + d + ); // New particle velocity diff --git a/src/lagrangian/intermediate/particleForces/particleForces.C b/src/lagrangian/intermediate/particleForces/particleForces.C index 7a07051def2cc55c69074e48a3758d99b3fcc8de..4f8fe1cfc7fb95196e11c22e496149b78302d42b 100644 --- a/src/lagrangian/intermediate/particleForces/particleForces.C +++ b/src/lagrangian/intermediate/particleForces/particleForces.C @@ -27,6 +27,8 @@ License #include "fvMesh.H" #include "volFields.H" #include "fvcGrad.H" +#include "mathematicalConstants.H" +#include "electromagneticConstants.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -41,16 +43,25 @@ Foam::particleForces::particleForces dict_(dict.subDict("particleForces")), g_(g), gradUPtr_(NULL), + gradHPtr_(NULL), gravity_(dict_.lookup("gravity")), virtualMass_(dict_.lookup("virtualMass")), Cvm_(0.0), pressureGradient_(dict_.lookup("pressureGradient")), - UName_(dict_.lookupOrDefault<word>("U", "U")) + paramagnetic_(dict_.lookup("paramagnetic")), + chi_(0.0), + UName_(dict_.lookupOrDefault<word>("U", "U")), + HName_(dict_.lookupOrDefault<word>("H", "H")) { if (virtualMass_) { dict_.lookup("Cvm") >> Cvm_; } + + if (paramagnetic_) + { + dict_.lookup("chi") >> chi_; + } } @@ -60,11 +71,15 @@ Foam::particleForces::particleForces(const particleForces& f) dict_(f.dict_), g_(f.g_), gradUPtr_(f.gradUPtr_), + gradHPtr_(f.gradHPtr_), gravity_(f.gravity_), virtualMass_(f.virtualMass_), Cvm_(f.Cvm_), pressureGradient_(f.pressureGradient_), - UName_(f.UName_) + paramagnetic_(f.paramagnetic_), + chi_(f.chi_), + UName_(f.UName_), + HName_(f.HName_) {} @@ -102,18 +117,42 @@ Foam::Switch Foam::particleForces::virtualMass() const } +Foam::scalar Foam::particleForces::Cvm() const +{ + return Cvm_; +} + + Foam::Switch Foam::particleForces::pressureGradient() const { return pressureGradient_; } +Foam::Switch Foam::particleForces::paramagnetic() const +{ + return paramagnetic_; +} + + +Foam::scalar Foam::particleForces::chi() const +{ + return chi_; +} + + const Foam::word& Foam::particleForces::UName() const { return UName_; } +const Foam::word& Foam::particleForces::HName() const +{ + return HName_; +} + + void Foam::particleForces::cacheFields(const bool store) { if (store && pressureGradient_) @@ -129,6 +168,20 @@ void Foam::particleForces::cacheFields(const bool store) gradUPtr_ = NULL; } } + + if (store && paramagnetic_) + { + const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_); + gradHPtr_ = fvc::grad(H).ptr(); + } + else + { + if (gradHPtr_) + { + delete gradHPtr_; + gradHPtr_ = NULL; + } + } } @@ -139,7 +192,8 @@ Foam::vector Foam::particleForces::calcCoupled const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const { vector Ftot = vector::zero; @@ -172,7 +226,8 @@ Foam::vector Foam::particleForces::calcNonCoupled const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const { vector Ftot = vector::zero; @@ -183,9 +238,34 @@ Foam::vector Foam::particleForces::calcNonCoupled Ftot += g_*(1.0 - rhoc/rho); } + // Magnetic field force + + if (paramagnetic_) + { + const volVectorField& H = mesh_.lookupObject<volVectorField>(HName_); + + const volTensorField& gradH = *gradHPtr_; + + Ftot += + 3.0*constant::electromagnetic::mu0.value()/rho + *chi_/(chi_ + 3) + *(H[cellI] & gradH[cellI]); + + // force is: + + // 4.0 + // *constant::mathematical::pi + // *constant::electromagnetic::mu0.value() + // *pow3(d/2) + // *chi/(chi + 3) + // *(H[cellI] & gradH[cellI]); + + // which is divided by mass ((4/3)*pi*r^3*rho) to produce + // acceleration + } + return Ftot; } // ************************************************************************* // - diff --git a/src/lagrangian/intermediate/particleForces/particleForces.H b/src/lagrangian/intermediate/particleForces/particleForces.H index d6d988c4f53f338a58e5af511c5cf506bcf9e09d..22299e34bd928f177bf60787c0e53384358f3cbd 100644 --- a/src/lagrangian/intermediate/particleForces/particleForces.H +++ b/src/lagrangian/intermediate/particleForces/particleForces.H @@ -69,6 +69,9 @@ class particleForces //- Velocity gradient field const volTensorField* gradUPtr_; + //- Magnetic field strength gradient field + const volTensorField* gradHPtr_; + // Forces to include in particle motion evaluation @@ -84,12 +87,21 @@ class particleForces //- Pressure gradient Switch pressureGradient_; + //- Paramagnetic force + Switch paramagnetic_; + + //- Magnetic susceptibility of particle + scalar chi_; + // Additional info - //- Name of velucity field - default = "U" + //- Name of velocity field - default = "U" const word UName_; + //- Name of magnetic field strength field - default = "H" + const word HName_; + public: @@ -128,7 +140,7 @@ public: Switch virtualMass() const; //- Return virtual mass force coefficient - Switch Cvm() const; + scalar Cvm() const; //- Return pressure gradient force activate switch Switch pressureGradient() const; @@ -136,6 +148,15 @@ public: //- Return name of velocity field const word& UName() const; + //- Return paramagnetic force activate switch + Switch paramagnetic() const; + + //- Return magnetic susceptibility + scalar chi() const; + + //- Return name of velocity field + const word& HName() const; + // Evaluation @@ -150,7 +171,8 @@ public: const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const; //- Calculate external forces applied to the particles @@ -161,7 +183,8 @@ public: const scalar rhoc, const scalar rho, const vector& Uc, - const vector& U + const vector& U, + const scalar d ) const; };