|
HG-MD
1
|
|
HG-MD
1
|
00001 #include "Chute.h" 00002 #include "Statistics.h" 00003 00010 class ChuteWithHopper : public Chute { 00011 protected: 00013 Mdouble HopperLength; 00015 Mdouble HopperHeight; 00017 Mdouble HopperAngle; 00019 Mdouble HopperExitLength; 00021 Mdouble HopperExitHeight; 00023 Mdouble shift; 00025 Mdouble lowerFillHeight; 00027 bool centerHopper; 00028 private: 00030 Mdouble lift; 00031 00033 unsigned int hopper_dim; 00034 00036 bool align_base; 00037 00039 Mdouble fill_percent; 00040 00041 00042 public: 00043 00045 ChuteWithHopper(Chute& other) : MD(other), Chute(other) {constructor();} 00046 ChuteWithHopper(HGRID_3D& other) : MD(other), Chute(other) {constructor();} 00047 ChuteWithHopper(HGRID_base& other) : MD(other), Chute(other) {constructor();} 00048 ChuteWithHopper(MD& other) : MD(other), Chute(other) {constructor();} 00049 00051 ChuteWithHopper() {constructor();} 00052 00054 void constructor() 00055 { 00056 lowerFillHeight=0.5; 00057 lift=0.0; 00058 set_Hopper(0.01, 0.01, 60.0, 0.08); 00059 shift = 0.0; 00060 hopper_dim=1; 00061 align_base=true; 00062 00063 fill_percent=50.0; 00064 centerHopper=false; 00065 00066 } 00067 00068 void set_HopperFillPercentage(Mdouble new_fill){fill_percent=new_fill;} 00069 00070 virtual void setup_particles_initial_conditions() 00071 { 00072 Chute::setup_particles_initial_conditions(); 00073 //cout << shift << " " << get_xmax() << " " << get_ParticleHandler().get_NumberOfParticles() << endl; 00074 add_hopper(); 00075 } 00076 00078 void add_hopper() { 00079 //hopper walls 00080 int n = get_NWall(); 00081 set_NWall(n+2); 00082 00083 //to create the finite hopper walls, we take vector between two wall points in xz-plane, then rotate clockwise and make unit length 00084 // A\ /A 00085 // \ / A,B,C denote three points on the left and right hopper walls which are used to construct the hopper 00086 // \ / shift denotes the space by which the chute has to be shifted to the right such that the hopper is in the domain 00087 // B| |B 00088 // | | 00089 // | | 00090 // C| |C 00091 00092 Vec3D A, B, C, temp, normal; 00093 Mdouble s = sin(get_ChuteAngle()); 00094 Mdouble c = cos(get_ChuteAngle()); 00095 Mdouble HopperLowestPoint = HopperExitHeight - HopperExitLength * tan(ChuteAngle); 00096 HopperHeight = HopperLowestPoint + 1.1 * 0.5*(HopperLength+HopperExitLength) / tan(HopperAngle); 00097 Mdouble HopperCornerHeight = HopperHeight - 0.5*(HopperLength-HopperExitLength) / tan(HopperAngle); 00098 if (HopperCornerHeight<=0.0) { HopperHeight += -HopperCornerHeight + P0.get_Radius(); HopperCornerHeight = P0.get_Radius(); } 00099 00100 //first we create the left hopper wall 00101 00102 //coordinates of A,B,C in (vertical parallel to flow,vertical normal to flow, horizontal) direction 00103 A = Vec3D(-0.5*(HopperLength-HopperExitLength), 0.0, HopperHeight); 00104 B = Vec3D(0.0, 0.0, HopperCornerHeight); 00105 C = Vec3D(0.0, 0.0, 0.0); 00106 00107 00108 00109 //now rotate the coordinates of A,B,C to be in (x,y,z) direction 00110 A = Vec3D(c*A.X-s*A.Z, 0.0, s*A.X+c*A.Z); 00111 B = Vec3D(c*B.X-s*B.Z, 0.0, s*B.X+c*B.Z); 00112 C = Vec3D(c*C.X-s*C.Z, 0.0, s*C.X+c*C.Z); 00113 // the position of A determines shift and zmax 00114 if (centerHopper) set_shift(-A.X+40); 00115 else set_shift(-A.X); 00116 set_zmax(A.Z); 00117 A.X +=shift; 00118 B.X +=shift; 00119 C.X +=shift; 00120 00121 //This lifts the hopper a distance above the chute 00122 A.Z+=lift; 00123 B.Z+=lift; 00124 C.Z+=lift; 00125 00126 //create a finite wall from B to A and from C to B 00127 temp = B-A; 00128 normal = Vec3D(temp.Z,0.0,-temp.X) / sqrt(temp.GetLength2()); 00129 Walls[n].add_finite_wall(normal, Dot(normal,A)); 00130 temp = C-B; 00131 normal = Vec3D(temp.Z,0.0,-temp.X) / sqrt(temp.GetLength2()); 00132 Walls[n].add_finite_wall(normal, Dot(normal,B)); 00133 temp = A-C; 00134 normal = Vec3D(temp.Z,0.0,-temp.X)/sqrt(temp.GetLength2()); 00135 Walls[n].add_finite_wall(normal,Dot(normal,C)); 00136 00137 00138 00139 //next, do the same for the right wall 00140 A = Vec3D(HopperLength-0.5*(HopperLength-HopperExitLength), 0.0, HopperHeight); 00141 B = Vec3D(0.5*(HopperLength+HopperExitLength)-0.5*(HopperLength-HopperExitLength), 0.0, HopperCornerHeight); 00142 C = Vec3D(0.5*(HopperLength+HopperExitLength)-0.5*(HopperLength-HopperExitLength), 0.0, HopperLowestPoint); 00143 00144 00145 00146 00147 //This rotates the right points 00148 A = Vec3D(c*A.X-s*A.Z+shift, 0.0, s*A.X+c*A.Z); 00149 B = Vec3D(c*B.X-s*B.Z+shift, 0.0, s*B.X+c*B.Z); 00150 C = Vec3D(c*C.X-s*C.Z+shift, 0.0, s*C.X+c*C.Z); 00151 00152 //This lifts the hopper a distance above the chute 00153 A.Z+=lift; 00154 B.Z+=lift; 00155 C.Z+=lift; 00156 00157 temp = A-B; 00158 normal = Vec3D(temp.Z,0.0,-temp.X) / sqrt(temp.GetLength2()); 00159 Walls[n+1].add_finite_wall(normal, Dot(normal,A)); 00160 temp = B-C; 00161 normal = Vec3D(temp.Z,0.0,-temp.X) / sqrt(temp.GetLength2()); 00162 Walls[n+1].add_finite_wall(normal, Dot(normal,B)); 00163 temp = C-A; 00164 normal = Vec3D(temp.Z,0.0,-temp.X) / sqrt(temp.GetLength2()); 00165 Walls[n+1].add_finite_wall(normal, Dot(normal,C)); 00166 00167 set_zmax(A.Z); 00168 00169 if (hopper_dim == 2) 00170 { 00171 00172 set_NWall(n+4); 00173 00174 00175 //coordinates of A,B,C in (vertical parallel to flow,vertical normal to flow, horizontal) direction 00176 A = Vec3D(0.0, (get_ymax()-get_ymin()-HopperLength)/2.0, HopperHeight); 00177 B = Vec3D(0.0, (get_ymax()-get_ymin()-HopperExitLength)/2.0, HopperCornerHeight); 00178 C = Vec3D(0.0, (get_ymax()-get_ymin()-HopperExitLength)/2.0, 0.0); 00179 00180 00181 00182 //now rotate the coordinates of A,B,C to be in (x,y,z) direction 00183 A = Vec3D(c*A.X-s*A.Z, A.Y, s*A.X+c*A.Z); 00184 B = Vec3D(c*B.X-s*B.Z, B.Y, s*B.X+c*B.Z); 00185 C = Vec3D(c*C.X-s*C.Z, C.Y, s*C.X+c*C.Z); 00186 // the position of A determines shift and zmax 00187 A.X +=shift; 00188 B.X +=shift; 00189 C.X +=shift; 00190 00191 //This lifts the hopper a distance above the chute 00192 A.Z+=lift; 00193 B.Z+=lift; 00194 C.Z+=lift; 00195 00196 00197 00198 //create a finite wall from B to A and from C to B 00199 temp = B-A; 00200 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00201 //normal = Vec3D(0.0,temp.Z,-temp.Y) / sqrt(temp.GetLength2()); 00202 Walls[n+2].add_finite_wall(normal, Dot(normal,A)); 00203 temp = C-B; 00204 //normal = Vec3D(0.0,temp.Z,-temp.Y) / sqrt(temp.GetLength2()); 00205 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00206 Walls[n+2].add_finite_wall(normal, Dot(normal,B)); 00207 temp = A-C; 00208 //normal = Vec3D(0.0,temp.Z,-temp.Y)/sqrt(temp.GetLength2()); 00209 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00210 Walls[n+2].add_finite_wall(normal,Dot(normal,C)); 00211 00212 00213 //Now for the right y-wall 00214 A = Vec3D(0.0, (get_ymax()-get_ymin()+HopperLength)/2.0,HopperHeight); 00215 B = Vec3D(0.0, (get_ymax()-get_ymin()+HopperExitLength)/2.0,HopperCornerHeight); 00216 C = Vec3D(0.0, (get_ymax()-get_ymin()+HopperExitLength)/2.0,0.0); 00217 00218 //now rotate the coordinates of A,B,C to be in (x,y,z) direction 00219 A = Vec3D(c*A.X-s*A.Z, A.Y, s*A.X+c*A.Z); 00220 B = Vec3D(c*B.X-s*B.Z, B.Y, s*B.X+c*B.Z); 00221 C = Vec3D(c*C.X-s*C.Z, C.Y, s*C.X+c*C.Z); 00222 // the position of A determines shift and zmax 00223 A.X +=shift; 00224 B.X +=shift; 00225 C.X +=shift; 00226 00227 //This lifts the hopper a distance above the chute 00228 A.Z+=lift; 00229 B.Z+=lift; 00230 C.Z+=lift; 00231 00232 //create a finite wall from B to A and from C to B 00233 temp = A-B; 00234 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00235 //normal = Vec3D(0.0,-temp.Z,temp.Y) / sqrt(temp.GetLength2()); 00236 Walls[n+3].add_finite_wall(normal, Dot(normal,A)); 00237 temp = B-C; 00238 //normal = Vec3D(0.0,-temp.Z,temp.Y) / sqrt(temp.GetLength2()); 00239 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00240 Walls[n+3].add_finite_wall(normal, Dot(normal,B)); 00241 temp = C-A; 00242 //normal = Vec3D(0.0,-temp.Z,temp.Y)/sqrt(temp.GetLength2()); 00243 normal=Cross(Vec3D(-c,0,-s),temp)/sqrt(temp.GetLength2()); 00244 Walls[n+3].add_finite_wall(normal,Dot(normal,C)); 00245 00246 00247 00248 00249 00250 00251 00252 00253 } 00254 00255 00256 00257 00258 00259 00260 00261 //now shift the fixed particles at the bottom so that they begin where the chute begins 00262 for (vector<Particle*>::iterator it= get_ParticleHandler().begin(); it!=get_ParticleHandler().end(); ++it) { 00263 (*it)->get_Position().X+=shift; 00264 #ifdef USE_SIMPLE_VERLET_INTEGRATION 00265 (*it)->PrevPosition.X +=shift; 00266 #endif 00267 } 00268 } 00269 00270 00276 virtual void create_inflow_particle() 00277 { 00278 00279 00280 //use this formula to obtain bidispersed particles 00281 //P0.Radius = random(0.0,1.0)<0.1?MinInflowParticleRadius:MaxInflowParticleRadius; 00282 00283 //the following formula yields polydispersed particle radii: 00284 00285 P0.set_Radius(random.get_RN(MinInflowParticleRadius,MaxInflowParticleRadius)); 00286 P0.compute_particle_mass(Species); 00287 00288 //Define a orthogonal coordinate system this is usful in the hopper, see diagram in html documentation for details. 00289 static Mdouble s = sin(get_ChuteAngle()); 00290 static Mdouble c = cos(get_ChuteAngle()); 00291 static Mdouble Ht = tan(HopperAngle); 00292 static Mdouble Hc = cos(HopperAngle); 00293 static Vec3D AB = Vec3D(c,0.0,s); 00294 static Vec3D AC = Vec3D(-s,0.0,c); 00295 static Vec3D AD = Vec3D(0.0,1.0,0.0); 00296 00297 //Point A is located in the centre of the hopper. 00298 static Vec3D A = Vec3D 00299 ( 00300 centerHopper?40:0.0, 00301 (get_ymax()-get_ymin())/2.0, 00302 s*(-0.5*(HopperLength-HopperExitLength)) + c*HopperHeight 00303 ) 00304 + AB*0.5*HopperLength 00305 + AC*(-0.5*HopperLength/Ht); 00306 00307 Mdouble gamma = random.get_RN((100.0-fill_percent)/100.0,1.0); 00308 00309 00310 // Mdouble gamma = random(lowerFillHeight,1.0); 00311 00312 Mdouble delta; 00313 00314 if (hopper_dim==1) 00315 { 00316 00317 00318 //For the one dimensional delta is a random distance between the two walls the -minus 2 particle radii is to stop 00322 //delta = random(ymin+P0.Radius,ymax-P0.Radius); 00323 00324 delta = random.get_RN(-0.5,0.5)*(get_ymax()-get_ymin()-2.0*P0.get_Radius()); 00325 00326 00327 } 00328 else 00329 { 00330 00331 delta= (random.get_RN(-1.0,1.0)*(0.5*gamma*HopperLength -P0.get_Radius()/Hc)); 00332 00333 } 00334 P0.set_Position( A 00335 + AC * (gamma*0.5*HopperLength/Ht) 00336 00337 + AB * (random.get_RN(-1.0,1.0)*(0.5*gamma*HopperLength - P0.get_Radius()/Hc)) 00338 + AD*delta); 00339 00340 00341 P0.get_Position().Z+=lift; 00342 00343 //P0.Position.Y = random(ymin+P0.Radius, ymax-P0.Radius); 00344 //P0.Position.Y=delta; 00345 00346 00347 00348 } 00349 00350 void set_Hopper(Mdouble ExitLength, Mdouble ExitHeight, Mdouble Angle, Mdouble Length){ 00351 if (ExitLength>=0.0) {HopperExitLength = ExitLength;} else cerr << "WARNING : Hopper exit length must be greater than or equal to zero" << endl; 00352 if (ExitHeight>=0.0) {HopperExitHeight = ExitHeight;} else cerr << "WARNING : Hopper exit height must be greater than or equal to zero" << endl; 00353 if (Angle>0.0&&Angle<90.0) {HopperAngle = Angle*constants::pi/180.0;} else cerr << "WARNING : Hopper angle must in (0,90)" << endl; 00354 if (Length>ExitLength) {HopperLength = Length;} else cerr << "WARNING : Hopper length must be greater than exit length" << endl; 00355 } 00356 00358 00359 Mdouble get_MaximumVelocityInducedByGravity(){ 00360 Mdouble height = HopperHeight+(get_xmax()-shift)*sin(ChuteAngle); 00361 00362 return sqrt(2.0*9.8*height); 00363 } 00364 00366 00367 Mdouble get_ChuteLength(){return get_xmax()-shift;} 00368 00370 void set_ChuteLength(Mdouble new_){if (new_>=0.0) {set_xmax(new_+shift); set_xmin(0.0);} else cerr << "WARNING : Chute length unchanged, value must be greater than or equal to zero" << endl;} 00371 00372 void set_centerHopper(bool new_){centerHopper=new_; } 00373 00374 void set_lowerFillHeight(Mdouble new_){lowerFillHeight=new_; } 00375 00376 void set_shift(Mdouble new_){if (new_>=0.0) {set_xmax(get_xmax()+new_-shift); shift = new_;} else cerr << "WARNING : Shift length unchanged, value must be greater than or equal to zero" << endl;} 00377 00379 virtual void read(std::istream& is) { 00380 Chute::read(is); 00381 is >> HopperExitLength >> HopperExitHeight >> HopperLength 00382 >> HopperAngle >> HopperHeight >> shift; 00383 } 00384 00386 virtual void write(std::ostream& os) { 00387 Chute::write(os); 00388 os << HopperExitLength << " " << HopperExitHeight << " " << HopperLength 00389 << " " << HopperAngle << " " << HopperHeight << " " << shift << " " << endl; 00390 } 00391 00392 virtual void print(std::ostream& os) { 00393 Chute::print(os); 00394 os 00395 << ", HopperExitLength:" << HopperExitLength 00396 << ", HopperExitHeight:" << HopperExitHeight 00397 << ", HopperLength:" << HopperLength 00398 << ", HopperAngle:" << HopperAngle 00399 << ", HopperHeight:" << HopperHeight 00400 << endl; 00401 } 00402 00404 void lift_hopper(Mdouble distance){lift=distance;} 00405 Mdouble get_lift_hopper(){return lift;} 00406 00407 void set_hopper_dim(Mdouble new_hopper_dim){hopper_dim=new_hopper_dim;} 00408 00409 void set_align_base(bool new_align){align_base=new_align;} 00410 00411 int readNextArgument(unsigned int& i, unsigned int argc, char *argv[]) { 00412 if (!strcmp(argv[i],"-hopperLength")) { 00413 HopperLength=(atof(argv[i+1])); 00414 } else if (!strcmp(argv[i],"-hopperHeight")) { 00415 HopperHeight=(atof(argv[i+1])); 00416 } else if (!strcmp(argv[i],"-hopperAngle")) { 00417 HopperAngle=(atof(argv[i+1])); 00418 } else if (!strcmp(argv[i],"-hopperExitLength")) { 00419 HopperExitLength=(atof(argv[i+1])); 00420 } else if (!strcmp(argv[i],"-hopperExitHeight")) { 00421 HopperExitHeight=(atof(argv[i+1])); 00422 } else if (!strcmp(argv[i],"-lowerFillHeight")) { 00423 lowerFillHeight=(atof(argv[i+1])); 00424 } else if (!strcmp(argv[i],"-centerHopper")) { 00425 centerHopper=(atoi(argv[i+1])); 00426 } else if (!strcmp(argv[i],"-alignBase")) { 00427 align_base=(atoi(argv[i+1])); 00428 } else if (!strcmp(argv[i],"-shift")) { 00429 shift=(atof(argv[i+1])); 00430 } else if (!strcmp(argv[i],"-lift")) { 00431 lift=(atof(argv[i+1])); 00432 } else return Chute::readNextArgument(i, argc, argv); //if argv[i] is not found, check the commands in Chute 00433 return true; //returns true if argv[i] is found 00434 } 00435 00436 };
1.7.6.1
