画海表粗糙度的迭代程序

liu541226

U0=zeros(81,1);
U0=0:80;
max=100;
for m=1:80
    k=1;
    Z0(k)=1.59*10^(-5);
    piancha=abs(Z0(k)-0.1);
    while(k<=max)
        ux=0.4*U0(m)/log(10/Z0(k));           %Z0=10/exp(0.4*U0/ux);
        Z0(k+1)=fun(ux);
        piancha=abs(Z0(k+1)-Z0(k));
        if(piancha<10^(-5))
            disp('successfully'),m;
            Z0_final(m)=Z0(k+1);
            break
        end
        k=k+1;
    end
    if(piancha>=10^(-5))
        disp('Failing');
        break
    end
end
plot(Z0_final);
hold on;
plot(Z0_final/2);
hold on;
plot(1.59E10^(-5));
hold on;
%-----------------------------

function Z0=fun(ux)
Zw=min(1,(ux/1.06)^0.3);
Z1=0.011*(ux^2/9.81)+1.59*10^-5;
Z2=10/exp(9.5*ux^(-1/3))+1.65*10^(-6)/max(ux,0.01);
Z0=max(1.27*10^(-7),min(Zw*Z2+(1-Zw)*Z1,2.85*10^-3));
end