46 typename std::deque<PARTICLETYPE<T> >::iterator p,
int pInt,
51 std::cout << _converter.getConversionFactorForce() << std::endl;
52 T fluidVel[3] = {0., 0., 0.};
53 std::cout << _physDensity <<
"phys density" <<std::endl;
54 std::cout << _dynVisc <<
"phys density" <<std::endl;
55 _getVel(fluidVel, &p->getPos()[0], p->getCuboid());
56 std::cout <<
"Fluidvel: " << fluidVel[0] <<
" " << fluidVel[1] <<
" " << fluidVel[2] << std::endl;
58 std::cout <<
"rel_velo " << rel_velo << std::endl;
61 T particleRe = 2. * p->getRad() * rel_velo * _physDensity / _dynVisc;
62std::cout <<
"radius " << p->getRad() << std::endl;
63 std::cout <<
"particleRe" << particleRe << std::endl;
65 if(particleRe < 0.000001){
68 coeffCD = 24./particleRe*(1. + p->getA()*
util::pow(particleRe, p->getB()))+ p->getC()/(1.+p->getD()/particleRe);
70 std::cout <<
"coeffCD " << coeffCD << std::endl;
73 for (
int i = 0; i < 3; i++) {
75 force[i] = -1. * 0.5* coeffCD *
M_PI * p->getRad()*p->getRad()*_physDensity*rel_velo*(p->getVel()[i]-fluidVel[i]);
76 std::cout <<
"force " << i <<
" " << force[i] << std::endl;
77 p->getForce()[i] += force[i] ;
80std::cout <<
"apply force finished " << std::endl;