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 "MixedDiffuseSpecular.H"
00027
00028
00029
00030 template <class CloudType>
00031 Foam::MixedDiffuseSpecular<CloudType>::MixedDiffuseSpecular
00032 (
00033 const dictionary& dict,
00034 CloudType& cloud
00035 )
00036 :
00037 WallInteractionModel<CloudType>(dict, cloud, typeName),
00038 diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction")))
00039 {}
00040
00041
00042
00043
00044 template <class CloudType>
00045 Foam::MixedDiffuseSpecular<CloudType>::~MixedDiffuseSpecular()
00046 {}
00047
00048
00049
00050
00051 template <class CloudType>
00052 void Foam::MixedDiffuseSpecular<CloudType>::correct
00053 (
00054 const wallPolyPatch& wpp,
00055 const label faceId,
00056 vector& U,
00057 scalar& Ei,
00058 label typeId
00059 )
00060 {
00061 label wppIndex = wpp.index();
00062
00063 label wppLocalFace = wpp.whichFace(faceId);
00064
00065 vector nw = wpp.faceAreas()[wppLocalFace];
00066
00067
00068 nw /= mag(nw);
00069
00070
00071 scalar U_dot_nw = U & nw;
00072
00073 CloudType& cloud(this->owner());
00074
00075 Random& rndGen(cloud.rndGen());
00076
00077 if (diffuseFraction_ > rndGen.scalar01())
00078 {
00079
00080
00081
00082 vector Ut = U - U_dot_nw*nw;
00083
00084 while (mag(Ut) < SMALL)
00085 {
00086
00087
00088
00089
00090 U = vector
00091 (
00092 U.x()*(0.8 + 0.2*rndGen.scalar01()),
00093 U.y()*(0.8 + 0.2*rndGen.scalar01()),
00094 U.z()*(0.8 + 0.2*rndGen.scalar01())
00095 );
00096
00097 U_dot_nw = U & nw;
00098
00099 Ut = U - U_dot_nw*nw;
00100 }
00101
00102
00103 vector tw1 = Ut/mag(Ut);
00104
00105
00106 vector tw2 = nw^tw1;
00107
00108 scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
00109
00110 scalar mass = cloud.constProps(typeId).mass();
00111
00112 scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
00113
00114 U =
00115 sqrt(CloudType::kb*T/mass)
00116 *(
00117 rndGen.GaussNormal()*tw1
00118 + rndGen.GaussNormal()*tw2
00119 - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
00120 );
00121
00122 U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
00123
00124 Ei = cloud.equipartitionInternalEnergy(T, iDof);
00125 }
00126 else
00127 {
00128
00129
00130 if (U_dot_nw > 0.0)
00131 {
00132 U -= 2.0*U_dot_nw*nw;
00133 }
00134 }
00135
00136 }
00137
00138
00139