/* turbulence model: RSM with TCL pressure strain model by Craft (1997), (see Paper Craft T.J., Launder B.E., 2001)*/ /* UDF (User Defined Function) Routine for FLUENT */ /* see also the paper: */ /* Implementation and comparison of different turbulence models for three dimensional wall jets with Fluent */ /* FLUENT CFD Forum 2005 Bad Nauheim, Germany */ /* Heschl Ch*., Sanz W.**, Klanatsky P.* */ /* * Fachhochschulstudiengänge Burgenland GmbH, Steinamangerstraße 21, A-7423 Pinkafeld */ /* ** TU Graz, Institut of Thermal Turbomachinery and Machine Dynamics, Inffeldgasse 25, A-8010 Graz */ /* Contact: christian.heschl@fh-burgenland.at, wolfgang.sanz@tugraz.at, peter.klanatsky@fh-burgenland.at */ /* Die Autoren übernhemen keine Verantwortung oder Haftung für den Inhalt der unten angeführten UDF-Routinen */ #include "udf.h" #include "math.h" /* Turbulence model constants */ DEFINE_SOURCE(uu_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh11_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; //C_2s_new=MIN(C_2s_new,0.); //C_2s_new=MAX(C_2s_new,8.); /* C_1_new=3.1*sqrt(MAX(A*A2,0.)); C_2_new=0.55; C_1s_new=1.2; C_2s_new=0.6; */ PHIh11_1=-C_1_new*C_D(c,t)*(a11+C_1s_new*(a11*a11+a12*a21+a13*a31-1.0/3.0*A2))-C_D(c,t)*a11; PHIh11_2=-0.6*(P11-1.0/3.0*Pk_Pk)+0.3*a11*Pk_Pk; PHIh11_2=PHIh11_2 -0.2*(u11*u11/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u21/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u31/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u11/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u21/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u31/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u11/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u21/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u31/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u11*C_DUDX(c,t)+u11*C_DUDX(c,t))-u12/C_K(c,t)*(u12*C_DUDX(c,t)+u12*C_DUDX(c,t))-u13/C_K(c,t)*(u13*C_DUDX(c,t)+u13*C_DUDX(c,t))- u21/C_K(c,t)*(u11*C_DUDY(c,t)+u11*C_DUDY(c,t))-u22/C_K(c,t)*(u12*C_DUDY(c,t)+u12*C_DUDY(c,t))-u23/C_K(c,t)*(u13*C_DUDY(c,t)+u13*C_DUDY(c,t))- u31/C_K(c,t)*(u11*C_DUDZ(c,t)+u11*C_DUDZ(c,t))-u32/C_K(c,t)*(u12*C_DUDZ(c,t)+u12*C_DUDZ(c,t))-u33/C_K(c,t)*(u13*C_DUDZ(c,t)+u13*C_DUDZ(c,t))); PHIh11_2=PHIh11_2 -C_2_new*(A2*(P11-D11)+3.0*(a11*a11*(P11-D11)+a11*a21*(P12-D12)+a11*a31*(P13-D13)+ a21*a11*(P21-D21)+a21*a21*(P22-D22)+a21*a31*(P23-D23)+ a31*a11*(P31-D31)+a31*a21*(P32-D32)+a31*a31*(P33-D33))); PHIh11_3=(7.0/15.0-A2/4.0)*(P11-1.0/3.0*Pk_Pk); PHIh11_3 = PHIh11_3 + 0.1*(a11-0.5*(a11*a11+a12*a21+a13*a31-1./3.*A2))*Pk_Pk; PHIh11_3 = PHIh11_3 - 0.05*a11*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh11_3 = PHIh11_3 + 0.1*((u11/C_K(c,t)*P11+u11/C_K(c,t)*P11)+(u12/C_K(c,t)*P21+u12/C_K(c,t)*P21) +(u13/C_K(c,t)*P31+u13/C_K(c,t)*P31)-2./(3.*C_K(c,t)) *(u11*P11+u12*P21+u13*P31+u21*P12+u22*P22+u23*P32+u31*P13+u32*P23+u33*P33)); PHIh11_3 = PHIh11_3 + 0.2/SQR(C_K(c,t))*(u11*u11*(D11-P11)+u11*u21*(D12-P12)+u11*u31*(D13-P13)+ u21*u11*(D21-P21)+u21*u21*(D22-P22)+u21*u31*(D23-P23)+ u31*u11*(D31-P31)+u31*u21*(D32-P32)+u31*u31*(D33-P33)); PHIh11_3 = PHIh11_3 + 0.1*(6./SQR(C_K(c,t)))*(D11*u11*u11+D12*u11*u21+D13*u11*u31+D21*u21*u11 +D22*u21*u21+D23*u21*u31+D31*u31*u11+D32*u31*u21+D33*u31*u31); PHIh11_3 = PHIh11_3 - 0.2/SQR(C_K(c,t)) *(D11*u11*u11+D12*u11*u21+D13*u11*u31+D21*u21*u11+D22*u21*u21+D23*u21*u31+D31*u31*u11+D32*u31*u21+D33*u31*u31 +D11*u12*u12+D12*u12*u22+D13*u12*u32+D21*u22*u12+D22*u22*u22+D23*u22*u32+D31*u32*u12+D32*u32*u22+D33*u32*u32 +D11*u13*u13+D12*u13*u23+D13*u13*u33+D21*u23*u13+D22*u23*u23+D23*u23*u33+D31*u33*u13+D32*u33*u23+D33*u33*u33); PHIh11_3 = PHIh11_3 + 1.3/C_K(c,t)*(u11*u11*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u21*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u31*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u11*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u21*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u31*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u11*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u21*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u31*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh11_3 = PHIh11_3 - 1.3/(3.*C_K(c,t)) *(u11*u11*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u21*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u31*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u11*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u21*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u31*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u11*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u21*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u31*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u12*u12*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u22*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u32*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u12*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u22*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u32*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u12*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u22*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u32*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u13*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u33*(C_DWDZ(c,t)+C_DWDZ(c,t)) ); PHIh11_3=PHIh11_3*C_2s_new; PHI11=PHIh11_1+PHIh11_2+PHIh11_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a11/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI11; return source; } DEFINE_SOURCE(vv_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh22_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; //C_2s_new=MIN(C_2s_new,0.); //C_2s_new=MAX(C_2s_new,8.); /* C_1_new=3.1*sqrt(MAX(A*A2,0.)); C_2_new=0.55; C_1s_new=1.2; C_2s_new=0.6; */ PHIh22_1=-C_1_new*C_D(c,t)*(a22+C_1s_new*(a21*a12+a22*a22+a23*a32-1.0/3.0*A2))-C_D(c,t)*a22; PHIh22_2=-0.6*(P22-1.0/3.0*Pk_Pk)+0.3*a22*Pk_Pk; PHIh22_2=PHIh22_2 -0.2*(u12*u12/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u22/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u32/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u12/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u22/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u32/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u12/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u22/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u32/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u21*C_DVDX(c,t)+u21*C_DVDX(c,t))-u12/C_K(c,t)*(u22*C_DVDX(c,t)+u22*C_DVDX(c,t))-u13/C_K(c,t)*(u23*C_DVDX(c,t)+u23*C_DVDX(c,t))- u21/C_K(c,t)*(u21*C_DVDY(c,t)+u21*C_DVDY(c,t))-u22/C_K(c,t)*(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))-u23/C_K(c,t)*(u23*C_DVDY(c,t)+u23*C_DVDY(c,t))- u31/C_K(c,t)*(u21*C_DVDZ(c,t)+u21*C_DVDZ(c,t))-u32/C_K(c,t)*(u22*C_DVDZ(c,t)+u22*C_DVDZ(c,t))-u33/C_K(c,t)*(u23*C_DVDZ(c,t)+u23*C_DVDZ(c,t))); PHIh22_2=PHIh22_2 -C_2_new*(A2*(P22-D22)+3.0*(a12*a12*(P11-D11)+a12*a22*(P12-D12)+a12*a32*(P13-D13)+ a22*a12*(P21-D21)+a22*a22*(P22-D22)+a22*a32*(P23-D23)+ a32*a12*(P31-D31)+a32*a22*(P32-D32)+a32*a32*(P33-D33))); PHIh22_3=(7.0/15.0-A2/4.0)*(P22-1.0/3.0*Pk_Pk); PHIh22_3=PHIh22_3 + 0.1*(a22-0.5*(a21*a12+a22*a22+a23*a32-1./3.*A2))*Pk_Pk; PHIh22_3=PHIh22_3 - 0.05*a22*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh22_3=PHIh22_3 + 0.1*((u21/C_K(c,t)*P12+u21/C_K(c,t)*P12)+(u22/C_K(c,t)*P22+u22/C_K(c,t)*P22) +(u23/C_K(c,t)*P23+u23/C_K(c,t)*P23) -2./(3.*C_K(c,t))*(u11*P11+u12*P21+u13*P31+u21*P12+u22*P22+u23*P32+u31*P13+u32*P23+u33*P33)); PHIh22_3=PHIh22_3 + 0.2/SQR(C_K(c,t))*(u12*u12*(D11-P11)+u12*u22*(D12-P12)+u12*u32*(D13-P13) +u22*u12*(D21-P21)+u22*u22*(D22-P22)+u22*u32*(D23-P23) +u32*u12*(D31-P31)+u32*u22*(D32-P32)+u32*u32*(D33-P33)); PHIh22_3=PHIh22_3 + 0.1*(6./SQR(C_K(c,t)))*(D11*u12*u12+D12*u12*u22+D13*u12*u32+D21*u22*u12+D22*u22*u22+D23*u22*u32+D31*u32*u12+D32*u32*u22+D33*u32*u32); PHIh22_3=PHIh22_3 - 0.2/SQR(C_K(c,t))*(D11*u11*u11+D12*u11*u21+D13*u11*u31+D21*u21*u11+D22*u21*u21+D23*u21*u31+D31*u31*u11+D32*u31*u21+D33*u31*u31+ D11*u12*u12+D12*u12*u22+D13*u12*u32+D21*u22*u12+D22*u22*u22+D23*u22*u32+D31*u32*u12+D32*u32*u22+D33*u32*u32+ D11*u13*u13+D12*u13*u23+D13*u13*u33+D21*u23*u13+D22*u23*u23+D23*u23*u33+D31*u33*u13+D32*u33*u23+D33*u33*u33); PHIh22_3=PHIh22_3 + 1.3/C_K(c,t)*(u12*u12*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u22*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u32*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u12*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u22*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u32*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u12*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u22*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u32*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh22_3=PHIh22_3 - 1.3/(3.*C_K(c,t))*(u11*u11*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u21*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u31*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u11*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u21*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u31*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u11*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u21*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u31*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u12*u12*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u22*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u32*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u12*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u22*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u32*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u12*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u22*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u32*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u13*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u33*(C_DWDZ(c,t)+C_DWDZ(c,t)) ); PHIh22_3=C_2s_new*PHIh22_3; PHI22=PHIh22_1+PHIh22_2+PHIh22_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a22/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI22; return source; } DEFINE_SOURCE(ww_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh33_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; //C_2s_new=MIN(C_2s_new,0.); //C_2s_new=MAX(C_2s_new,8.); /* C_1_new=3.1*sqrt(MAX(A*A2,0.)); C_2_new=0.55; C_1s_new=1.2; C_2s_new=0.6; */ PHIh33_1=-C_1_new*C_D(c,t)*(a33+C_1s_new*(a31*a13+a32*a23+a33*a33-1.0/3.0*A2))-C_D(c,t)*a33; PHIh33_2=-0.6*(P33-1.0/3.0*Pk_Pk)+0.3*a33*Pk_Pk; PHIh33_2=PHIh33_2 -0.2*(u13*u13/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u23/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u33/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u13/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u23/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u33/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u13/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u23/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u33/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u31*C_DWDX(c,t)+u31*C_DWDX(c,t))-u32/C_K(c,t)*(u12*C_DWDX(c,t)+u32*C_DWDX(c,t))-u33/C_K(c,t)*(u13*C_DWDX(c,t)+u33*C_DWDX(c,t))- u21/C_K(c,t)*(u31*C_DWDY(c,t)+u31*C_DWDY(c,t))-u32/C_K(c,t)*(u12*C_DWDY(c,t)+u32*C_DWDY(c,t))-u33/C_K(c,t)*(u13*C_DWDY(c,t)+u33*C_DWDY(c,t))- u31/C_K(c,t)*(u31*C_DWDZ(c,t)+u31*C_DWDZ(c,t))-u32/C_K(c,t)*(u12*C_DWDZ(c,t)+u32*C_DWDZ(c,t))-u33/C_K(c,t)*(u13*C_DWDZ(c,t)+u33*C_DWDZ(c,t))); PHIh33_2=PHIh33_2 -C_2_new*(A2*(P33-D33)+3.0*(a13*a13*(P11-D11)+a13*a23*(P12-D12)+a13*a33*(P13-D13)+ a23*a13*(P21-D21)+a23*a23*(P22-D22)+a23*a33*(P23-D23)+ a33*a13*(P31-D31)+a33*a23*(P32-D32)+a33*a33*(P33-D33))); PHIh33_3=(7.0/15.0-A2/4.0)*(P33-1.0/3.0*Pk_Pk); PHIh33_3=PHIh33_3 + 0.1*(a33-0.5*(a31*a13+a32*a23+a33*a33-1./3.*A2))*Pk_Pk; PHIh33_3=PHIh33_3 - 0.05*a33*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh33_3=PHIh33_3 + 0.1*((u31/C_K(c,t)*P13+u31/C_K(c,t)*P13)+(u32/C_K(c,t)*P23+u32/C_K(c,t)*P23) +(u33/C_K(c,t)*P33+u33/C_K(c,t)*P33)-2./(3.*C_K(c,t)) *(u11*P11+u12*P21+u13*P31+u21*P12+u22*P22+u23*P32+u31*P13+u32*P23+u33*P33)); PHIh33_3=PHIh33_3 + 0.2/SQR(C_K(c,t))*(u13*u13*(D11-P11)+u13*u23*(D12-P12)+u13*u33*(D13-P13)+ u23*u13*(D21-P21)+u23*u23*(D22-P22)+u23*u33*(D23-P23)+ u33*u13*(D31-P31)+u33*u23*(D32-P32)+u33*u33*(D33-P33)); PHIh33_3=PHIh33_3 + 0.1*(6./SQR(C_K(c,t))) *(D11*u13*u13+D12*u13*u23+D13*u13*u33+D21*u23*u13+D22*u23*u23+D23*u23*u33+D31*u33*u13+D32*u33*u23+D33*u33*u33); PHIh33_3=PHIh33_3 - 0.2/SQR(C_K(c,t)) *(D11*u11*u11+D12*u11*u21+D13*u11*u31+D21*u21*u11+D22*u21*u21+D23*u21*u31+D31*u31*u11+D32*u31*u21+D33*u31*u31+ D11*u12*u12+D12*u12*u22+D13*u12*u32+D21*u22*u12+D22*u22*u22+D23*u22*u32+D31*u32*u12+D32*u32*u22+D33*u32*u32+ D11*u13*u13+D12*u13*u23+D13*u13*u33+D21*u23*u13+D22*u23*u23+D23*u23*u33+D31*u33*u13+D32*u33*u23+D33*u33*u33); PHIh33_3=PHIh33_3 + 1.3/C_K(c,t)*(u13*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u33*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh33_3=PHIh33_3 - 1.3/(3.*C_K(c,t)) *( u11*u11*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u21*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u31*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u11*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u21*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u31*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u11*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u21*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u31*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u12*u12*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u22*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u32*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u12*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u22*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u32*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u12*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u22*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u32*(C_DWDZ(c,t)+C_DWDZ(c,t))+ u13*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u33*(C_DWDZ(c,t)+C_DWDZ(c,t)) ); PHIh33_3=C_2s_new*PHIh33_3; PHI33=PHIh33_1+PHIh33_2+PHIh33_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a33/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI33; return source; } DEFINE_SOURCE(uv_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real PHI12, PHIh12_3, PHI13, PHIh13_3, PHI23, PHIh23_3; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh11_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; PHIh12_1=-C_1_new*C_D(c,t)*(a12+C_1s_new*(a11*a12+a12*a22+a13*a32))-C_D(c,t)*a12; PHIh12_2=-0.6*P12+0.3*a12*Pk_Pk; PHIh12_2=PHIh12_2 -0.2*(u12*u11/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u21/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u31/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u11/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u21/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u31/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u11/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u21/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u31/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u11*C_DVDX(c,t)+u21*C_DUDX(c,t))-u12/C_K(c,t)*(u12*C_DVDX(c,t)+u22*C_DUDX(c,t))-u13/C_K(c,t)*(u13*C_DVDX(c,t)+u23*C_DUDX(c,t))- u21/C_K(c,t)*(u11*C_DVDY(c,t)+u21*C_DUDY(c,t))-u22/C_K(c,t)*(u12*C_DVDY(c,t)+u22*C_DUDY(c,t))-u23/C_K(c,t)*(u13*C_DVDY(c,t)+u23*C_DUDY(c,t))- u31/C_K(c,t)*(u11*C_DVDZ(c,t)+u21*C_DUDZ(c,t))-u32/C_K(c,t)*(u12*C_DVDZ(c,t)+u22*C_DUDZ(c,t))-u33/C_K(c,t)*(u13*C_DVDZ(c,t)+u23*C_DUDZ(c,t))); PHIh12_2=PHIh12_2 -C_2_new*(A2*(P12-D12)+3.0*(a11*a12*(P11-D11)+a11*a22*(P12-D12)+a11*a32*(P13-D13)+ a21*a12*(P21-D21)+a21*a22*(P22-D22)+a21*a32*(P23-D23)+ a31*a12*(P31-D31)+a31*a22*(P32-D32)+a31*a32*(P33-D33))); PHIh12_3=(7.0/15.0-A2/4.0)*P12; PHIh12_3 = PHIh12_3 + 0.1*(a12-0.5*(a11*a12+a12*a22+a13*a32))*Pk_Pk; PHIh12_3 = PHIh12_3 - 0.05*a12*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh12_3 = PHIh12_3 + 0.1*((u11/C_K(c,t)*P12+u21/C_K(c,t)*P11)+(u12/C_K(c,t)*P22+u22/C_K(c,t)*P21)+(u13/C_K(c,t)*P32+u23/C_K(c,t)*P31)); PHIh12_3 = PHIh12_3 + 0.2/SQR(C_K(c,t))*(u11*u12*(D11-P11)+u11*u22*(D12-P12)+u11*u32*(D13-P13)+ u21*u12*(D21-P21)+u21*u22*(D22-P22)+u21*u32*(D23-P23)+ u31*u12*(D31-P31)+u31*u22*(D32-P32)+u31*u32*(D33-P33)); PHIh12_3 = PHIh12_3 + 0.1*(6./SQR(C_K(c,t)))*(D11*u11*u12+D12*u11*u22+D13*u11*u32+D21*u21*u12+D22*u21*u22+D23*u21*u32+D31*u31*u12+D32*u31*u22+D33*u31*u32); PHIh12_3 = PHIh12_3 + 1.3/C_K(c,t)*(u11*u12*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u22*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u32*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u12*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u22*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u32*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u12*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u22*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u32*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh12_3=PHIh12_3*C_2s_new; PHI12=PHIh12_1+PHIh12_2+PHIh12_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a12/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI12; return source; } DEFINE_SOURCE(uw_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real PHI12, PHIh12_3, PHI13, PHIh13_3, PHI23, PHIh23_3; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh11_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; PHIh13_1=-C_1_new*C_D(c,t)*(a13+C_1s_new*(a11*a13+a12*a23+a13*a33))-C_D(c,t)*a13; PHIh13_2=-0.6*P13+0.3*a13*Pk_Pk; PHIh13_2=PHIh13_2 -0.2*(u13*u11/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u21/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u31/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u11/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u21/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u31/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u11/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u21/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u31/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u11*C_DWDX(c,t)+u31*C_DUDX(c,t))-u12/C_K(c,t)*(u12*C_DWDX(c,t)+u32*C_DUDX(c,t))-u13/C_K(c,t)*(u13*C_DWDX(c,t)+u33*C_DUDX(c,t))- u21/C_K(c,t)*(u11*C_DWDY(c,t)+u31*C_DUDY(c,t))-u22/C_K(c,t)*(u12*C_DWDY(c,t)+u32*C_DUDY(c,t))-u23/C_K(c,t)*(u13*C_DWDY(c,t)+u33*C_DUDY(c,t))- u31/C_K(c,t)*(u11*C_DWDZ(c,t)+u31*C_DUDZ(c,t))-u32/C_K(c,t)*(u12*C_DWDZ(c,t)+u32*C_DUDZ(c,t))-u33/C_K(c,t)*(u13*C_DWDZ(c,t)+u33*C_DUDZ(c,t))); PHIh13_2=PHIh13_2 -C_2_new*(A2*(P13-D13)+3.0*(a11*a13*(P11-D11)+a11*a23*(P12-D12)+a11*a33*(P13-D13)+ a21*a13*(P21-D21)+a21*a23*(P22-D22)+a21*a33*(P23-D23)+ a31*a13*(P31-D31)+a31*a23*(P32-D32)+a31*a33*(P33-D33))); PHIh13_3=(7.0/15.0-A2/4.0)*P13; PHIh13_3 = PHIh13_3 + 0.1*(a13-0.5*(a11*a13+a12*a23+a13*a33))*Pk_Pk; PHIh13_3 = PHIh13_3 - 0.05*a13*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh13_3 = PHIh13_3 + 0.1*((u11/C_K(c,t)*P13+u31/C_K(c,t)*P11)+(u12/C_K(c,t)*P23+u32/C_K(c,t)*P21)+(u13/C_K(c,t)*P33+u33/C_K(c,t)*P31)); PHIh13_3 = PHIh13_3 + 0.2/SQR(C_K(c,t))*(u11*u13*(D11-P11)+u11*u23*(D12-P12)+u11*u33*(D13-P13)+ u21*u13*(D21-P21)+u21*u23*(D22-P22)+u21*u33*(D23-P23)+ u31*u13*(D31-P31)+u31*u23*(D32-P32)+u31*u33*(D33-P33)); PHIh13_3 = PHIh13_3 + 0.1*(6./SQR(C_K(c,t)))*(D11*u11*u13+D12*u11*u23+D13*u11*u33+D21*u21*u13+D22*u21*u23+D23*u21*u33+D31*u31*u13+D32*u31*u23+D33*u31*u33); PHIh13_3 = PHIh13_3 + 1.3/C_K(c,t)*(u11*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u11*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u11*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u21*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u21*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u21*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u31*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u31*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u31*u33*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh13_3=PHIh13_3*C_2s_new; PHI13=PHIh13_1+PHIh13_2+PHIh13_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a13/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI13; return source; } DEFINE_SOURCE(vw_source, c, t, dS, eqn) { real source; real A, A2, A3; real P11, P12, P13, P21, P22, P23, P31, P32, P33; real u11, u12, u13, u21, u22, u23, u31, u32, u33; real S11, S12, S13, S21, S22, S23, S31, S32, S33; real W11, W12, W13, W21, W22, W23, W31, W32, W33; real uk_uk, Pk_Pk; real a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, b21, b22, b23, b31, b32, b33; real Re_t; real C_1_new, C_2_new, C_1s_new, C_2s_new, PHIw11, PHIw22, PHIw33, PHI11, PHI22, PHI33; real PHIh11_1, PHIh12_1,PHIh13_1,PHIh21_1,PHIh22_1,PHIh23_1,PHIh31_1,PHIh32_1,PHIh33_1; real PHIh11_2, PHIh12_2,PHIh13_2,PHIh21_2,PHIh22_2,PHIh23_2,PHIh31_2,PHIh32_2,PHIh33_2; real PHI12, PHIh12_3, PHI13, PHIh13_3, PHI23, PHIh23_3; real D11, D12, D13, D21, D22, D23, D31, D32, D33, S, W, S_I, PHIh22_3; u11=C_RUU(c,t); u12=C_RUV(c,t); u13=C_RUW(c,t); u21=C_RUV(c,t); u22=C_RVV(c,t); u23=C_RVW(c,t); u31=C_RUW(c,t); u32=C_RVW(c,t); u33=C_RWW(c,t); S11 = ( C_DUDX(c,t)+C_DUDX(c,t) ); S12 = ( C_DUDY(c,t)+C_DVDX(c,t) ); S13 = ( C_DUDZ(c,t)+C_DWDX(c,t) ); S21 = ( C_DVDX(c,t)+C_DUDY(c,t) ); S22 = ( C_DVDY(c,t)+C_DVDY(c,t) ); S23 = ( C_DVDZ(c,t)+C_DWDY(c,t) ); S31 = ( C_DWDX(c,t)+C_DUDZ(c,t) ); S32 = ( C_DWDY(c,t)+C_DVDZ(c,t) ); S33 = ( C_DWDZ(c,t)+C_DWDZ(c,t) ); W11 = ( C_DUDX(c,t)-C_DUDX(c,t) ); W12 = ( C_DUDY(c,t)-C_DVDX(c,t) ); W13 = ( C_DUDZ(c,t)-C_DWDX(c,t) ); W21 = ( C_DVDX(c,t)-C_DUDY(c,t) ); W22 = ( C_DVDY(c,t)-C_DVDY(c,t) ); W23 = ( C_DVDZ(c,t)-C_DWDY(c,t) ); W31 = ( C_DWDX(c,t)-C_DUDZ(c,t) ); W32 = ( C_DWDY(c,t)-C_DVDZ(c,t) ); W33 = ( C_DWDZ(c,t)-C_DWDZ(c,t) ); P11=C_R(c,t)*(-u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t) -u11*C_DUDX(c,t)-u12*C_DUDY(c,t)-u13*C_DUDZ(c,t)); P12=C_R(c,t)*(-u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t) -u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t)); P13=C_R(c,t)*(-u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t) -u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t)); P21=C_R(c,t)*(-u21*C_DUDX(c,t)-u22*C_DUDY(c,t)-u23*C_DUDZ(c,t) -u11*C_DVDX(c,t)-u12*C_DVDY(c,t)-u13*C_DVDZ(c,t)); P22=C_R(c,t)*(-u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t) -u21*C_DVDX(c,t)-u22*C_DVDY(c,t)-u23*C_DVDZ(c,t)); P23=C_R(c,t)*(-u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t) -u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t)); P31=C_R(c,t)*(-u31*C_DUDX(c,t)-u32*C_DUDY(c,t)-u33*C_DUDZ(c,t) -u11*C_DWDX(c,t)-u12*C_DWDY(c,t)-u13*C_DWDZ(c,t)); P32=C_R(c,t)*(-u31*C_DVDX(c,t)-u32*C_DVDY(c,t)-u33*C_DVDZ(c,t) -u21*C_DWDX(c,t)-u22*C_DWDY(c,t)-u23*C_DWDZ(c,t)); P33=C_R(c,t)*(-u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t) -u31*C_DWDX(c,t)-u32*C_DWDY(c,t)-u33*C_DWDZ(c,t)); uk_uk=u11+u22+u33; Pk_Pk=P11+P22+P33; a11=-((-u11+2./3.*C_K(c,t))/C_K(c,t)); a12=-(-u12/C_K(c,t)); a13=-(-u13/C_K(c,t)); a21=-(-u21/C_K(c,t)); a22=-((-u22+2./3.*C_K(c,t))/C_K(c,t)); a23=-(-u23/C_K(c,t)); a31=-(-u31/C_K(c,t)); a32=-(-u32/C_K(c,t)); a33=-((-u33+2./3.*C_K(c,t))/C_K(c,t)); A2=a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+ a31*a13+a32*a23+a33*a33; A3=a11*a11*a11+a12*a21*a11+a13*a31*a11+ a21*a11*a12+a22*a21*a12+a23*a31*a12+ a31*a11*a13+a32*a21*a13+a33*a31*a13+ a11*a12*a21+a12*a22*a21+a13*a32*a21+ a21*a12*a22+a22*a22*a22+a23*a32*a22+ a31*a12*a23+a32*a22*a23+a33*a32*a23+ a11*a13*a31+a12*a23*a31+a13*a33*a31+ a21*a13*a32+a22*a23*a32+a23*a33*a32+ a31*a13*a33+a32*a23*a33+a33*a33*a33; A=(1.0-9./8.*(A2-A3)); A=MAX(A,1.e-12); b11=-((-u11+2./3.*C_K(c,t))/(2.*C_K(c,t))); b12=-(-u12/(2.*C_K(c,t))); b13=-(-u13/(2.*C_K(c,t))); b21=-(-u21/(2.*C_K(c,t))); b22=-((-u22+2./3.*C_K(c,t))/(2.*C_K(c,t))); b23=-(-u23/(2.*C_K(c,t))); b31=-(-u31/(2.*C_K(c,t))); b32=-(-u32/(2.*C_K(c,t))); b33=-((-u33+2./3.*C_K(c,t))/(2.*C_K(c,t))); /*S=S11*S11+S12*S21+S13*S31+ S21*S12+S22*S22+S23*S32+ S31*S13+S32*S23+S33*S33;*/ /*W=W11*W11+W12*W21+W13*W31+ W21*W12+W22*W22+W23*W32+ W31*W13+W32*W23+W33*W33;*/ S=S11*S11+S12*S12+S13*S13+ S21*S21+S22*S22+S23*S23+ S31*S31+S32*S32+S33*S33; W=W11*W11+W12*W12+W13*W13+ W21*W21+W22*W22+W23*W23+ W31*W31+W32*W32+W33*W33; //S=MAX(S,0.0); //W=MAX(W,0.0); //S=S11*S11+S22*S22+S33*S33; //W=W11*W11+W22*W22+W33*W33; S=C_K(c,t)/C_D(c,t)*sqrt(S/2.0); W=C_K(c,t)/C_D(c,t)*sqrt(W/2.0); S_I=(S11*S11*S11+S11*S12*S21+S11*S13*S31+ S12*S21*S11+S12*S22*S21+S12*S23*S31+ S13*S31*S11+S13*S32*S21+S13*S33*S31+ S21*S11*S12+S21*S12*S22+S21*S13*S32+ S22*S21*S12+S22*S22*S22+S22*S23*S32+ S23*S31*S12+S23*S32*S22+S23*S33*S32+ S31*S11*S13+S31*S12*S23+S31*S13*S33+ S32*S21*S13+S32*S22*S23+S32*S23*S33+ S33*S31*S13+S33*S32*S23+S33*S33*S33)/ pow( (S11*S11+S12*S12+S13*S13+S21*S21+S22*S22+S23*S23+S31*S31+S32*S32+S33*S33)/2. , 1.5 ); D11=-( (u11*C_DUDX(c,t)+u11*C_DUDX(c,t))+(u12*C_DVDX(c,t)+u12*C_DVDX(c,t))+(u13*C_DWDX(c,t)+u13*C_DWDX(c,t)) ); D12=-( (u11*C_DUDY(c,t)+u21*C_DUDX(c,t))+(u12*C_DVDY(c,t)+u22*C_DVDX(c,t))+(u13*C_DWDY(c,t)+u23*C_DWDX(c,t)) ); D13=-( (u11*C_DUDZ(c,t)+u31*C_DUDX(c,t))+(u12*C_DVDZ(c,t)+u32*C_DVDX(c,t))+(u13*C_DWDZ(c,t)+u33*C_DWDX(c,t)) ); D21=-( (u21*C_DUDX(c,t)+u11*C_DUDY(c,t))+(u22*C_DVDX(c,t)+u12*C_DVDY(c,t))+(u23*C_DWDX(c,t)+u13*C_DWDY(c,t)) ); D22=-( (u21*C_DUDY(c,t)+u21*C_DUDY(c,t))+(u22*C_DVDY(c,t)+u22*C_DVDY(c,t))+(u23*C_DWDY(c,t)+u23*C_DWDY(c,t)) ); D23=-( (u21*C_DUDZ(c,t)+u31*C_DUDY(c,t))+(u22*C_DVDZ(c,t)+u32*C_DVDY(c,t))+(u23*C_DWDZ(c,t)+u33*C_DWDY(c,t)) ); D31=-( (u31*C_DUDX(c,t)+u11*C_DUDZ(c,t))+(u32*C_DVDX(c,t)+u12*C_DVDZ(c,t))+(u33*C_DWDX(c,t)+u13*C_DWDZ(c,t)) ); D32=-( (u31*C_DUDY(c,t)+u21*C_DUDZ(c,t))+(u32*C_DVDY(c,t)+u22*C_DVDZ(c,t))+(u33*C_DWDY(c,t)+u23*C_DWDZ(c,t)) ); D33=-( (u31*C_DUDZ(c,t)+u31*C_DUDZ(c,t))+(u32*C_DVDZ(c,t)+u32*C_DVDZ(c,t))+(u33*C_DWDZ(c,t)+u33*C_DWDZ(c,t)) ); Re_t=C_R(c,t)*SQR(C_K(c,t))/(C_MU_L(c,t)*C_D(c,t)); C_1_new=3.1*pow(MAX(A*A2,0.),0.5); C_2_new=MIN(0.55*(1.0-exp(-pow(MAX(A,0.0),1.5)*Re_t/100.0)),3.2*MAX(A,0.)/(1.+S)); C_1s_new=1.1; C_2s_new=MIN(0.6,MAX(A,0.))+3.5*(S-W)/(3.0+S+W)-2.0*S_I; PHIh23_1=-C_1_new*C_D(c,t)*(a23+C_1s_new*(a21*a13+a22*a23+a23*a33))-C_D(c,t)*a23; PHIh23_2=-0.6*P23+0.3*a23*Pk_Pk; PHIh23_2=PHIh23_2 -0.2*(u13*u12/C_K(c,t)*(C_DUDX(c,t)+C_DUDX(c,t))+u13*u22/C_K(c,t)*(C_DUDY(c,t)+C_DVDX(c,t))+u13*u32/C_K(c,t)*(C_DUDZ(c,t)+C_DWDX(c,t))+ u23*u12/C_K(c,t)*(C_DVDX(c,t)+C_DUDY(c,t))+u23*u22/C_K(c,t)*(C_DVDY(c,t)+C_DVDY(c,t))+u23*u32/C_K(c,t)*(C_DVDZ(c,t)+C_DWDY(c,t))+ u33*u12/C_K(c,t)*(C_DWDX(c,t)+C_DUDZ(c,t))+u33*u22/C_K(c,t)*(C_DWDY(c,t)+C_DVDZ(c,t))+u33*u32/C_K(c,t)*(C_DWDZ(c,t)+C_DWDZ(c,t))- u11/C_K(c,t)*(u21*C_DWDX(c,t)+u31*C_DVDX(c,t))-u12/C_K(c,t)*(u22*C_DWDX(c,t)+u32*C_DVDX(c,t))-u13/C_K(c,t)*(u23*C_DWDX(c,t)+u33*C_DVDX(c,t))- u21/C_K(c,t)*(u21*C_DWDY(c,t)+u31*C_DVDY(c,t))-u22/C_K(c,t)*(u22*C_DWDY(c,t)+u32*C_DVDY(c,t))-u23/C_K(c,t)*(u23*C_DWDY(c,t)+u33*C_DVDY(c,t))- u31/C_K(c,t)*(u21*C_DWDZ(c,t)+u31*C_DVDZ(c,t))-u32/C_K(c,t)*(u22*C_DWDZ(c,t)+u32*C_DVDZ(c,t))-u33/C_K(c,t)*(u23*C_DWDZ(c,t)+u33*C_DVDZ(c,t))); PHIh23_2=PHIh23_2 -C_2_new*(A2*(P23-D23)+3.0*(a12*a13*(P11-D11)+a12*a23*(P12-D12)+a12*a33*(P13-D13)+ a22*a13*(P21-D21)+a22*a23*(P22-D22)+a22*a33*(P23-D23)+ a32*a13*(P31-D31)+a32*a23*(P32-D32)+a32*a33*(P33-D33))); PHIh23_3=(7.0/15.0-A2/4.0)*P23; PHIh23_3=PHIh23_3 + 0.1*(a23-0.5*(a21*a13+a22*a23+a23*a33))*Pk_Pk; PHIh23_3=PHIh23_3 - 0.05*a23*(a11*P11+a12*P21+a13*P31+a21*P12+a22*P22+a23*P32+a31*P13+a32*P23+a33*P33); PHIh23_3=PHIh23_3 + 0.1*((u21/C_K(c,t)*P13+u31/C_K(c,t)*P12)+(u22/C_K(c,t)*P23+u32/C_K(c,t)*P22)+(u23/C_K(c,t)*P33+u33/C_K(c,t)*P32)); PHIh23_3=PHIh23_3 + 0.2/SQR(C_K(c,t))*(u12*u13*(D11-P11)+u12*u23*(D12-P12)+u12*u33*(D13-P13) +u22*u13*(D21-P21)+u22*u23*(D22-P22)+u22*u33*(D23-P23) +u32*u13*(D31-P31)+u32*u23*(D32-P32)+u32*u33*(D33-P33)); PHIh23_3=PHIh23_3 + 0.1*(6./SQR(C_K(c,t)))*(D11*u12*u13+D12*u12*u23+D13*u12*u33+D21*u22*u13+D22*u22*u23+D23*u22*u33+D31*u32*u13+D32*u32*u23+D33*u32*u33); PHIh23_3=PHIh23_3 + 1.3/C_K(c,t)*(u12*u13*(C_DUDX(c,t)+C_DUDX(c,t))+u12*u23*(C_DUDY(c,t)+C_DVDX(c,t))+u12*u33*(C_DUDZ(c,t)+C_DWDX(c,t))+ u22*u13*(C_DVDX(c,t)+C_DUDY(c,t))+u22*u23*(C_DVDY(c,t)+C_DVDY(c,t))+u22*u33*(C_DVDZ(c,t)+C_DWDY(c,t))+ u32*u13*(C_DWDX(c,t)+C_DUDZ(c,t))+u32*u23*(C_DWDY(c,t)+C_DVDZ(c,t))+u32*u33*(C_DWDZ(c,t)+C_DWDZ(c,t))); PHIh23_3=C_2s_new*PHIh23_3; PHI23=PHIh23_1+PHIh23_2+PHIh23_3; dS[eqn]=-C_1_new*C_D(c,t)*(1./C_K(c,t)+C_1s_new*(2.*a23/C_K(c,t)))-C_D(c,t)/C_K(c,t); source=PHI23; return source; }