00001
00002
00003 #ifndef LBM_SEPPLEMENTALBOUNDARY_H
00004 #define LBM_SEPPLEMENTALBOUNDARY_H
00005
00019
00020
00021
00022
00023
00024 class BoundaryConditionsSupplemental : public LBMBoundaryConditions<LBMType,DIM> {
00025 typedef LBMBoundaryConditions<LBMType,DIM> base;
00026 public:
00027 typedef base::MicroType MicroType;
00028 typedef base::MacroType MacroType;
00029 BoundaryConditionsSupplemental(LBMType &lbm) : base(lbm) {}
00030 enum BCSupp { EqInletRamp=int(LBMType::NoSlipWallEntranceExit+1), oneSeventhInletRamp, parabolicInletRamp, slidingWallRamp };
00031
00032
00033 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00034 base::SetupData(gh,ghosts);
00035 const double* bndry = base::GH().wholebndry();
00036 const int nbndry = base::GH().nbndry();
00037 DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00038 if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00039 DataType omega = LBM().Omega(LBM().TimeScale());
00040 DataType u_p_avg=Side(0).aux[0];
00041 double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00042 double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00043 ( comm_service::log() << "l0_p=" << l0_p << " u_p_avg=" << u_p_avg
00044 << " Re_p=" << Re_p << " Re_lbm=" << Re_lbm
00045 << "\n omega=" << omega
00046 << " SpeedUp=" << LBM().SpeedUp()
00047 << " u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00048 << " SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00049 << " LBM().TO()=" << LBM().T0()
00050 << "\n Re_p-Re_lbm=" << Re_p-Re_lbm
00051 << std::endl ).flush();
00052 }
00053
00054 virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00055 const int &dir, const double& time) {
00056
00057 base::SetBndry(fvec,level,bb,dir,time);
00058 BBox wholebbox = base::GH().wholebbox();
00059 BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00060 if (!inside.empty())
00061 LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00062
00063 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00064 int b = fvec.bottom();
00065 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00066 && fvec.stepsize(2)==bb.stepsize(2));
00067 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00068 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00069 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00070 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00071 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00072 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00073 MicroType *f = (MicroType *)fvec.databuffer(), ftmp;
00074 MacroType q;
00075
00077 if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00078 Side(0).Type==EqInletRamp) {
00079
00080 DataType rho=Side(0).aux[0];
00081 DataType a0=Side(0).aux[1];
00082 DataType a1=Side(0).aux[2];
00083 DataType a2=Side(0).aux[3];
00084 DataType T_ramp=Side(0).aux[4];
00085 DataType uinit=Side(0).aux[5];
00086 MacroType q, q_;
00087 DataType ramp;
00088
00089 if (Side(0).Scaling==LBMType::Physical) {
00090 rho /= LBM().DensityScale();
00091 a0 /= LBM().VelocityScale();
00092 a1 /= LBM().VelocityScale();
00093 a2 /= LBM().VelocityScale();
00094 uinit /= LBM().VelocityScale();
00095 }
00096
00097 q_(0) = rho;
00098 q_(1) = a0;
00099 q_(2) = a1;
00100 q_(3) = a2;
00101
00102 if ( time < T_ramp ) {
00103 ramp = time/T_ramp;
00104 q(0) = q_(0);
00105 q(1) = uinit+(q_(1)-uinit)*ramp;
00106 q(2) = q_(2)*ramp;
00107 q(3) = q_(3)*ramp;
00108 }
00109 else
00110 q = q_;
00111 for (register int k=mzs; k<=mze; k++) {
00112 for (register int j=mys; j<=mye; j++) {
00113 q(0) = LBM().MacroVariables(f[b+LBM().idx(mxe+1,j,k,Nx,Ny)])(0);
00114 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00115 }
00116 }
00117 }
00118
00120 if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00121 Side(0).Type==oneSeventhInletRamp) {
00122
00123 DataType rho=Side(0).aux[0];
00124 DataType a0=Side(0).aux[1];
00125 DataType a1=Side(0).aux[2];
00126 DataType a2=Side(0).aux[3];
00127 DataType elevation=Side(0).aux[4];
00128
00129 MacroType q, q_;
00130 DataType profile, height;
00131 DCoords dx = base::GH().worldStep(fvec.stepsize());
00132 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), 0 };
00133
00134 if (Side(0).Scaling==LBMType::Physical) {
00135 rho /= LBM().DensityScale();
00136 a0 /= LBM().VelocityScale();
00137 a1 /= LBM().VelocityScale();
00138 a2 /= LBM().VelocityScale();
00139 }
00140
00141 q_(0) = rho;
00142 q_(1) = a0;
00143 q_(2) = a1;
00144 q_(3) = a2;
00145
00146 for (register int k=mzs; k<=mze; k++) {
00147 lc_indx[2] = k*fvec.stepsize(2);
00148 DCoords wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00149 height = wc(2)+dx(2)/2.;
00150
00151 if ( height < elevation ) {
00152 profile = std::pow((height/elevation),(1./7.));
00153 q(0) = q_(0);
00154 q(1) = q_(1)*profile;
00155 q(2) = q_(2)*profile;
00156 q(3) = q_(3)*profile;
00157 }
00158 else
00159 q = q_;
00160
00161 for (register int j=mys; j<=mye; j++) {
00162 ftmp = LBM().Equilibrium(q);
00163 for (register int qi=0; qi<19; qi++)
00164 f[b+LBM().idx(mxe,j,k,Nx,Ny)](qi) = ftmp(qi);
00165 }
00166
00167 }
00168 }
00169
00171 if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00172 Side(0).Type==parabolicInletRamp) {
00173
00174 DataType rho=Side(0).aux[0];
00175 DataType a0=Side(0).aux[1];
00176 DataType a1=Side(0).aux[2];
00177 DataType a2=Side(0).aux[3];
00178 DataType T_ramp=Side(0).aux[4];
00179 DataType H=Side(0).aux[5];
00180
00181 MacroType q, q_;
00182 DataType profile, height;
00183 DataType ramp, pi=std::asin(1.);
00184 DCoords dx = base::GH().worldStep(fvec.stepsize());
00185 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), 0 };
00186
00187 if (Side(0).Scaling==LBMType::Physical) {
00188 rho /= LBM().DensityScale();
00189 a0 /= LBM().VelocityScale();
00190 a1 /= LBM().VelocityScale();
00191 a2 /= LBM().VelocityScale();
00192 }
00193
00194 q_(0) = rho;
00195 q_(1) = a0;
00196 q_(2) = a1;
00197 q_(3) = a2;
00198
00199 if ( time < T_ramp ) {
00200 ramp = std::tanh(time/T_ramp*4.*pi);
00201 q_(0) = q_(0);
00202 q_(1) *= ramp;
00203 q_(2) *= ramp;
00204 q_(3) *= ramp;
00205 }
00206 else
00207 q = q_;
00208
00209 for (register int k=mzs; k<=mze; k++) {
00210 lc_indx[2] = k*fvec.stepsize(2);
00211 DCoords wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00212 height = wc(2)+dx(2)/2.+H/2.;
00213 if ( height < 0.0) height = 0.0;
00214 profile = 4.*height*(H-height)/(H*H);
00215 q(0) = q_(0);
00216 q(1) = q_(1)*profile;
00217 q(2) = q_(2)*profile;
00218 q(3) = q_(3)*profile;
00219 for (register int j=mys; j<=mye; j++) {
00220 f[b+LBM().idx(mxe,j,k,Nx,Ny)] = LBM().Equilibrium(q);
00221 }
00222 }
00223 }
00224
00226 if (dir==LBMType::Top && bb.lower(1)-wholebbox.stepsize(1)==wholebbox.upper(1) &&
00227 Side(3).Type==slidingWallRamp) {
00228
00229 DataType rho=Side(3).aux[0];
00230 DataType a0=Side(3).aux[1];
00231 DataType a1=Side(3).aux[2];
00232 DataType a2=Side(3).aux[3];
00233 DataType T_ramp=Side(3).aux[4];
00234 MacroType q, q_;
00235 DataType ramp, pi=std::asin(1.);
00236
00237 if (Side(3).Scaling==LBMType::Physical) {
00238 rho /= LBM().DensityScale();
00239 a0 /= LBM().VelocityScale();
00240 a1 /= LBM().VelocityScale();
00241 a2 /= LBM().VelocityScale();
00242 }
00243
00244 q_(0) = rho;
00245 q_(1) = a0;
00246 q_(2) = a1;
00247 q_(3) = a2;
00248
00249 if ( time < T_ramp ) {
00250 ramp = std::tanh(time/T_ramp*4.*pi);
00251 q(0) = q_(0);
00252 q(1) = q_(1)*ramp;
00253 q(2) = q_(2)*ramp;
00254 q(3) = q_(3)*ramp;
00255 }
00256 else
00257 q = q_;
00258
00259 for (register int k=mzs; k<=mze; k++)
00260 for (register int i=mxs; i<=mxe; i++) {
00261 MacroType ql = LBM().MacroVariables(f[b+LBM().idx(i,mys-1,k,Nx,Ny)]);
00262 DataType rd1 = f[b+LBM().idx(i,mys-1,k,Nx,Ny)](7);
00263 DataType rd2 = f[b+LBM().idx(i,mys-1,k,Nx,Ny)](10);
00264 DataType qw = (ql(0)/DataType(6.)*q(1))/(rd1-rd2);
00265 DataType pw = DataType(1.)-qw;
00266 if (i<mxe) f[b+LBM().idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
00267 f[b+LBM().idx(i,mys,k,Nx,Ny)](4) = f[b+LBM().idx(i,mys-1,k,Nx,Ny)](3);
00268 if (i>mxs) f[b+LBM().idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
00269 if (k>mzs) f[b+LBM().idx(i,mys,k,Nx,Ny)](12) = f[b+LBM().idx(i,mys-1,k-1,Nx,Ny)](11);
00270 if (k<mze) f[b+LBM().idx(i,mys,k,Nx,Ny)](14) = f[b+LBM().idx(i,mys-1,k+1,Nx,Ny)](13);
00271 }
00272 }
00273 }
00274
00275 };
00276
00277 #endif