Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "ETAB.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/mathematicalConstants.H>
00029
00030
00031
00032 namespace Foam
00033 {
00034
00035
00036
00037 defineTypeNameAndDebug(ETAB, 0);
00038
00039 addToRunTimeSelectionTable
00040 (
00041 breakupModel,
00042 ETAB,
00043 dictionary
00044 );
00045
00046
00047
00048
00049 ETAB::ETAB
00050 (
00051 const dictionary& dict,
00052 spray& sm
00053 )
00054 :
00055 breakupModel(dict, sm),
00056 coeffsDict_(dict.subDict(typeName + "Coeffs")),
00057 Cmu_(readScalar(coeffsDict_.lookup("Cmu"))),
00058 Comega_(readScalar(coeffsDict_.lookup("Comega"))),
00059 k1_(readScalar(coeffsDict_.lookup("k1"))),
00060 k2_(readScalar(coeffsDict_.lookup("k2"))),
00061 WeCrit_(readScalar(coeffsDict_.lookup("WeCrit"))),
00062 WeTransition_(readScalar(coeffsDict_.lookup("WeTransition"))),
00063 AWe_(0.0)
00064 {
00065 scalar k21 = k2_/k1_;
00066 AWe_ = (k21*sqrt(WeTransition_) - 1.0)/pow(WeTransition_, 4.0);
00067 }
00068
00069
00070
00071
00072 ETAB::~ETAB()
00073 {}
00074
00075
00076
00077
00078 void ETAB::breakupParcel
00079 (
00080 parcel& p,
00081 const scalar deltaT,
00082 const vector& Ug,
00083 const liquidMixture& fuels
00084 ) const
00085 {
00086
00087 scalar T = p.T();
00088 scalar pc = spray_.p()[p.cell()];
00089 scalar r = 0.5*p.d();
00090 scalar r2 = r*r;
00091 scalar r3 = r*r2;
00092
00093 scalar rho = fuels.rho(pc, T, p.X());
00094 scalar sigma = fuels.sigma(pc, T, p.X());
00095 scalar mu = fuels.mu(pc, T, p.X());
00096
00097
00098 scalar rtd = 0.5*Cmu_*mu/(rho*r2);
00099
00100
00101 scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
00102
00103 if (omega2 > 0)
00104 {
00105 scalar omega = sqrt(omega2);
00106 scalar romega = 1.0/omega;
00107
00108 scalar rhog = spray_.rho()[p.cell()];
00109 scalar We = p.We(Ug, rhog, sigma);
00110 scalar Wetmp = We/WeCrit_;
00111
00112 scalar y1 = p.dev() - Wetmp;
00113 scalar y2 = p.ddev()*romega;
00114
00115 scalar a = sqrt(y1*y1 + y2*y2);
00116
00117
00118 if (a+Wetmp > 1.0)
00119 {
00120 scalar phic = y1/a;
00121
00122
00123 phic = max(min(phic, 1), -1);
00124
00125 scalar phit = acos(phic);
00126 scalar phi = phit;
00127 scalar quad = -y2/a;
00128 if (quad < 0)
00129 {
00130 phi = 2*mathematicalConstant::pi - phit;
00131 }
00132
00133 scalar tb = 0;
00134
00135 if (mag(p.dev()) < 1.0)
00136 {
00137 scalar theta = acos((1.0 - Wetmp)/a);
00138
00139 if (theta < phi)
00140 {
00141 if (2*mathematicalConstant::pi-theta >= phi)
00142 {
00143 theta = -theta;
00144 }
00145 theta += 2*mathematicalConstant::pi;
00146 }
00147 tb = (theta-phi)*romega;
00148
00149
00150 if (deltaT > tb)
00151 {
00152 p.dev() = 1.0;
00153 p.ddev() = -a*omega*sin(omega*tb + phi);
00154 }
00155 }
00156
00157
00158 if (deltaT > tb)
00159 {
00160 scalar sqrtWe = AWe_*pow(We, 4.0) + 1.0;
00161 scalar Kbr = k1_*omega*sqrtWe;
00162
00163 if (We > WeTransition_)
00164 {
00165 sqrtWe = sqrt(We);
00166 Kbr =k2_*omega*sqrtWe;
00167 }
00168
00169 scalar rWetmp = 1.0/Wetmp;
00170 scalar cosdtbu = max(-1.0, min(1.0, 1.0-rWetmp));
00171 scalar dtbu = romega*acos(cosdtbu);
00172 scalar decay = exp(-Kbr*dtbu);
00173
00174 scalar rNew = decay*r;
00175 if (rNew < r)
00176 {
00177 p.d() = 2*rNew;
00178 p.dev() = 0;
00179 p.ddev() = 0;
00180 }
00181 }
00182 }
00183 }
00184 else
00185 {
00186
00187 p.dev() = 0;
00188 p.ddev() = 0;
00189 }
00190 }
00191
00192
00193
00194 }
00195
00196