Hi, I am having problems iterating to find drag coefficient of falling particles. I know the terminal velocity (v) and can therefore calculate the Reynolds number (Re). Terminal velocity is a function of drag coefficient CD:
v = ((4/3)*(d/CD)*g*((rhoP-rho)/rho))^0.5
where rhoP is particle density, rho is fluid density, g is the gravitational constant and d is particle diameter.
There is an analytical solution for CD for falling spheres. Since I know that CD for my particles is higher than that of spheres, I am using this equation as an initial guess and then iterating, increasing CD by a set increment every time, until the calculated terminal velocity is equal to known terminal velocity (within a tolerance).
My code looks like this:
CD = C1 + (24/Re) + (C2/(1+(sqrt(Re)))); % analytical solution based on spherical particletolerance = 0.00001;deltav = 1; % difference between calculated and known velocitywhile deltav > tolerancedeltavOLD = deltav;vNEW = ((4/3)*(d/CD)*g*((rhoP-rho)/rho))^0.5; % calculate new velocity using drag coefficientdeltavNEW = sqrt((vNEW - v)^2); % calculate difference between new and old velocitiesdeltav = deltavNEW;CD = CD+int; % increase drag coefficient by 0.00001 for next iterationend
The problem is, this is very computationally expensive when the drag coefficient differs by many orders of magnitude from the initial estimate. Also, there are cases where increasing CD by the chosen interval does not result in a velocity solution to within tolerance.
Does anyone know of an iteration method which would be better suited to solving this problem?