Commit 8eb14f35 by Dorchies David

### #22 Validation des tests scilab / code matlab

parent c6f3597f
 function [res]=find_Q_nat2(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg) //fonctin pour ax=ay 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; ... ... @@ -10,79 +9,87 @@ tolx=1e-16; g=9.81; kappa=0.41; U0=Q./L./pf; if(bDbg) then print_r("U0"); end; 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 %if Cd0==2 % fFr=(min(2.5,Frg^(-4/3)));% %else fFr=(min(coeff_contraction./(1-(Frg.^2)/4),Frg.^(-2/3))).^2; //end %end if Frg>1.5 fFr=(Frg.^(-2/3)).^2; end if(bDbg) then print_r("fFr"); end; alpha=1-(1.*C).^0.5-1/2*sigma.*C; //Cd=Cd0.*(1+1./(pf./D).^2).*fFr; 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; %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; 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; CdCh=Cd0.*(1+0.4./(pf./D).^2).*C.*htilde; u0=(2*g.*S.*D.*R./(Cd0.*C)).^0.5; // [P]=fminbnd(@(alphai) resolve_alpha(alphai,CdCh,R,u0,hstar,h,C,D,Cd0,ustar,choixturb),1e-5*h,h,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); fprintf('ustar=%f\n',ustar) Cd1=Cd0.*(1+0.4./(pf./D).^2).*fFr; CdCh=Cd1.*C.*htilde; fprintf('CdCh=%f\n',CdCh) //[P]=fminbnd(@(alphai) resolve_alpha(alphai,CdCh,R,u0,hstar,h,C,D,Cd0,ustar,choixturb),1e-5*h,h); [alpha fval] = fminsearch(list(resolve_alpha, CdCh,R,u0,hstar,h,C,D,Cd0,ustar), 1e-5*h) 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; beta2=(beta2).^0.5; a1=beta2*(hstar-1)/(cosh(beta2)); beta=(beta2).^0.5; a1=beta*(hstar-1)/(cosh(beta)); c=1; UhU0=(a1*sinh(beta2)+c)^0.5; Uh=UhU0*u0; 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))); //calcul intégrale dans la canopée---- U(1)=u0; U(2)=u0.*(beta2.*Rh./h.*sinh(beta2*0.1)./cosh(beta2)+c).^0.5; U(3)=u0.*(beta2.*Rh./h.*sinh(beta2*0.2)./cosh(beta2)+c).^0.5; U(4)=u0.*(beta2.*Rh./h.*sinh(beta2*0.3)./cosh(beta2)+c).^0.5; U(5)=u0.*(beta2.*Rh./h.*sinh(beta2*0.4)./cosh(beta2)+c).^0.5; U(6)=u0.*(beta2.*Rh./h.*sinh(beta2*0.5)./cosh(beta2)+c).^0.5; U(7)=u0.*(beta2.*Rh./h.*sinh(beta2*0.6)./cosh(beta2)+c).^0.5; U(8)=u0.*(beta2.*Rh./h.*sinh(beta2*0.7)./cosh(beta2)+c).^0.5; U(9)=u0.*(beta2.*Rh./h.*sinh(beta2*0.8)./cosh(beta2)+c).^0.5; U(10)=u0.*(beta2.*Rh./h.*sinh(beta2*0.9)./cosh(beta2)+c).^0.5; U(11)=Uh; Ub=zeros(U); Ub(1:\$-1)=U(2:\$); qinf=sum((U(1:\$-1)+Ub(1:\$-1))/2*0.1.*h); 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.2; PI=0; delta=1; Umax=ustar./kappa*(log((delta*(h-h)-dhp*h)/(z0hp*h))+2*PI); Umoy=qtot./pf; res=abs(U0-Umoy); ... ... @@ -97,22 +104,12 @@ else Cf=2/(5.1*log10(pf/ks-1)+6)^2; end if(bDbg) then print_r("Cf"); end; //[u res]=fminsearch(@(U0i) find_U0_complet(U0i,pf,C,D,sigma,Cd0,Cf,coeff_contraction,S),U0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); [u res] = fminsearch ( list(find_U0_R0,pf,C,D,sigma,Cd0,Cf,coeff_contraction,S) , U0 ) N= (alpha.*Cf)./(pf./D.*Cd.*C); // N= (alpha.*Cf)./(pf./D.*Cd.*C); // res=abs(U0-(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5); res=abs(U0-u); res=abs(U0-(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5); end endfunction fprintf('res=%f\n',res)
 function [res]= resolve_alpha(alpha,CdCh,R,U0,hstar,hp,C,D,Cd,ustar,choixturb) if 1==1 fprintf('resolve_alpha(alpha=%f,CdCh=%f,R=%f,U0=%f,hstar=%f,hp=%f,C=%f,D=%f,Cd=%f,ustar=%f)\n',alpha,CdCh,R,U0,hstar,hp,C,D,Cd,ustar) end g=9.81; kappa=0.41; L=D*(1/C^0.5-1); beta2=hp.*CdCh./alpha./R; beta=(beta2).^0.5; a1=beta*(hstar-1)/(cosh(beta)); c=1; UhU0=(a1*sinh(beta)+c)^0.5; Uh=UhU0*U0; %choix du modele de turbulence switch choixturb case 1 L1=min(L,0.15*hp); %L1=L;%0.15*hp; % case 2 % % L1=L.*hp./(1/0.41*L+hp); % % case 3 % L1=0.15*min(hp,D./C./Cd);%; % case 4 % if L/hp<1 % L1=(1-L/hp).*L; % if L1<0.05*hp && L/hp>0.8;L1=0.05*hp;end % else % L1=0.41*((1-hp./L))*hp; % if L1<0.05*hp && L/hp<1.2;L1=0.05*hp;end % end end res=abs(alpha*Uh-L1*ustar); fprintf('resolve_alpha res=%f\n',res) end \ No newline at end of file
 function varargout = warning_4a1_1(varargin) % WARNING_4A1_1 MATLAB code for warning_4a1_1.fig % WARNING_4A1_1, by itself, creates a new WARNING_4A1_1 or raises the existing % singleton*. % % H = WARNING_4A1_1 returns the handle to a new WARNING_4A1_1 or the handle to % the existing singleton*. % % WARNING_4A1_1('CALLBACK',hObject,eventData,handles,...) calls the local % function named CALLBACK in WARNING_4A1_1.M with the given input arguments. % % WARNING_4A1_1('Property','Value',...) creates a new WARNING_4A1_1 or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before warning_4a1_1_OpeningFcn gets called. An % unrecognized property name or invalid value makes property application % stop. All inputs are passed to warning_4a1_1_OpeningFcn via varargin. % % *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)". % % See also: GUIDE, GUIDATA, GUIHANDLES % Edit the above text to modify the response to help warning_4a1_1 % Last Modified by GUIDE v2.5 08-Oct-2014 13:48:09 % Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @warning_4a1_1_OpeningFcn, ... 'gui_OutputFcn', @warning_4a1_1_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT % --- Executes just before warning_4a1_1 is made visible. function warning_4a1_1_OpeningFcn(hObject, eventdata, handles, varargin) % This function has no output args, see OutputFcn. % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to warning_4a1_1 (see VARARGIN) % Choose default command line output for warning_4a1_1 handles.output = hObject; % Update handles structure guidata(hObject, handles); % UIWAIT makes warning_4a1_1 wait for user response (see UIRESUME) % uiwait(handles.figure1); % --- Outputs from this function are returned to the command line. function varargout = warning_4a1_1_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % Get default command line output from handles structure varargout{1} = handles.output; % --- Executes on button press in pushbutton1. function pushbutton1_Callback(hObject, eventdata, handles) % hObject handle to pushbutton1 (see GCBO) % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) close
