An error occurred while loading the file. Please try again.
-
Dorchies David authored8eb14f35
function [res]=find_Q_nat(cote_bas,ks,D,h,Cd0,S,L,pf,C,Q,sigma)
fprintf('*************************************\n')
fprintf('find_Q_nat Q=%f\n',Q)
maxfun=5000;
maxiter=5000;
tolfun=1e-16;
tolx=1e-16;
g=9.81;
kappa=0.41;
U0=Q./L./pf;
Fr=U0./(9.81*pf).^0.5;
Frg=Fr/(1-C^0.5);
if ks==0;ks=1e-34;end
coeff_contraction=0.4*Cd0+0.7;
%if Cd0==2
% fFr=(min(2.5,Frg^(-4/3)));%
%else
fFr=(min(coeff_contraction./(1-(Frg.^2)/4),Frg.^(-2/3))).^2;
%end
if Frg>1.5
fFr=(Frg.^(-2/3)).^2;
end
alpha=1-(1.*C).^0.5-1/2*sigma.*C;
Cd=Cd0.*(1+1./(pf./D).^2).*fFr;
%Cd=Cd0.*(0.8-2*C).*(1+0.4./(pf./D).^2).*fFr;
R=(1-sigma*C);%.*(1-C.^0.5).^2;
if pf/h>1.1; %fFr=1;
choixturb=1;
htilde=h./D;
hstar=pf./h;
Rh=h.*(hstar-1);
ustar=(g.*S.*Rh).^0.5;
fprintf('ustar=%f\n',ustar)
Cd1=Cd0.*(1+0.4./(pf./D).^2).*fFr;
CdCh=Cd1.*C.*htilde;
fprintf('CdCh=%f\n',CdCh)
Cf=2./(5.1.*log10(h./ks)+6).^2;
U0b=(2*g.*S.*R./(Cd1.*C.*h./D+Cf.*R).*h).^0.5;
fprintf('U0b=%f\n',U0b)
% U0=(2*g.*S.*D.*R./(Cd1.*C)).^0.5;
[P]=fminbnd(@(alphai) resolve_alpha(alphai,CdCh,R,U0b,hstar,h,C,D,Cd1,ustar,choixturb),1e-5*h,h,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx));
alpha=P(1);
fprintf('alpha=%f\n',alpha)
beta2=h.*CdCh./alpha./R;
beta=(beta2).^0.5;
a1=beta*(hstar-1)/(cosh(beta));
c=1;
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116
UhU0=(a1*sinh(beta)+c)^0.5;
Uh=UhU0*U0b;
fprintf('Uh=%f\n',Uh)
dhp=1-1/kappa*alpha./h.*Uh./ustar;
z0hp=(1-dhp).*exp(-1*(kappa*Uh./ustar));
qsup=ustar./kappa.*h.*((hstar-dhp).*(log((hstar-dhp)./z0hp) - 1)-((1-dhp).*(log((1-dhp)./z0hp) - 1)));
fprintf('qsup=%f\n',qsup)
%calcul intgrale dans la canope----
dzinf=0.01;
Zinf=(0:dzinf:1);
Uinf=U0b .*(beta.*Rh./h.*sinh(beta*Zinf)./cosh(beta)+1).^0.5;
Ub=zeros(size(Uinf));
Ub(1:end-1)=Uinf(2:end);
qinf=sum((Uinf(1:end-1)+Ub(1:end-1))/2*dzinf.*h);
fprintf('qinf=%f\n',qinf)
qtot=qinf+qsup;
PI=0;
delta=1;
Umoy=qtot./pf;
res=abs(U0-Umoy);
else
hstar=pf/D;
Re=U0.*pf/1e-6;
if ks==0
Cf=0.3164/4.*Re.^(-0.25);
else
Cf=2/(5.1*log10(pf/ks-1)+6)^2;
end
N= (alpha.*Cf)./(pf./D.*Cd.*C);
res=abs(U0-(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5);
end
fprintf('res=%f\n',res)