MATLAB: Iteration to convergence for drag coefficient – is there a better method

convergencedragfluid dynamicsiterationparticle settlingterminal velocity

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

Why don't you simply solve your equation for CD ?
CD = 4/3*d/v^2*g*(rhoP/rho-1)
Best wishes
Torsten.
Related Question