This diff is collapsed.
 function [res]=find_Q_nat(Q,ks,D,h,Cd0,S,L,pf,C,sigma,bDbg) if bDbg then printf('*************************************\n') print_r('Q') end //fonction pour ax=ay maxfun=5000; maxiter=5000; tolfun=1e-16; tolx=1e-16; opt = optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'TolFun',tolfun,'TolX',tolx) g=9.81; kappa=0.41; U0=Q./L./pf; if(bDbg) then print_r("U0"); end; 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 if(bDbg) then print_r("fFr"); 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; if bDbg then print_r('ustar'); end; Cd1=Cd0.*(1+0.4./(pf./D).^2).*fFr; CdCh=Cd1.*C.*htilde; if bDbg then print_r('CdCh'); end; Cf=2./(5.1.*log10(h./ks)+6).^2; U0b=(2*g.*S.*R./(Cd1.*C.*h./D+Cf.*R).*h).^0.5; if bDbg then print_r('U0b'); end; // 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); [alpha fval] = fminsearch(list(resolve_alpha, CdCh,R,U0b,hstar,h,C,D,Cd1,ustar), h/10, opt) if bDbg then print_r('alpha'); end; beta2=h.*CdCh./alpha./R; beta=(beta2).^0.5; a1=beta*(hstar-1)/(cosh(beta)); c=1; UhU0=(a1*sinh(beta)+c)^0.5; Uh=UhU0*U0b; if bDbg then print_r('Uh'); end; 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))); if bDbg then print_r('qsup'); end; //calcul intégrale dans la canopée---- dzinf=0.01; Zinf=(0:dzinf:1); Uinf = U0b.*(beta.*Rh./h.*sinh(beta*Zinf)./cosh(beta)+1).^0.5; Ub=zeros(Uinf); Ub(1:\$-1)=Uinf(2:\$); qinf=sum((Uinf(1:\$-1)+Ub(1:\$-1))/2*dzinf.*h); if bDbg then print_r('qinf'); end; 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 endfunction
 function [res]= resolve_alpha(alpha,CdCh,R,U0,hstar,hp,C,D,Cd,ustar) if bDbg then printf('resolve_alpha(alpha=%f,CdCh=%f,R=%f,U0=%f,hstar=%f,hp=%f,C=%f,D=%f,Cd=%f,ustar=%f)\n',alpha,CdCh,R,U0,hstar,hp,C,D,Cd,ustar) end g=9.81; kappa=0.41; ... ... @@ -41,5 +43,8 @@ Uh=UhU0*U0; res=abs(alpha*Uh-L1*ustar); if bDbg then printf('resolve_alpha res=%f\n',res) end endfunction
 ... ... @@ -11,40 +11,33 @@ function macrorugo_resultComp(z_amont, S, long, Q, L, pf, C, Cd, h, D) Vg=Q/L/pf/(1-C^0.5); Fr=Vg./(g.*pf).^0.5; print_r("Fr"); if Cd==2 Vmax=Vg.*(min(1.6./(1-(Fr.^2)/4),Fr.^(-2/3)));% else Vmax=Vg.*(min(1.2./(1-(Fr.^2)/4),Fr.^(-2/3))); end print_r("Vmax"); P=1000*g*Q/L*S; print_r("P"); if pf/h<1 flowcond = 'emergent' elseif pf/h<1.1 & pf/h>=1 elseif pf/h<1.1 && pf/h>=1 flowcond = 'quasi emergent' elseif pf/h>1.2 else flowcond = 'immerge' end print_r("flowcond"); if pf/h>1.1 q_technique=0.955*(pf/h)^2.282*S^0.466*C^(-0.23)*(9.81*D)^0.5.*h*L; q_technique=0.955*(pf/h)^2.282*S^0.466*C^(-0.23)*(9.81*h)^0.5.*h*L; else coeff_contraction=0.4*Cd+0.7; if Cd==2 q_technique= 0.648*(pf/D)^1.084*S^0.56*C^(-0.456)*(9.81*D)^0.5.*D*L; V_technique=3.35*(pf/D)^0.27*S^0.53*(9.81*D)^0.5; Vmax=Vg.*(min(coeff_contraction/(1-(Fr.^2)/4),Fr.^(-2/3))); else q_technique=0.815*(pf/D)^1.45*S^0.557*C^(-0.456)*(9.81*D)^0.5.*D*L; V_technique=4.54*(pf/D)^0.32*S^0.56*(9.81*D)^0.5; Vmax=Vg.*(min(coeff_contraction/(1-(Fr.^2)/4),Fr.^(-2/3))); end print_r("Vmax"); print_r("V_technique"); end print_r("q_technique"); ... ...
 function macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg) printf("ks=%f\n", ks) printf("D=%f\n", D) printf("k=%f\n", k) printf("Cd0=%f\n", Cd0) printf("S=%f\n", S) printf("B=%f\n", B) printf("h=%f\n", h) printf("C=%f\n", C) maxfun=5000; maxiter=5000; tolfun=1e-24; tolx=1e-24; if bDbg then maxiter=1 end opt = optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'TolFun',tolfun,'TolX',tolx) if Cd0==2; sigma=1; else sigma=%pi/4; end g=9.81 N=0; q0=(2*g.*S.*D.*(1-(sigma*C))/(Cd0.*C.*(1+N))).^0.5*h*B; [Q fVal, exitflag, outputs] = fminsearch(list(find_Q_nat, ks,D,k,Cd0,S,B,h,C,sigma,bDbg), q0, opt); printf("Q=%f fVal=%f\n",Q, fVal); macrorugo_resultComp(z_amont, S, long, Q, B, h, C, Cd0, k, D) endfunction
 ... ... @@ -3,46 +3,39 @@ clear sCurrentPath = get_absolute_file_path("main_macrorugo.sce"); getd(sCurrentPath); bDbg = %f; // Tests parameters ks = 0.01 Cd0=1.5 S = 0.05 C = 0.05; D = 0.5; sigma = 1; B = 1 z_amont = 12.5; long = 6; ks = 0.01 // Rugosité de fond (m) Cd0=1.5 // Forme S = 0.05 // Pente C = 0.13; // Concentration D = 0.5; // Diamètre B = 1 // Largeur z_amont = 12.5; // Cote amont (m) long = 6; // Longueur rampe (m) // ***************************************************************************** printf("\n*** Emergent conditions ***\n") printf("\n*** Emergent conditions Cd=1.5***\n") // ***************************************************************************** h = 0.6 k = 0.8 // Test de find_U0_R0 Cf=2/(5.1*log10(h/ks-1)+6)^2 coeff_contraction=0.4*Cd0+0.7 res = find_U0_R0(2.58, h, C, D, sigma, Cd0, Cf, coeff_contraction, S) printf("find_U0_R0=%f\n",res); // Test de find_Q_nat2 [Q fVal] = fminsearch(list(find_Q_nat2, ks,D,k,Cd0,S,B,h,C,sigma,%f), 0.1); printf("Qemerg=%f fVal=%f\n",Q, fVal); res = find_Q_nat2(Q,ks,D,k,Cd0,S,B,h,C,sigma,%t) printf("find_Q_nat2=%f\n",res); macrorugo_resultComp(z_amont, S, long, Q, B, h, C, Cd0, k, D) k = 0.7 Cd0 = 1.5 macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg) // ***************************************************************************** printf("\n*** Emergent conditions Cd=2***\n") // ***************************************************************************** h = 0.6 k = 0.7 Cd0 = 2 macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg) // ***************************************************************************** printf("\n*** Submerged conditions ***\n") // ***************************************************************************** k = 0.6 k = 0.7 h = 0.8 [Q fVal] = fminsearch(list(find_Q_nat2, ks,D,k,Cd0,S,B,h,C,sigma,%f), 0.1); printf("Qsubmerg=%f fVal=%f\n",Q, fVal); res = find_Q_nat2(Q,ks,D,k,Cd0,S,B,h,C,sigma,%t) printf("find_Q_nat2=%f\n",res); macrorugo_resultComp(z_amont, S, long, Q, B, h, C, Cd0, k, D) Cd0 = 1.5 macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg)
 function print_r(s) e = eval(s) e = evstr(s) if typeof(e) == "string" then pat = "s" else ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!