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 "TAB.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <OpenFOAM/mathematicalConstants.H>
00029
00030
00031
00032 namespace Foam
00033 {
00034
00035
00036
00037 defineTypeNameAndDebug(TAB, 0);
00038
00039 addToRunTimeSelectionTable
00040 (
00041 breakupModel,
00042 TAB,
00043 dictionary
00044 );
00045
00046
00047
00048
00049 TAB::TAB
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 WeCrit_(readScalar(coeffsDict_.lookup("WeCrit")))
00060 {
00061
00062
00063 const scalar xx0 = 12.0;
00064 const scalar rrd100 = 1.0/(1.0-exp(-xx0)*(1+xx0+pow(xx0, 2)/2+pow(xx0, 3)/6));
00065
00066 for(label n=0; n<100; n++)
00067 {
00068 scalar xx = 0.12*(n+1);
00069 rrd_[n] = (1-exp(-xx)*(1 + xx + pow(xx, 2)/2 + pow(xx, 3)/6))*rrd100;
00070 }
00071 }
00072
00073
00074
00075
00076 TAB::~TAB()
00077 {}
00078
00079
00080
00081
00082 void TAB::breakupParcel
00083 (
00084 parcel& p,
00085 const scalar deltaT,
00086 const vector& Ug,
00087 const liquidMixture& fuels
00088 ) const
00089 {
00090
00091 scalar T = p.T();
00092 scalar pc = spray_.p()[p.cell()];
00093 scalar r = 0.5*p.d();
00094 scalar r2 = r*r;
00095 scalar r3 = r*r2;
00096
00097 scalar rho = fuels.rho(pc, T, p.X());
00098 scalar sigma = fuels.sigma(pc, T, p.X());
00099 scalar mu = fuels.mu(pc, T, p.X());
00100
00101
00102 scalar rtd = 0.5*Cmu_*mu/(rho*r2);
00103
00104
00105 scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;
00106
00107 if (omega2 > 0)
00108 {
00109 scalar omega = sqrt(omega2);
00110 scalar rhog = spray_.rho()[p.cell()];
00111 scalar We = p.We(Ug, rhog, sigma);
00112 scalar Wetmp = We/WeCrit_;
00113
00114 scalar y1 = p.dev() - Wetmp;
00115 scalar y2 = p.ddev()/omega;
00116
00117 scalar a = sqrt(y1*y1 + y2*y2);
00118
00119
00120 if (a+Wetmp > 1.0)
00121 {
00122 scalar phic = y1/a;
00123
00124
00125 phic = max(min(phic, 1), -1);
00126
00127 scalar phit = acos(phic);
00128 scalar phi = phit;
00129 scalar quad = -y2/a;
00130 if (quad < 0)
00131 {
00132 phi = 2*mathematicalConstant::pi - phit;
00133 }
00134
00135 scalar tb = 0;
00136
00137 if (mag(p.dev()) < 1.0)
00138 {
00139 scalar coste = 1.0;
00140 if
00141 (
00142 (Wetmp - a < -1)
00143 && (p.ddev() < 0)
00144 )
00145 {
00146 coste = -1.0;
00147 }
00148
00149 scalar theta = acos((coste-Wetmp)/a);
00150
00151 if (theta < phi)
00152 {
00153 if (2*mathematicalConstant::pi-theta >= phi)
00154 {
00155 theta = -theta;
00156 }
00157 theta += 2*mathematicalConstant::pi;
00158 }
00159 tb = (theta-phi)/omega;
00160
00161
00162 if (deltaT > tb)
00163 {
00164 p.dev() = 1.0;
00165 p.ddev() = -a*omega*sin(omega*tb + phi);
00166 }
00167
00168 }
00169
00170
00171 if (deltaT > tb)
00172 {
00173 scalar rs = r/
00174 (
00175 1
00176 + (4.0/3.0)*pow(p.dev(), 2)
00177 + rho*r3/(8*sigma)*pow(p.ddev(), 2)
00178 );
00179
00180 label n = 0;
00181 bool found = false;
00182 scalar random = rndGen_.scalar01();
00183 while (!found && (n<99))
00184 {
00185 if (rrd_[n]>random)
00186 {
00187 found = true;
00188 }
00189 n++;
00190
00191 }
00192 scalar rNew = 0.04*n*rs;
00193 if (rNew < r)
00194 {
00195 p.d() = 2*rNew;
00196 p.dev() = 0;
00197 p.ddev() = 0;
00198 }
00199
00200 }
00201
00202 }
00203
00204 }
00205 else
00206 {
00207
00208 p.dev() = 0;
00209 p.ddev() = 0;
00210 }
00211
00212 }
00213
00214
00215
00216 }
00217
00218