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 particle`

tolerance = 0.00001; deltav = 1; % difference between calculated and known velocity

while deltav > tolerance deltavOLD = deltav; vNEW = ((4/3)*(d/CD)*g*((rhoP-rho)/rho))^0.5; % calculate new velocity using drag coefficient

deltavNEW = sqrt((vNEW - v)^2); % calculate difference between new and old velocities

deltav = deltavNEW; CD = CD+int; % increase drag coefficient by 0.00001 for next iteration

end

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?

## Best Answer