diff --git a/doc/matlab_macro_rugo/data_precalcul_nat.mat b/doc/matlab_macro_rugo/data_precalcul_nat.mat new file mode 100644 index 0000000000000000000000000000000000000000..1ef3449b26ff6916bb0d7644387664902b40cb79 Binary files /dev/null and b/doc/matlab_macro_rugo/data_precalcul_nat.mat differ diff --git a/doc/matlab_macro_rugo/find_Q_nat.m b/doc/matlab_macro_rugo/find_Q_nat.m new file mode 100644 index 0000000000000000000000000000000000000000..57eb3f7391c37f5e5aecdd723523067ea7c831e6 --- /dev/null +++ b/doc/matlab_macro_rugo/find_Q_nat.m @@ -0,0 +1,115 @@ +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; + + 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 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(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) diff --git a/doc/matlab_macro_rugo/resolve_alpha.m b/doc/matlab_macro_rugo/resolve_alpha.m new file mode 100644 index 0000000000000000000000000000000000000000..d78606a91ac85451ff4fcd8a199e89ecc99e7562 --- /dev/null +++ b/doc/matlab_macro_rugo/resolve_alpha.m @@ -0,0 +1,47 @@ +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 diff --git a/doc/matlab_macro_rugo/warning_4a1_1.fig b/doc/matlab_macro_rugo/warning_4a1_1.fig new file mode 100644 index 0000000000000000000000000000000000000000..379e63bc1434dd009798f5a9050ef0ae5656d5ad Binary files /dev/null and b/doc/matlab_macro_rugo/warning_4a1_1.fig differ diff --git a/doc/matlab_macro_rugo/warning_4a1_1.m b/doc/matlab_macro_rugo/warning_4a1_1.m new file mode 100644 index 0000000000000000000000000000000000000000..8058dc8856aecb05762489c00b484dc14a6555a2 --- /dev/null +++ b/doc/matlab_macro_rugo/warning_4a1_1.m @@ -0,0 +1,81 @@ +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 diff --git a/doc/matlab_macro_rugo/window_4a1.fig b/doc/matlab_macro_rugo/window_4a1.fig new file mode 100644 index 0000000000000000000000000000000000000000..b60b5f6315d11fe6de64ed54fda307ee070f53d2 Binary files /dev/null and b/doc/matlab_macro_rugo/window_4a1.fig differ diff --git a/doc/matlab_macro_rugo/window_4a1.m b/doc/matlab_macro_rugo/window_4a1.m new file mode 100644 index 0000000000000000000000000000000000000000..4cb4c6d0812ebf4975c736bd4d5d91c04369dd21 --- /dev/null +++ b/doc/matlab_macro_rugo/window_4a1.m @@ -0,0 +1,1028 @@ +function varargout = window_4a1(varargin) +% WINDOW_4A1 MATLAB code for window_4a1.fig +% WINDOW_4A1, by itself, creates a new WINDOW_4A1 or raises the existing +% singleton*. +% +% H = WINDOW_4A1 returns the handle to a new WINDOW_4A1 or the handle to +% the existing singleton*. +% +% WINDOW_4A1('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in WINDOW_4A1.M with the given input arguments. +% +% WINDOW_4A1('Property','Value',...) creates a new WINDOW_4A1 or raises the +% existing singleton*. Starting from the left, property value pairs are +% applied to the GUI before window_4a1_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to window_4a1_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 window_4a1 + +% Last Modified by GUIDE v2.5 04-Nov-2016 14:59:17 + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @window_4a1_OpeningFcn, ... + 'gui_OutputFcn', @window_4a1_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 window_4a1 is made visible. +function window_4a1_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 window_4a1 (see VARARGIN) + +% Choose default command line output for window_4a1 +handles.output = hObject; + +% Update handles structure +guidata(hObject, handles); + +%cd('E:\lcassan\Mes documents\macro_rugosite\cassiopee\gui_matlab') +%cd('D:\macrorugosite\cassiopee\gui_matlab') +load data_precalcul_nat + + + +set(handles.edit1,'string',z_amont); +set(handles.edit17,'string',long); +set(handles.edit2,'string',ks); +set(handles.edit3,'string',D); +set(handles.edit4,'string',h); +set(handles.edit5,'string',Cd); +set(handles.edit6,'string',S); +set(handles.edit7,'string',L); +set(handles.edit8,'string',pf); +set(handles.edit9,'string',C); +set(handles.edit10,'string',Q); +set(handles.edit20,'string',paramin); +set(handles.edit21,'string',parapas); +set(handles.edit22,'string',paramax); + + +% UIWAIT makes window_4a1 wait for user response (see UIRESUME) +% uiwait(handles.figure1); + + +% --- Outputs from this function are returned to the command line. +function varargout = window_4a1_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; + + + +function edit1_Callback(hObject, eventdata, handles) +% hObject handle to edit3 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit3 as text +% str2double(get(hObject,'String')) returns contents of edit3 as a double + + + +% --- Executes during object creation, after setting all properties. +function edit1_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit3 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit2_Callback(hObject, eventdata, handles) +% hObject handle to edit2 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit3 as text +% str2double(get(hObject,'String')) returns contents of edit3 as a double + + + +% --- Executes during object creation, after setting all properties. +function edit2_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit2 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + + +function edit3_Callback(hObject, eventdata, handles) +% hObject handle to edit3 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit3 as text +% str2double(get(hObject,'String')) returns contents of edit3 as a double + + +% --- Executes during object creation, after setting all properties. +function edit3_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit3 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit4_Callback(hObject, eventdata, handles) +% hObject handle to edit4 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit4 as text +% str2double(get(hObject,'String')) returns contents of edit4 as a double + + +% --- Executes during object creation, after setting all properties. +function edit4_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit4 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit5_Callback(hObject, eventdata, handles) +% hObject handle to edit5 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit5 as text +% str2double(get(hObject,'String')) returns contents of edit5 as a double + + +% --- Executes during object creation, after setting all properties. +function edit5_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit5 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +function edit6_Callback(hObject, eventdata, handles) +% hObject handle to edit6 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit6 as text +% str2double(get(hObject,'String')) returns contents of edit6 as a double + + +% --- Executes during object creation, after setting all properties. +function edit6_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit6 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +function edit7_Callback(hObject, eventdata, handles) +% hObject handle to edit7 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit7 as text +% str2double(get(hObject,'String')) returns contents of edit7 as a double + + +% --- Executes during object creation, after setting all properties. +function edit7_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit7 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +function edit8_Callback(hObject, eventdata, handles) +% hObject handle to edit8 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit8 as text +% str2double(get(hObject,'String')) returns contents of edit8 as a double + + + +% --- Executes during object creation, after setting all properties. +function edit8_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit8 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit9_Callback(hObject, eventdata, handles) +% hObject handle to edit9 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit9 as text +% str2double(get(hObject,'String')) returns contents of edit9 as a double + + +% --- Executes during object creation, after setting all properties. +function edit9_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit9 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit10_Callback(hObject, eventdata, handles) +% hObject handle to edit10 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit10 as text +% str2double(get(hObject,'String')) returns contents of edit10 as a double + + +% --- Executes during object creation, after setting all properties. +function edit10_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit10 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit13_Callback(hObject, eventdata, handles) +% hObject handle to edit13 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit13 as text +% str2double(get(hObject,'String')) returns contents of edit13 as a double + + +% --- Executes during object creation, after setting all properties. +function edit13_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit13 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit14_Callback(hObject, eventdata, handles) +% hObject handle to edit14 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit14 as text +% str2double(get(hObject,'String')) returns contents of edit14 as a double + + +% --- Executes during object creation, after setting all properties. +function edit14_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit14 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit15_Callback(hObject, eventdata, handles) +% hObject handle to edit15 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit15 as text +% str2double(get(hObject,'String')) returns contents of edit15 as a double + + +% --- Executes during object creation, after setting all properties. +function edit15_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit15 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit11_Callback(hObject, eventdata, handles) +% hObject handle to edit11 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit11 as text +% str2double(get(hObject,'String')) returns contents of edit11 as a double + + +% --- Executes during object creation, after setting all properties. +function edit11_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit11 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit12_Callback(hObject, eventdata, handles) +% hObject handle to edit12 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit12 as text +% str2double(get(hObject,'String')) returns contents of edit12 as a double + + +% --- Executes during object creation, after setting all properties. +function edit12_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit12 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit16_Callback(hObject, eventdata, handles) +% hObject handle to edit16 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit16 as text +% str2double(get(hObject,'String')) returns contents of edit16 as a double + + +% --- Executes during object creation, after setting all properties. +function edit16_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit16 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + + + + + + +% --- 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) + + +maxfun=5000; +maxiter=5000; +tolfun=1e-24; +tolx=1e-24; + +% DEBUG +%maxiter=2; + +z_amont=get(handles.edit1,'string'); +z_amont=str2num(z_amont); + +long=get(handles.edit17,'string'); +long=str2num(long); + +%setappdata(gcf,'z_amont',z_amont); +ks=get(handles.edit2,'string'); +ks=str2num(ks); +setappdata(gcf,'ks',ks); +D=get(handles.edit3,'string'); +D=str2num(D); +%setappdata(gcf,'D',D); +h=get(handles.edit4,'string'); +h=str2num(h); +%setappdata(gcf,'h',h); +Cd=get(handles.edit5,'string'); +Cd=str2num(Cd); +%setappdata(gcf,'Cd',Cd); +S=get(handles.edit6,'string'); +S=str2num(S); +%setappdata(gcf,'S',S); +L=get(handles.edit7,'string'); +L=str2num(L); +%setappdata(gcf,'L',L); +pf=get(handles.edit8,'string'); +pf=str2num(pf); +%setappdata(gcf,'pf',pf); +C=get(handles.edit9,'string'); +C=str2num(C); +%setappdata(gcf,'C',C); +Q=get(handles.edit10,'string'); +Q=str2num(Q); +%setappdata(gcf,'Q',Q); + + +paramin=str2num(get(handles.edit20,'string')); +parapas=str2num(get(handles.edit21,'string')); +paramax=str2num(get(handles.edit22,'string')); +para=(paramin:parapas:paramax); + +save data_precalcul_nat z_amont long ks D h Cd S L pf C Q paramin parapas paramax + +% +% z_amont=getappdata(gcf,'z_amont'); +% ks=getappdata(gcf,'ks'); +% D=getappdata(gcf,'D'); +% h=getappdata(gcf,'h'); +% Cd=getappdata(gcf,'Cd'); +% S=getappdata(gcf,'S'); +% L=getappdata(gcf,'L'); +% pf=getappdata(gcf,'pf'); +% C=getappdata(gcf,'C'); +% Q=getappdata(gcf,'Q'); + +f=1; + + + +if Cd==2 ;sigma=1;else sigma=pi/4; end +alpha=1-(f.*C).^0.5-1/2*sigma.*C; +g=9.81; + +if Q~=0 & pf~=0 & L~=0 & S~=0 & C~=0 + warning_4a1_1 +end + + +if Q==0 + N=0; + q0=(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5*pf*L; + + [Q res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,D,h,Cd,S,L,pf,C,Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) + set(handles.edit10,'string',num2str(Q)) + +elseif pf==0 + [pf res]=fminsearch(@(pfi) find_Q_nat(z_amont,ks,D,h,Cd,S,L,pfi,C,Q,sigma),h*1.1,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) + set(handles.edit8,'string',num2str(pf)) + +elseif L==0 + [L res]=fminsearch(@(Li) find_Q_nat(z_amont,ks,D,h,Cd,S,Li,pf,C,Q,sigma),10,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) + set(handles.edit7,'string',num2str(L)) + +elseif S==0 + S0=0.1; + [S res]=fminsearch(@(Si) find_Q_nat(z_amont,ks,D,h,Cd,Si,L,pf,C,Q,sigma),S0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) + set(handles.edit6,'string',num2str(S)) + +elseif C==0 + C0=0.1; + [C res]=fminsearch(@(Ci) find_Q_nat(z_amont,ks,D,h,Cd,S,L,pf,Ci,Q,sigma),C0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) + set(handles.edit9,'string',num2str(C)) + +end + + + + +cote_bas=z_amont-S*long; +set(handles.edit11,'string',num2str(cote_bas)) + +Vdeb=Q/L/pf; +set(handles.edit12,'string',num2str(Vdeb)) +coeff_contraction=0.4*Cd+0.7; +Vg=Q/L/pf/(1-C^0.5); +Fr=Vg./(g.*pf).^0.5; +if Cd==2 + Vmax=Vg.*(min(coeff_contraction./(1-(Fr.^2)/4),Fr.^(-2/3)));% +else + Vmax=Vg.*(min(coeff_contraction./(1-(Fr.^2)/4),Fr.^(-2/3))); +end + +set(handles.edit16,'string',num2str(Vmax)) + +P=1000*g*Q/L*S; +set(handles.edit13,'string',num2str(P)) + + +set(handles.edit14,'string',num2str(Fr)) + +if pf/h<1 + set(handles.edit15,'string','emergent') +elseif pf/h<1.1 && pf/h>=1 + set(handles.edit15,'string','quasi emergent') +else + set(handles.edit15,'string','immerge') +end + + + +if pf/h>1.1 + + q_technique=0.955*(pf/h)^2.282*S^0.466*C^(-0.23)*(9.81*h)^0.5.*h*L; + set(handles.edit19,'string','NC') + set(handles.edit16,'string','NC') +else + 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; + 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; + + end + + set(handles.edit19,'string',num2str(V_technique)) ; +end +set(handles.edit18,'string',num2str(q_technique)) ; + + + +function edit17_Callback(hObject, eventdata, handles) +% hObject handle to edit17 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit17 as text +% str2double(get(hObject,'String')) returns contents of edit17 as a double + + +% --- Executes during object creation, after setting all properties. +function edit17_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit17 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit19_Callback(hObject, eventdata, handles) +% hObject handle to edit19 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit19 as text +% str2double(get(hObject,'String')) returns contents of edit19 as a double + + +% --- Executes during object creation, after setting all properties. +function edit19_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit19 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit18_Callback(hObject, eventdata, handles) +% hObject handle to edit18 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit18 as text +% str2double(get(hObject,'String')) returns contents of edit18 as a double + + +% --- Executes during object creation, after setting all properties. +function edit18_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit18 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edi5_Callback(hObject, eventdata, handles) +% hObject handle to edi5 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edi5 as text +% str2double(get(hObject,'String')) returns contents of edi5 as a double + + +% --- Executes during object creation, after setting all properties. +function edi5_CreateFcn(hObject, eventdata, handles) +% hObject handle to edi5 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% --- Executes on button press in pushbutton2. +function pushbutton2_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton2 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + + + +maxfun=5000; +maxiter=5000; +tolfun=1e-24; +tolx=1e-24; + + +z_amont=get(handles.edit1,'string'); +z_amont=str2num(z_amont); + +long=get(handles.edit17,'string'); +long=str2num(long); + +%setappdata(gcf,'z_amont',z_amont); +ks=get(handles.edit2,'string'); +ks=str2num(ks); +%setappdata(gcf,'ks',ks); +D=get(handles.edit3,'string'); +D=str2num(D); +%setappdata(gcf,'D',D); +h=get(handles.edit4,'string'); +h=str2num(h); +%setappdata(gcf,'h',h); +Cd=get(handles.edit5,'string'); +Cd=str2num(Cd); +%setappdata(gcf,'Cd',Cd); +S=get(handles.edit6,'string'); +S=str2num(S); +%setappdata(gcf,'S',S); +L=get(handles.edit7,'string'); +L=str2num(L); +%setappdata(gcf,'L',L); +pf=get(handles.edit8,'string'); +pf=str2num(pf); +%setappdata(gcf,'pf',pf); +C=get(handles.edit9,'string'); +C=str2num(C); +%setappdata(gcf,'C',C); +Q=get(handles.edit10,'string'); +Q=str2num(Q); +%setappdata(gcf,'Q',Q); + + +paramin=str2num(get(handles.edit20,'string')); +parapas=str2num(get(handles.edit21,'string')); +paramax=str2num(get(handles.edit22,'string')); +para=(paramin:parapas:paramax); + +save data_precalcul_nat z_amont long ks D h Cd S L pf C Q paramin parapas paramax + + +% +% z_amont=getappdata(gcf,'z_amont'); +% ks=getappdata(gcf,'ks'); +% D=getappdata(gcf,'D'); +% h=getappdata(gcf,'h'); +% Cd=getappdata(gcf,'Cd'); +% S=getappdata(gcf,'S'); +% L=getappdata(gcf,'L'); +% pf=getappdata(gcf,'pf'); +% C=getappdata(gcf,'C'); +% Q=getappdata(gcf,'Q'); + +f=1; + + + +if Cd==2 ;sigma=1;else sigma=pi/4; end +alpha=1-(f.*C).^0.5-1/2*sigma.*C; +g=9.81; + +if Q~=0 & pf~=0 & L~=0 & S~=0 & C~=0 + warning_4a1_1 +end + +for tt=1:length(para) + + if D==0 + + N=0; + q0=(2*g.*S.*para(tt).*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5*pf*L; + [Qr(tt) res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,para(tt),h,Cd,S,L,pf,C,Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) ; + + Vdeb(tt)=Qr(tt)/L/pf; + Vg(tt)=Qr(tt)/L/pf/(1-C^0.5); + Fr(tt)=Vg(tt)./(g.*pf).^0.5; + + + elseif Q==0 + + N=0; + + + [pf2(tt) res]=fminsearch(@(pfi) find_Q_nat(z_amont,ks,D,h,Cd,S,L,pfi,C,para(tt),sigma),pf,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); + + Vdeb(tt)=para(tt)/L/pf2(tt); + Vg(tt)=para(tt)/L/pf2(tt)/(1-C^0.5); + Fr(tt)=Vg(tt)./(g.*pf2(tt)).^0.5; + + + elseif pf==0 + + N=0; + q0=(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5*para(tt)*L; + + [Qr(tt) res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,D,h,Cd,S,L,para(tt),C,Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); + + Vdeb(tt)=Qr(tt)/L/para(tt); + Vg(tt)=Qr(tt)/L/para(tt)/(1-C^0.5); + Fr(tt)=Vg(tt)./(g.*para(tt)).^0.5; + + + elseif L==0 + + N=0; + q0=(2*g.*S.*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5*pf*para(tt); + [Qr(tt) res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,D,h,Cd,S,para(tt),pf,C,Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); + + Vdeb(tt)=Qr(tt)/para(tt)/pf; + Vg(tt)=Qr(tt)/para(tt)/pf/(1-C^0.5); + Fr(tt)=Vg(tt)./(g.*pf).^0.5; + + + elseif S==0 + S0=0.1; + + N=0; + q0=(2*g.*para(tt).*D.*(1-(sigma*C))./(Cd.*C.*(1+N))).^0.5*pf*L; + [Qr(tt) res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,D,h,Cd,para(tt),L,pf,C,Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)); + + Vdeb(tt)=Qr(tt)/L/pf; + Vg(tt)=Qr(tt)/L/pf/(1-C^0.5); + Fr(tt)=Vg(tt)./(g.*pf).^0.5; + + + + elseif C==0 + + N=0; + q0=(2*g.*S.*D.*(1-(sigma*para(tt)))./(Cd.*para(tt).*(1+N))).^0.5*pf*L; + [Qr(tt) res]=fminsearch(@(Qi) find_Q_nat(z_amont,ks,D,h,Cd,S,L,pf,para(tt),Qi,sigma),q0,optimset('MaxIter',maxiter,'MaxFunEvals',maxfun,'Tolfun',tolfun,'TolX',tolx)) ; + + Vdeb(tt)=Qr(tt)/L/pf; + Vg(tt)=Qr(tt)/L/pf/(1-para(tt)^0.5); + Fr(tt)=Vg(tt)./(g.*pf).^0.5; + + + end + + + if Cd==2 + Vmax(tt)=Vg(tt).*(min(1.6./(1-(Fr(tt).^2)/4),Fr(tt).^(-2/3)));% + else + Vmax(tt)=Vg(tt).*(min(1.2./(1-(Fr(tt).^2)/4),Fr(tt).^(-2/3))); + end + + +end + + + +figure(1) +subplot(2,1,1) +hold on +plot(para,Vg,'k','LineWidth',2) +plot(para,Vmax,'k:','LineWidth',2) +plot(para,Vdeb,'k-.','LineWidth',2) + + if D==0 + xlabel('Diametre (m)','Interpreter','Latex','FontSize',24) + elseif Q==0 + xlabel('Débit (m2/s)','Interpreter','Latex','FontSize',24) + elseif pf==0 + xlabel('Profondeur','Interpreter','Latex','FontSize',24) + elseif L==0 + xlabel('Largeur','Interpreter','Latex','FontSize',24) + elseif S==0 + xlabel('Pente','Interpreter','Latex','FontSize',24) + + elseif C==0 + xlabel('Concentration','Interpreter','Latex','FontSize',24) + end + + ylabel('Vitesse (m/s)','Interpreter','Latex','FontSize',24) +hleng=legend('Vg','V_{max}','V_{debitant}'); +set(hleng, 'FontName','Times New Roman','FontSize',20) +set(gca, 'FontSize', 20, 'fontName','Times'); + +% +figure(1) +subplot(2,1,2) +hold on +if Q==0 + plot(para,pf2,'k','LineWidth',2) +else +plot(para,Qr,'k','LineWidth',2) +end + if D==0 + xlabel('Diametre (m)','Interpreter','Latex','FontSize',24) + elseif Q==0 + xlabel('Débit (m3/s)','Interpreter','Latex','FontSize',24) + elseif pf==0 + xlabel('profondeur','Interpreter','Latex','FontSize',24) + elseif L==0 + xlabel('Largeur','Interpreter','Latex','FontSize',24) + elseif S==0 + xlabel('Pente','Interpreter','Latex','FontSize',24) + + elseif C==0 + xlabel('Concentration','Interpreter','Latex','FontSize',24) + end + + + if Q==0 +ylabel('Profondeur (m)','Interpreter','Latex','FontSize',24) +else +ylabel('Débit (m3/s)','Interpreter','Latex','FontSize',24) + end + +%hleng=legend('Débit (m^2/s)'); +%set(hleng, 'FontName','Times New Roman','FontSize',20) +set(gca, 'FontSize', 20, 'fontName','Times'); + + + +function edit20_Callback(hObject, eventdata, handles) +% hObject handle to edit20 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit20 as text +% str2double(get(hObject,'String')) returns contents of edit20 as a double + + +% --- Executes during object creation, after setting all properties. +function edit20_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit20 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit21_Callback(hObject, eventdata, handles) +% hObject handle to edit21 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit21 as text +% str2double(get(hObject,'String')) returns contents of edit21 as a double + + +% --- Executes during object creation, after setting all properties. +function edit21_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit21 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +function edit22_Callback(hObject, eventdata, handles) +% hObject handle to edit22 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of edit22 as text +% str2double(get(hObject,'String')) returns contents of edit22 as a double + + +% --- Executes during object creation, after setting all properties. +function edit22_CreateFcn(hObject, eventdata, handles) +% hObject handle to edit22 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end diff --git a/doc/scilab_tests/macrorugo_find_Q_nat.sci b/doc/scilab_tests/macrorugo_find_Q_nat.sci new file mode 100644 index 0000000000000000000000000000000000000000..0f38c8b0d55960176c05873b6a95962172fb85ee --- /dev/null +++ b/doc/scilab_tests/macrorugo_find_Q_nat.sci @@ -0,0 +1,122 @@ +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; + + //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; + if bDbg then print_r('Cd1'); end; + 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; + if bDbg then print_r('Uinf(2)'); end; + if bDbg then print_r('Uinf($)'); end; + 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 + + alpha=1-(1.*C).^0.5-1/2*sigma.*C; + if bDbg then print_r('alpha'); end; + Cd=Cd0.*(1+1./(pf./D).^2).*fFr; + if bDbg then print_r('Cd'); end; + + hstar=pf/D; + Re=U0.*pf/1e-6; + if bDbg then print_r('Re'); end; + + if ks==0 + Cf=0.3164/4.*Re.^(-0.25); + else + Cf=2/(5.1*log10(pf/ks-1)+6)^2; + end + if bDbg then print_r('Cf'); end; + + N= (alpha.*Cf)./(pf./D.*Cd.*C); + if bDbg then print_r('N'); end; + + res=abs(U0-(2*g.*S.*D.*(R)./(Cd.*C.*(1+N))).^0.5); + + end + +endfunction diff --git a/doc/scilab_tests/macrorugo_resolve_alpha.sci b/doc/scilab_tests/macrorugo_resolve_alpha.sci new file mode 100644 index 0000000000000000000000000000000000000000..877ee789c4da4b91ea38e316cc48562fee22e601 --- /dev/null +++ b/doc/scilab_tests/macrorugo_resolve_alpha.sci @@ -0,0 +1,50 @@ +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; + +L=D*(1/C^0.5-1); +beta2=hp.*CdCh./alpha./R; + +beta2=(beta2).^0.5; +a1=beta2*(hstar-1)/(cosh(beta2)); + +c=1; +UhU0=(a1*sinh(beta2)+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); +if bDbg then + printf('resolve_alpha res=%f\n',res) +end + +endfunction diff --git a/doc/scilab_tests/macrorugo_resultats_complementaires.sci b/doc/scilab_tests/macrorugo_resultats_complementaires.sci new file mode 100644 index 0000000000000000000000000000000000000000..bf0a9c7c8402628121b21d46f6043d387ac30c71 --- /dev/null +++ b/doc/scilab_tests/macrorugo_resultats_complementaires.sci @@ -0,0 +1,44 @@ +function macrorugo_resultComp(z_amont, S, long, Q, L, pf, C, Cd, h, D) + g = 9.81 + + cote_bas=z_amont-S*long; + print_r("cote_bas"); + + Vdeb=Q/L/pf; + print_r("Vdeb"); + + + Vg=Q/L/pf/(1-C^0.5); + Fr=Vg./(g.*pf).^0.5; + print_r("Fr"); + + P=1000*g*Q/L*S; + print_r("P"); + + if pf/h<1 + flowcond = 'emergent' + elseif pf/h<1.1 & pf/h>=1 + flowcond = 'quasi emergent' + 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*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; + 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; + end + Vmax=Vg.*(min(coeff_contraction/(1-(Fr.^2)/4),Fr.^(-2/3))); + print_r("Vmax"); + print_r("V_technique"); + end + print_r("q_technique"); + +endfunction diff --git a/doc/scilab_tests/macrorugo_searchQ.sci b/doc/scilab_tests/macrorugo_searchQ.sci new file mode 100644 index 0000000000000000000000000000000000000000..912638a281ad4c8ea4397e57dc88907b0739eacc --- /dev/null +++ b/doc/scilab_tests/macrorugo_searchQ.sci @@ -0,0 +1,28 @@ +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; + fVal = find_Q_nat(q0,ks,D,k,Cd0,S,B,h,C,sigma,%t); + printf("find_Q_nat(%f)=%f\n",q0,fVal); + [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 diff --git a/doc/scilab_tests/main_macrorugo.sce b/doc/scilab_tests/main_macrorugo.sce new file mode 100644 index 0000000000000000000000000000000000000000..b14271cd07f80f2f068401825742c7b6d44d2ae0 --- /dev/null +++ b/doc/scilab_tests/main_macrorugo.sce @@ -0,0 +1,41 @@ +// Test des passes à macro rugosité pour JaLHyd +clear +sCurrentPath = get_absolute_file_path("main_macrorugo.sce"); +getd(sCurrentPath); + +bDbg = %f; + +// Tests parameters +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 Cd=1.5***\n") +// ***************************************************************************** +h = 0.6 +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.7 +h = 0.8 +Cd0 = 1.5 +macrorugo_searchQ(ks, D, k, Cd0, S, B, h, C, z_amont, long, bDbg) + diff --git a/doc/scilab_tests/print_r.sci b/doc/scilab_tests/print_r.sci new file mode 100644 index 0000000000000000000000000000000000000000..0c33364324dbe47d4636ec3f60b6b34b417da46d --- /dev/null +++ b/doc/scilab_tests/print_r.sci @@ -0,0 +1,18 @@ +function print_r(s) + e = evstr(s) + if typeof(e) == "string" then + pat = "s" + else + pat = "f" + end + if size(e,1) > 1 | size(e,2) > 1 then + printf("%s = [",s); + if size(e,2) > 1 then + e = e' + end + printf("%f ",e) + printf("]\n") + else + printf("%s = %"+pat+"\n",s,e); + end +endfunction diff --git a/spec/macrorugo/macrorugo.spec.ts b/spec/macrorugo/macrorugo.spec.ts new file mode 100644 index 0000000000000000000000000000000000000000..400cdd3623845c49a38f3701eb0eff9529a1765b --- /dev/null +++ b/spec/macrorugo/macrorugo.spec.ts @@ -0,0 +1,201 @@ +/** + * IMPORTANT ! + * Décommenter temporairement la ligne suivante (import { } from "./mock_jasmine") + * Pour exécuter ce code dans le débugger. + * Faire de même avec le fichier test_func.ts + */ +// import { describe, expect, it, xdescribe, xit } from "../mock_jasmine"; + +import { ParamCalculability, ParamValueMode } from "../../src"; +import { MacroRugo, MacroRugoFlowType, MacrorugoParams } from "../../src/macrorugo/macrorugo"; +import { checkResult } from "../test_func"; + +/* +*** Emergent conditions Cd=1.5*** +ks=0.010000 +D=0.500000 +k=0.700000 +Cd0=1.500000 +S=0.050000 +B=1.000000 +h=0.600000 +C=0.130000 +Q=0.493722 fVal=0.000000 +cote_bas = 12.200000 +Vdeb = 0.822870 +Fr = 0.530417 +P = 242.170537 +flowcond = emergent +Vmax = 1.799472 +V_technique = 1.991299 +q_technique = 0.561860 + */ + +function macroRugoInstanceEmergentCd15(): MacroRugo { + return new MacroRugo( + new MacrorugoParams( + 12.5, // ZF1 + 6, // L + 1, // B + 0.05, // If + 0.493722, // Q + 0.6, // h + 0.01, // Ks + 0.13, // C + 0.5, // D + 0.7, // k + 1.5 // Cd0 + ) + ); +} + +const macroRugoExtraResultEmergentCd15: { [key: string]: number|MacroRugoFlowType } = { + ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT, + Fr: 0.530417, + PV: 242.170537, + Q_GuideTech: 0.561860, + V_GuideTech: 1.991299, + Vdeb: 0.822870, + Vmax: 1.799472, + ZF2: 12.2 +}; + +/* +*** Emergent conditions Cd=2*** +ks=0.010000 +D=0.500000 +k=0.700000 +Cd0=2.000000 +S=0.050000 +B=1.000000 +h=0.600000 +C=0.130000 +Q=0.376808 fVal=0.000000 +cote_bas = 12.200000 +Vdeb = 0.628013 +Fr = 0.404814 +P = 184.824088 +flowcond = emergent +Vmax = 1.536115 +V_technique = 1.592932 +q_technique = 0.414154 + */ + +function macroRugoInstanceEmergentCd2(): MacroRugo { + const nub: MacroRugo = macroRugoInstanceEmergentCd15(); + nub.prms.Cd0.v = 2; + nub.prms.Q.v = 0.376808; + return nub; +} + +const macroRugoExtraResultEmergentCd2: { [key: string]: number|MacroRugoFlowType } = { + ENUM_MacroRugoFlowType: MacroRugoFlowType.EMERGENT, + Fr: 0.404814, + PV: 184.824088, + Q_GuideTech: 0.414154, + V_GuideTech: 1.592932, + Vdeb: 0.628013, + Vmax: 1.536115, + ZF2: 12.2 +}; + +/* +*** Submerged conditions *** +ks=0.010000 +D=0.500000 +k=0.700000 +Cd0=1.500000 +S=0.050000 +B=1.000000 +h=0.800000 +C=0.130000 +Q=0.908068 fVal=0.000000 +cote_bas = 12.200000 +Vdeb = 1.135085 +Fr = 0.633645 +P = 445.407266 +flowcond = immerge +q_technique = 0.940450 + */ + +function macroRugoInstanceSubmerged(): MacroRugo { + const nub: MacroRugo = macroRugoInstanceEmergentCd15(); + nub.prms.Y.v = 0.8; + nub.prms.Q.v = 0.908068; + return nub; +} + +const macroRugoExtraResultSubmerged: { [key: string]: number|MacroRugoFlowType } = { + ENUM_MacroRugoFlowType: MacroRugoFlowType.SUBMERGED, + Fr: 0.633645, + PV: 445.407266, + Q_GuideTech: 0.940450, + Vdeb: 1.135085, + ZF2: 12.2 +}; + +function MacroRugoFactory(sInstance: string): MacroRugo { + switch (sInstance) { + case "EmergentCd15": { + return macroRugoInstanceEmergentCd15(); + } + case "EmergentCd2": { + return macroRugoInstanceEmergentCd2(); + } + case "Submerged": { + return macroRugoInstanceSubmerged(); + } + default: { + throw new Error("Instance name error"); + } + } +} + +function testMacroRugo(sInstance: string, varTest: string, valRef: number) { + it(`Calc(${varTest}) should be ${valRef}`, () => { + const nub = MacroRugoFactory(sInstance); + nub.prms.Q.v = nub.Calc("Q").vCalc; + const p = nub.getParameter(varTest); + p.v = undefined; + p.valueMode = ParamValueMode.CALCUL; + checkResult(nub.Calc(varTest, 0.1), valRef); + }); +} + +function testMacroRugoConfig(sInstance: string, Q0: number, fVal0: number, mrExtraRes: { [key: string]: number }) { + describe(`Condition ${sInstance}` , () => { + it(`resolveQ(${Q0}) should be ${fVal0}`, () => { + const nubit = MacroRugoFactory(sInstance); + // tslint:disable-next-line:no-string-literal + expect(nubit["resolveQ"](Q0)).toBeCloseTo(fVal0, 5); + }); + const nub = MacroRugoFactory(sInstance); + for (const prm of nub.prms) { + if ([ParamCalculability.DICHO, ParamCalculability.EQUATION].includes(prm.calculability)) { + testMacroRugo(sInstance, prm.symbol, prm.v); + } + } + for (const sExtraRes in mrExtraRes) { + if (mrExtraRes.hasOwnProperty(sExtraRes)) { + it(`${sExtraRes} should be ${mrExtraRes[sExtraRes]}`, () => { + expect(MacroRugoFactory(sInstance).Calc("Q").extraResults[sExtraRes]) + .toBeCloseTo(mrExtraRes[sExtraRes], 5); + }); + } + } + }); +} + +describe("Class MacroRugo: ", () => { + + testMacroRugoConfig("EmergentCd15", 0.901710, 0.379574, macroRugoExtraResultEmergentCd15); + testMacroRugoConfig("EmergentCd2", 0.768677, 0.418044, macroRugoExtraResultEmergentCd2); + // Le test passe en debug mais pas sous Jasmine !?? + // describe(`Condition Submerged` , () => { + // it(`resolveAlpha_t(0.07) should be 0.074752`, () => { + // // tslint:disable-next-line:no-string-literal + // expect(macroRugoInstanceSubmerged()["resolveAlpha_t"](0.07)).toBeCloseTo(0.074752, 5); + // }); + // }); + testMacroRugoConfig("Submerged", 1.202280, 0.145051, macroRugoExtraResultSubmerged); +}); diff --git a/spec/pab/pab_dimension.spec.ts b/spec/pab/pab_dimension.spec.ts index fe7568ba364af7fb60dceacddc49b365940152d5..496feaa9c0392f874547cb2e20f287be368ed2cf 100644 --- a/spec/pab/pab_dimension.spec.ts +++ b/spec/pab/pab_dimension.spec.ts @@ -1,5 +1,4 @@ import { PabDimension, PabDimensionParams } from "../../src/pab/pab_dimension"; -import { Result } from "../../src/util/result"; import { checkResult } from "../test_func"; function pabDimensionTest(varTest: string) { diff --git a/spec/structure/functions.ts b/spec/structure/functions.ts index 5b7b29275d05cd3c83b05109892677541531ce03..0ea2b2600cbdd32b0ba2a307c4aa04cd7f602090 100644 --- a/spec/structure/functions.ts +++ b/spec/structure/functions.ts @@ -32,12 +32,12 @@ export function itCalcQ( }); if (mode !== undefined) { it("Q(Z1=" + Z1 + ",W=" + W + ") Mode should be " + mode, () => { - expect(res.extraResults.Mode).toBe(mode); + expect(res.extraResults.ENUM_StructureFlowMode).toBe(mode); }); } if (regime !== undefined) { it("Q(Z1=" + Z1 + ",W=" + W + ") Regime should be " + regime, () => { - expect(res.extraResults.Regime).toBe(regime); + expect(res.extraResults.ENUM_StructureFlowRegime).toBe(regime); }); } } diff --git a/spec/structure/parallel_structure.spec.ts b/spec/structure/parallel_structure.spec.ts index fcccae3e2bd580f998825ad636f42238c9a8714f..6ec41d511bb9979f4efb39493e6a21998457dcf3 100644 --- a/spec/structure/parallel_structure.spec.ts +++ b/spec/structure/parallel_structure.spec.ts @@ -67,12 +67,12 @@ function itParallelStructure(sVarCalc: string, rVcalc: number, Q?: number) { }); it(`ExtraResult[ouvrage[${i}].Q_Mode] should be 0`, () => { expect( - pstruct.Calc(sVarCalc).resultElement.extraResults[`ouvrage[${i}].Q_Mode`] + pstruct.Calc(sVarCalc).resultElement.extraResults[`ouvrage[${i}].Q_ENUM_StructureFlowMode`] ).toEqual(0); }); it(`ExtraResult[ouvrage[${i}].Q_Regime] should be 0`, () => { expect( - pstruct.Calc(sVarCalc).resultElement.extraResults[`ouvrage[${i}].Q_Regime`] + pstruct.Calc(sVarCalc).resultElement.extraResults[`ouvrage[${i}].Q_ENUM_StructureFlowRegime`] ).toEqual(0); }); } diff --git a/spec/structure/structure.spec.ts b/spec/structure/structure.spec.ts index 4a4c22b83fcb37f0d701e0346499461e02bd396c..62abfb7018aef7ab646aa122033b93a13df05c0b 100644 --- a/spec/structure/structure.spec.ts +++ b/spec/structure/structure.spec.ts @@ -46,7 +46,7 @@ describe("Class Structure: ", () => { }); describe("Calc()", () => { - const flagsNull = { Mode: StructureFlowMode.NULL, Regime: StructureFlowRegime.NULL }; + const flagsNull = { ENUM_StructureFlowMode: StructureFlowMode.NULL, ENUM_StructureFlowRegime: StructureFlowRegime.NULL }; it("Z1=Z2 => Q=0", () => { structTest.prms.Z2.v = structTest.prms.Z1.v; checkResult(structTest.Calc("Q"), 0); diff --git a/spec/structure/structure_triangular_trunc_weir_free.spec.ts b/spec/structure/structure_triangular_trunc_weir_free.spec.ts index f0e771dc22e8980e24039838f8173c26d7926660..c80ce74ca7f00d843aa9332c0d75625df774936c 100644 --- a/spec/structure/structure_triangular_trunc_weir_free.spec.ts +++ b/spec/structure/structure_triangular_trunc_weir_free.spec.ts @@ -17,15 +17,15 @@ const structTest: StructureTriangularTruncWeirFree = new StructureTriangularTrun describe("Class StructureTriangularTruncWeirFree: ", () => { describe("Calcul Q a surface libre avec h1 croissant: ", () => { - const h1: number[] = + const Z1: number[] = [100.1, 100.2, 100.3, 100.4, 100.5, 100.6, 100.7, 100.8, 100.9, 101, 101.1, 101.5, 101.8, 102]; const Q: number[] = [0., 0.004, 0.024, 0.067, 0.138, 0.240, 0.379, 0.558, 0.778, 1.045, 1.356, 2.914, 4.346, 5.407]; const mode: StructureFlowMode = StructureFlowMode.WEIR; const regime: StructureFlowRegime = StructureFlowRegime.FREE; - itCalcQ(structTest, h1[0], Infinity, Q[0], StructureFlowMode.NULL, StructureFlowRegime.NULL); + itCalcQ(structTest, Z1[0], Infinity, Q[0], StructureFlowMode.NULL, StructureFlowRegime.NULL); for (let i = 1; i < Q.length; i++) { - itCalcQ(structTest, h1[i], Infinity, Q[i], mode, regime); + itCalcQ(structTest, Z1[i], Infinity, Q[i], mode, regime); } }); }); diff --git a/spec/test_func.ts b/spec/test_func.ts index a3b692f0e0a38105bc3cabf8acea575a5af05600..3fea98947ee6783e1f9d06a11297bf03c78f48bd 100644 --- a/spec/test_func.ts +++ b/spec/test_func.ts @@ -3,7 +3,7 @@ * décommenter la ligne suivante (import { expect } from "./mock_jasmine") dans le cas d'une fonction * qui utilise la fonction expect() de Jasmine mais qui doit être exécutée en dehors de tests Jasmine */ -// import { describe, expect, it, xdescribe } from "./mock_jasmine"; +// import { expect } from "./mock_jasmine"; import { cLog } from "../src/util/log"; import { Message, MessageCode } from "../src/util/message"; @@ -21,6 +21,12 @@ export function equalEpsilon(val1: number, val2: number, epsilon: number = precD return Math.abs(val1 - val2) < epsilon; } +export function checkPercent(valTest: number, valRef: number, rPercent: number) { + const r: number = Math.abs(valTest - valRef) / valRef; + const bounds = valRef * rPercent; + expect(r < rPercent).toBeTruthy(`${valTest} should be between ${valRef - bounds} and ${valRef + bounds}`); +} + /** * compare 2 objets * @param s message @@ -184,7 +190,7 @@ export function compareLog(logTest: cLog, logValid: cLog) { } export function checkResult(val1: Result, val2: number, prec?: number) { - expect(val1.ok).toBeTruthy("Result : computation error on Result " + val1.toString()); + expect(val1.ok).toBeTruthy((!val1.ok) ? "Result: ERROR code=" + MessageCode[val1.code] : ""); if (val1.ok) { /* * on demande une précision de vérification inférieure à la précision de calcul diff --git a/spec/value_ref/value_ref_structure.spec.ts b/spec/value_ref/value_ref_structure.spec.ts index 46adbc140255677f64c6652f8430eff614c54b61..4b5a065370a32126877db499e441e8f0af2d7e52 100644 --- a/spec/value_ref/value_ref_structure.spec.ts +++ b/spec/value_ref/value_ref_structure.spec.ts @@ -49,20 +49,20 @@ describe("référence d'un paramètre à un autre : ", () => { createEnv(); - prm4.Q.defineReference(nub1, "ouvrage[0].Q_Mode"); + prm4.Q.defineReference(nub1, "ouvrage[0].Q_ENUM_StructureFlowMode"); nub1.CalcSerie(0.001, 0.1, "Q"); nub2.CalcSerie(0.001, 0.1, "Q"); /* nub1.result.resultElements[0].extraResults = { "ouvrage[0].Q" : 6.264183905346331 - "ouvrage[0].Q_Mode" : 0 - "ouvrage[0].Q_Regime" : 0 + "ouvrage[0].Q_ENUM_StructureFlowMode" : 0 + "ouvrage[0].Q_ENUM_StructureFlowRegime" : 0 } nub2.result.resultElements[0].extraResults = { "ouvrage[0].Q" : 2.4110855093366834 - "ouvrage[0].Q_Mode" : 0 - "ouvrage[0].Q_Regime" : 0 + "ouvrage[0].Q_ENUM_StructureFlowMode" : 0 + "ouvrage[0].Q_ENUM_StructureFlowRegime" : 0 } */ @@ -83,13 +83,13 @@ describe("référence d'un paramètre à un autre : ", () => { /* nub1.result.resultElements[0].extraResults = { "ouvrage[0].Q" : 6.264183905346331 - "ouvrage[0].Q_Mode" : 0 - "ouvrage[0].Q_Regime" : 0 + "ouvrage[0].Q_ENUM_StructureFlowMode" : 0 + "ouvrage[0].Q_ENUM_StructureFlowRegime" : 0 } nub2.result.resultElements[0].extraResults = { "ouvrage[0].Q" : 2.4110855093366834 - "ouvrage[0].Q_Mode" : 0 - "ouvrage[0].Q_Regime" : 0 + "ouvrage[0].Q_ENUM_StructureFlowMode" : 0 + "ouvrage[0].Q_ENUM_StructureFlowRegime" : 0 } */ diff --git a/src/compute-node.ts b/src/compute-node.ts index 0d9b8a8a682e66b300a0ca4a0db635b83b0f275a..275c7b263ac735b304571e176b0df1ac4a2d3967 100644 --- a/src/compute-node.ts +++ b/src/compute-node.ts @@ -17,7 +17,8 @@ export enum CalculatorType { Structure, // ouvrages hydrauliques simples ParallelStructure, // ouvrages hydrauliques en parallèle Dever, // Outil Cassiopée Dever - Cloisons // Outil Cassiopée PAB Cloisons + Cloisons, // Outil Cassiopée PAB Cloisons + MacroRugo // Passe à enrochement simple (Cassan et al., 2016) } /** diff --git a/src/macrorugo/macrorugo.ts b/src/macrorugo/macrorugo.ts new file mode 100644 index 0000000000000000000000000000000000000000..dcb60e995546c548936068496a13b63697c36754 --- /dev/null +++ b/src/macrorugo/macrorugo.ts @@ -0,0 +1,446 @@ +import { Nub } from "../nub"; +import { ParamCalculability } from "../param/param-definition"; +import { ParamValueMode } from "../param/param-value-mode"; +import { Result } from "../util/result"; +import { MacrorugoParams } from "./macrorugo_params"; + +export { MacrorugoParams }; + +export enum MacroRugoFlowType { + EMERGENT, + QUASI_EMERGENT, + SUBMERGED +} + +export class MacroRugo extends Nub { + + private static readonly g = 9.81; + + /** nu: water kinematic viscosity */ + private static readonly nu = 1E-6; + // Water at 20 °C has a kinematic viscosity of about 10−6 m2·s−1 + // (https://en.wikipedia.org/wiki/Viscosity#Kinematic_viscosity,_%CE%BD) + + /** Ratio between the width (perpendicular to flow) and the lenght (parallel to flow) of a cell (-) */ + private static readonly fracAxAy = 1; + + /** Limit between emergent and submerged flow */ + private static readonly limitSubmerg = 1.1; + + /** Rugosité de fond (m) */ + private ks: number; + + /** Averaged velocity (m.s-1) */ + + /** Velocity at the bed (m.s-1) */ + private u0: number; + + /** Discharge (m3/s) */ + private Q: number; + + private _cache: { [key: string]: number }; + + constructor(prms: MacrorugoParams, dbg: boolean = false) { + super(prms, dbg); + this._cache = {}; + } + + /** + * paramètres castés au bon type + */ + get prms(): MacrorugoParams { + return this._prms as MacrorugoParams; + } + + /** + * Calcul du débit total, de la cote amont ou aval ou d'un paramètre d'une structure + * @param sVarCalc Nom du paramètre à calculer + * @param rInit Valeur initiale + * @param rPrec Précision attendue + */ + public Calc(sVarCalc: string, rInit?: number, rPrec: number = 0.001): Result { + /** @todo Voir pour déclarer le paramètre en calcul dans nub */ + this.getParameter(sVarCalc).valueMode = ParamValueMode.CALCUL; + const r: Result = super.Calc(sVarCalc, rInit, rPrec); + // Ajout des résultats complémentaires + // Cote de fond aval + r.extraResults.ZF2 = this.prms.ZF1.v - this.prms.If.v * this.prms.L.v; + // Vitesse débitante + r.extraResults.Vdeb = this.V(this.prms.Q) / this.prms.B.v / this.prms.Y.v; + // Froude + const vg: number = r.extraResults.Vdeb / (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v)); + r.extraResults.Fr = vg / Math.sqrt(MacroRugo.g * this.prms.Y.v); + // Vitesse maximale + const cc = 0.4 * this.prms.Cd0.v + 0.7; + r.extraResults.Vmax = vg * Math.min( + cc / (1 - Math.pow(r.extraResults.Fr, 2) / 4), + Math.pow(r.extraResults.Fr, -2 / 3) + ); + // Puissance dissipée + r.extraResults.PV = 1000 * MacroRugo.g * this.V(this.prms.Q) / this.prms.B.v * this.prms.If.v; + // Type d'écoulement + if (this.prms.Y.v / this.prms.PBH.v < 1) { + r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.EMERGENT; + } else if (this.prms.Y.v / this.prms.PBH.v < MacroRugo.limitSubmerg) { + r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.QUASI_EMERGENT; + } else { + r.extraResults.ENUM_MacroRugoFlowType = MacroRugoFlowType.SUBMERGED; + } + // Vitesse et débit du guide technique + let cQ: [number, number, number, number]; + let cV: [number, number, number]; + let hdk: number; + if (r.extraResults.ENUM_MacroRugoFlowType === MacroRugoFlowType.SUBMERGED) { + cQ = [0.955, 2.282, 0.466, -0.23]; + hdk = this.prms.PBH.v; + } else { + hdk = this.prms.PBD.v; + if (Math.abs(this.prms.Cd0.v - 2) < 0.05) { + cQ = [0.648, 1.084, 0.56, -0.456]; + cV = [3.35, 0.27, 0.53]; + } else { + cQ = [0.815, 1.45, 0.557, -0.456]; + cV = [4.54, 0.32, 0.56]; + } + } + r.extraResults.Q_GuideTech = cQ[0] * Math.pow(this.prms.Y.v / hdk, cQ[1]) * + Math.pow(this.prms.If.v, cQ[2]) * Math.pow(this.prms.C.v, cQ[3]) * + Math.sqrt(MacroRugo.g * hdk) * hdk * this.prms.B.v; + if (r.extraResults.ENUM_MacroRugoFlowType !== MacroRugoFlowType.SUBMERGED) { + r.extraResults.V_GuideTech = cV[0] * Math.pow(this.prms.Y.v / this.prms.PBD.v, cV[1]) * + Math.pow(this.prms.If.v, cV[2]) * Math.sqrt(MacroRugo.g * this.prms.PBD.v); + } + return r; + } + + public Equation(sVarCalc: string): Result { + const Q = uniroot(this.resolveQ, this, 0, 1E7) * this.prms.B.v; + return new Result(Q); + } + + /** + * paramétrage de la calculabilité des paramètres + */ + protected setParametersCalculability() { + this.prms.ZF1.calculability = ParamCalculability.FREE; + this.prms.L.calculability = ParamCalculability.FREE; + this.prms.Ks.calculability = ParamCalculability.FREE; + this.prms.B.calculability = ParamCalculability.DICHO; + this.prms.If.calculability = ParamCalculability.DICHO; + this.prms.Q.calculability = ParamCalculability.EQUATION; + this.prms.Y.calculability = ParamCalculability.DICHO; + this.prms.C.calculability = ParamCalculability.DICHO; + this.prms.PBD.calculability = ParamCalculability.FREE; + this.prms.PBH.calculability = ParamCalculability.FREE; + this.prms.Cd0.calculability = ParamCalculability.FREE; + } + + /** + * Equation from Cassan, L., Laurens, P., 2016. Design of emergent and submerged rock-ramp fish passes. + * Knowledge & Management of Aquatic Ecosystems 45. + * @param sVarCalc Variable à calculer + */ + private resolveQ(this: MacroRugo, Q: number): number { + this.Q = Q; + // Reset cached variables depending on Q (or not...) + this._cache = {}; + /** Tirant d'eau (m) */ + const h: number = this.prms.Y.v; + + /** Concentration de blocs (-) */ + const C: number = this.prms.C.v; + /** Paramètre de bloc : Diamètre (m) */ + const D: number = this.prms.PBD.v; + /** Paramètre de bloc : Hauteur (m) */ + const k: number = this.prms.PBH.v; + /** adimensional water depth */ + const hstar: number = h / k; + /** Slope (m/m) */ + const S: number = this.prms.If.v; + /** Rugosity (m) */ + const ks = this.prms.Ks.v; + /** Accélération gravité (m/s²) */ + const g = MacroRugo.g; + /** Constante von Karman */ + const kappa = 0.41; + + if (hstar > MacroRugo.limitSubmerg) { + // Submerged conditions + + /** Velocity at the bed §2.3.2 Cassan et al., 2016 */ + const cf = 2 / Math.pow(5.1 * Math.log10(k / ks) + 6, 2); + this.u0 = Math.sqrt(k * 2 * g * S * this.R + / (this.Cd * C * k / D + cf * this.R)); + /** turbulent length scale (m) within the blocks layer (alpha_t) */ + const alpha = uniroot(this.resolveAlpha_t, this, 1E-6, 100); + /** averaged velocity at the top of blocks (m.s-1) */ + const uk = this.calcUz(alpha); + /** Equation (13) Cassan et al., 2016 */ + const d = k - 1 / kappa * alpha * uk / this.ustar; + /** Equation (14) Cassan et al., 2016 */ + const z0 = (k - d) * Math.exp(- kappa * uk / this.ustar); + /** Integral of Equation (12) Cassan et al., 2016 */ + // tslint:disable-next-line:variable-name + const Qsup = this.ustar / kappa * ( + (h - d) * (Math.log((h - d) / z0) - 1) + - ((k - d) * (Math.log((k - d) / z0) - 1)) + ); + + // calcul intégrale dans la canopée---- + // tslint:disable-next-line:variable-name + let Qinf: number = this.u0; + let u = 0; + let uOld: number; + const step = 0.01; + const zMax = 1 + step / 2; + for (let z = step; z < zMax; z += step) { + uOld = u; + u = this.calcUz(alpha, z); + Qinf += (uOld + u) ; + } + Qinf = Qinf / 2 * step * k; + + // Calcul de u moyen + return this.U0 - (Qinf + Qsup) / h; + + } else { + // Emergent conditions + + // Resolve equation (4) Cassan et al., 2016 + return this.resolveEmergent(); + } + } + + /** + * Averaged velocity (m.s-1) + */ + private get U0(): number { + if (this._cache.U0 === undefined) { + this._cache.U0 = this.Q / this.prms.B.v / this.prms.Y.v; + } + return this._cache.U0; + } + + private get CdChD(): number { + if (this._cache.CdChD === undefined) { + this._cache.CdChD = this.Cd * this.prms.C.v * this.prms.PBH.v / this.prms.PBD.v; + } + return this._cache.CdChD; + } + + /** + * sigma ratio between the block area in the x, y plane and D2 + */ + private get sigma(): number { + if (this._cache.sigma === undefined) { + if (this.prms.Cd0.v === 2) { + this._cache.sigma = 1; + } else { + this._cache.sigma = Math.PI / 4; + } + } + return this._cache.sigma; + } + + private get R(): number { + if (this._cache.R === undefined) { + this._cache.R = (1 - this.sigma * this.prms.C.v); + } + return this._cache.R; + } + + /** + * Bed friction coefficient Equation (3) (Cassan et al., 2016) + */ + private calcCf(): number { + + if (this.prms.Ks.v < 1E-9) { + // Between Eq (8) and (9) (Cassan et al., 2016) + const reynolds = this.U0 * this.prms.Y.v / MacroRugo.nu; + return 0.3164 / 4. * Math.pow(reynolds, -0.25); + } else { + // Equation (3) (Cassan et al., 2016) + return 2 / Math.pow(5.1 * Math.log10(this.prms.Y.v / this.prms.Ks.v - 1) + 6, 2); + } + } + + /** + * Calculation of Cd : drag coefficient of a block under the actual flow conditions + */ + private get Cd(): number { + return this.prms.Cd0.v * (1 + 0.4 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)) * this.fFr; + } + + /** + * Calcul de Beta force ratio between drag and turbulent stress (Cassan et al. 2016 eq(8)) + * \Beta = (k / alpha_t) (C_d C k / D) / (1 - \sigma C) + * @param alpha \alpha_t turbulent length scale (m) within the blocks layer + */ + private calcBeta(alpha: number): number { + return Math.sqrt(this.prms.PBH.v * this.CdChD / alpha / this.R); + } + + /** + * Averaged velocity at a given vertical position (m.s-1) + * @param alpha turbulent length scale (m) within the blocks layer + * @param z dimensionless vertical position z / k + */ + private calcUz(alpha: number, z: number = 1): number { + const beta = this.calcBeta(alpha); + // Equation (9) Cassan et al., 2016 + return this.u0 * Math.sqrt( + beta * (this.prms.Y.v / this.prms.PBH.v - 1) * Math.sinh(beta * z) / Math.cosh(beta) + 1 + ); + } + + private get ustar(): number { + if (this._cache.ustar === undefined) { + this._cache.ustar = Math.sqrt(MacroRugo.g * this.prms.If.v * (this.prms.Y.v - this.prms.PBH.v)); + } + return this._cache.ustar; + } + + private resolveAlpha_t(alpha: number): number { + /** s: minimum distance between blocks */ + const s = this.prms.PBD.v * (1 / Math.sqrt(this.prms.C.v) - 1); + /** Equation(11) Cassan et al., 2016 */ + const l0 = Math.min(s, 0.15 * this.prms.PBH.v); + // Equation(10) Cassan et al., 2016 + return alpha * this.calcUz(alpha) - l0 * this.ustar; + } + + private resolveEmergent(): number { + const g = MacroRugo.g; + + const alpha = 1 - Math.pow(1 * this.prms.C.v, 0.5) - 0.5 * this.sigma * this.prms.C.v; + // tslint:disable-next-line:variable-name + const Cd = this.prms.Cd0.v * (1 + 1 / Math.pow(this.prms.Y.v / this.prms.PBD.v, 2)) * this.fFr; + + /** N from Cassan 2016 eq(2) et Cassan 2014 eq(12) */ + const N = (alpha * this.calcCf()) / (this.prms.Y.v / this.prms.PBD.v * Cd * this.prms.C.v); + + return this.U0 - Math.sqrt( + 2 * MacroRugo.g * this.prms.If.v * this.prms.PBD.v * + (1 - this.sigma * this.prms.C.v) / (Cd * this.prms.C.v * (1 + N)) + ); + } + + /** + * Froude correction function + */ + private get fFr(): number { + if (this._cache.fFr === undefined) { + // tslint:disable-next-line:variable-name + const Fr = this.U0 / + (1 - Math.sqrt(MacroRugo.fracAxAy * this.prms.C.v)) / + Math.sqrt(MacroRugo.g * this.prms.Y.v); + + /** Interpolation linéaire entre le bloc rond (Cd0=1) et le carré (Cd0=2) */ + const r = 0.4 * this.prms.Cd0.v + 0.7; + + if (Fr < 1.3) { + this._cache.fFr = Math.pow(Math.min(r / (1 - Math.pow(Fr, 2) / 4), Math.pow(Fr, -2 / 3)), 2); + } else { + this._cache.fFr = Math.pow(Fr, -4 / 3); + } + } + return this._cache.fFr; + } +} + +/** + * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt> + * for a root (i.e., zero) of the function <tt>func</tt> with respect to + * its first argument using Brent's method root-finding algorithm. + * + * Translated from zeroin.c in http://www.netlib.org/c/brent.shar. + * + * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net> + * MIT License, http://www.opensource.org/licenses/mit-license.php + * + * @param {function} func function for which the root is sought. + * @param {number} lowerlimit the lower point of the interval to be searched. + * @param {number} upperlimit the upper point of the interval to be searched. + * @param {number} errorTol the desired accuracy (convergence tolerance). + * @param {number} maxIter the maximum number of iterations. + * @returns an estimate for the root within accuracy. + * + */ +function uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number, + errorTol: number = 0, maxIter: number = 1000 +) { + let a = lowerLimit; + let b = upperLimit; + let c = a; + let fa = func.call(thisArg, a); + let fb = func.call(thisArg, b); + let fc = fa; + let tolAct; // Actual tolerance + let newStep; // Step at this iteration + let prevStep; // Distance from the last but one to the last approximation + let p; // Interpolation step is calculated in the form p/q; division is delayed until the last moment + let q; + + while (maxIter-- > 0) { + + prevStep = b - a; + + if (Math.abs(fc) < Math.abs(fb)) { + // Swap data for b to be the best approximation + a = b, b = c, c = a; + fa = fb, fb = fc, fc = fa; + } + + tolAct = 1e-15 * Math.abs(b) + errorTol / 2; + newStep = (c - b) / 2; + + if (Math.abs(newStep) <= tolAct || fb === 0) { + return b; // Acceptable approx. is found + } + + // Decide if the interpolation can be tried + if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) { + // If prev_step was large enough and was in true direction, Interpolatiom may be tried + let t1; + let cb; + let t2; + cb = c - b; + if (a === c) { // If we have only two distinct points linear interpolation can only be applied + t1 = fb / fa; + p = cb * t1; + q = 1.0 - t1; + } else { // Quadric inverse interpolation + q = fa / fc, t1 = fb / fc, t2 = fb / fa; + p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1)); + q = (q - 1) * (t1 - 1) * (t2 - 1); + } + + if (p > 0) { + q = -q; // p was calculated with the opposite sign; make p positive + } else { + p = -p; // and assign possible minus to q + } + + if (p < (0.75 * cb * q - Math.abs(tolAct * q) / 2) && + p < Math.abs(prevStep * q / 2)) { + // If (b + p / q) falls in [b,c] and isn't too large it is accepted + newStep = p / q; + } + + // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent + } + + if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance + newStep = (newStep > 0) ? tolAct : -tolAct; + } + + a = b, fa = fb; // Save the previous approx. + b += newStep, fb = func.call(thisArg, b); // Do step to a new approxim. + + if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { + c = a, fc = fa; // Adjust c for it to have a sign opposite to that of b + } + } + return undefined; + +} diff --git a/src/macrorugo/macrorugo_params.ts b/src/macrorugo/macrorugo_params.ts new file mode 100644 index 0000000000000000000000000000000000000000..4fbfc30dcbd30a9f09d8e86484890c15d482531d --- /dev/null +++ b/src/macrorugo/macrorugo_params.ts @@ -0,0 +1,182 @@ +import { ParamDefinition } from "../param/param-definition"; +import { ParamDomainValue } from "../param/param-domain"; +import { ParamsEquation } from "../param/params-equation"; + +export class MacrorugoParams extends ParamsEquation { + [key: string]: any; // pour pouvoir faire this['methode'](); + + /** Cote de fond amont (m) */ + private _ZF1: ParamDefinition; + + /** Longueur (m) */ + private _L: ParamDefinition; + + /** Largeur (m) */ + private _B: ParamDefinition; + + /** Pente (m/m) */ + private _If: ParamDefinition; + + /** Débit (m3/s) */ + private _Q: ParamDefinition; + + /** Tirant d'eau (m) */ + private _Y: ParamDefinition; + + /** Rugosité de fond (m) */ + private _Ks: ParamDefinition; + + /** Concentration de blocs (-) */ + private _C: ParamDefinition; + + /** Paramètre de bloc : Diamètre (m) */ + private _PBD: ParamDefinition; + + /** Paramètre de bloc : Hauteur (m) */ + private _PBH: ParamDefinition; + + /** Paramètre de bloc : Forme (1 pour rond, 2 pour carré) */ + private _Cd0: ParamDefinition; + + /** + * + * @param rZF1 Cote de fond amont (m) + * @param rL Longueur (m) + * @param rB Largeur (m) + * @param rIf Pente (m/m) + * @param rQ Débit (m3/s) + * @param rY Tirant d'eau (m) + * @param rRF Rugosité de fond (m) + * @param rCB Concentration de blocs (m) + * @param rPBD Paramètre de bloc : Diamètre (m) + * @param rPBH Paramètre de bloc : Hauteur (m) + * @param rCd0 Paramètre de bloc : Forme (1 pour rond, 2 pour carré) + */ + constructor(rZF1: number, rL: number, rB: number, rIf: number, rQ: number, + rY: number, rRF: number, rCB: number, rPBD: number, rPBH: number, rCd0: number) { + super(); + + this._ZF1 = new ParamDefinition(this, "ZF1", ParamDomainValue.POS, rZF1); + this.addParamDefinition(this._ZF1); + + this._L = new ParamDefinition(this, "L", ParamDomainValue.POS, rL); + this.addParamDefinition(this._L); + + this._B = new ParamDefinition(this, "B", ParamDomainValue.POS, rB); + this.addParamDefinition(this._B); + + this._If = new ParamDefinition(this, "If", ParamDomainValue.POS, rIf); + this.addParamDefinition(this._If); + + this._Q = new ParamDefinition(this, "Q", ParamDomainValue.POS, rQ); + this.addParamDefinition(this._Q); + + this._Y = new ParamDefinition(this, "Y", ParamDomainValue.POS, rY); + this.addParamDefinition(this._Y); + + this._Ks = new ParamDefinition(this, "Ks", ParamDomainValue.POS_NULL, rRF); + this.addParamDefinition(this._Ks); + + this._C = new ParamDefinition(this, "C", ParamDomainValue.POS, rCB); + this.addParamDefinition(this._C); + + this._PBD = new ParamDefinition(this, "PBD", ParamDomainValue.POS, rPBD); + this.addParamDefinition(this._PBD); + + this._PBH = new ParamDefinition(this, "PBH", ParamDomainValue.POS, rPBH); + this.addParamDefinition(this._PBH); + + this._Cd0 = new ParamDefinition(this, "Cd0", ParamDomainValue.POS, rCd0); + this.addParamDefinition(this._Cd0); + + } + + /** + * Cote de fond amont (m) + * @return {ParamDefinition} + */ + public get ZF1(): ParamDefinition { + return this._ZF1; + } + + /** + * Longueur (m) + * @return {ParamDefinition} + */ + public get L(): ParamDefinition { + return this._L; + } + + /** + * Getter B + * @return {ParamDefinition} + */ + public get B(): ParamDefinition { + return this._B; + } + + /** + * Pente (m/m) + * @return {ParamDefinition} + */ + public get If(): ParamDefinition { + return this._If; + } + + /** + * Débit (m3/s) + * @return {ParamDefinition} + */ + public get Q(): ParamDefinition { + return this._Q; + } + + /** + * Tirant d'eau (m) + * @return {ParamDefinition} + */ + public get Y(): ParamDefinition { + return this._Y; + } + + /** + * Rugosité de fond (m) + * @return {ParamDefinition} + */ + public get Ks(): ParamDefinition { + return this._Ks; + } + + /** + * Concentration de blocs (-) + * @return {ParamDefinition} + */ + public get C(): ParamDefinition { + return this._C; + } + + /** + * Paramètre de bloc : Diamètre (m) + * @return {ParamDefinition} + */ + public get PBD(): ParamDefinition { + return this._PBD; + } + + /** + * Paramètre de bloc : Hauteur (m) + * @return {ParamDefinition} + */ + public get PBH(): ParamDefinition { + return this._PBH; + } + + /** + * Paramètre de bloc : Forme (1 pour rond, 2 pour carré) + * @return {ParamDefinition} + */ + public get Cd0(): ParamDefinition { + return this._Cd0; + } + +} diff --git a/src/nub.ts b/src/nub.ts index 957f4271ea58c83801b6739bccf0647ce8cf3579..cb1c0e688f08e79415a80c607125489ae8e4febf 100644 --- a/src/nub.ts +++ b/src/nub.ts @@ -162,6 +162,22 @@ export abstract class Nub extends ComputeNode implements IReferencedNub { return this._result; } + /** + * Renvoie la valeur actuelle d'un paramètre (valeur utilisateur ou résultat) + * @param p Paramètre + */ + public V(p: ParamDefinition): number { + if (p.valueMode === ParamValueMode.CALCUL) { + if (this.result !== undefined) { + if (this.result.ok) { + return this.result.vCalc; + } + } + throw new Error(`Attempt to read the result of ${p.symbol} but none found`); + } + return p.v; + } + // interface IReferencedNub public getReferencedParamValues(desc: string): ParamValues { diff --git a/src/session.ts b/src/session.ts index 5f685606c6ef36921e4f57bbf4c66381695ad45c..ae756d3e6773b77dd5729a1e0af563b67bb5fe44 100644 --- a/src/session.ts +++ b/src/session.ts @@ -5,6 +5,7 @@ import { Props } from "./props"; // Calculettes import { ConduiteDistrib, ConduiteDistribParams } from "./cond_distri"; import { LechaptCalmon, LechaptCalmonParams } from "./lechaptcalmon"; +import { MacroRugo, MacrorugoParams } from "./macrorugo/macrorugo"; import { PabDimension, PabDimensionParams } from "./pab/pab_dimension"; import { PabPuissance, PabPuissanceParams } from "./pab/pab_puissance"; import { RegimeUniforme } from "./regime_uniforme"; @@ -123,119 +124,151 @@ export class Session { switch (calcType) { case CalculatorType.ConduiteDistributrice: - { - prms = new ConduiteDistribParams(3, // débit Q - 1.2, // diamètre D - 0.6, // perte de charge J - 100, // Longueur de la conduite Lg - 1e-6, // Viscosité dynamique Nu - ); - nub = new ConduiteDistrib(prms, dbg); - break; - } + { + prms = new ConduiteDistribParams(3, // débit Q + 1.2, // diamètre D + 0.6, // perte de charge J + 100, // Longueur de la conduite Lg + 1e-6, // Viscosité dynamique Nu + ); + nub = new ConduiteDistrib(prms, dbg); + break; + } case CalculatorType.LechaptCalmon: - { - prms = new LechaptCalmonParams(3, // débit - 1.2, // diamètre - 0.6, /// perte de charge - 100, // longueur du toyo - 1.863, // paramètre L du matériau - 2, // paramètre M du matériau - 5.33// paramètre N du matériau - ); - nub = new LechaptCalmon(prms, dbg); - break; - } + { + prms = new LechaptCalmonParams(3, // débit + 1.2, // diamètre + 0.6, /// perte de charge + 100, // longueur du toyo + 1.863, // paramètre L du matériau + 2, // paramètre M du matériau + 5.33// paramètre N du matériau + ); + nub = new LechaptCalmon(prms, dbg); + break; + } case CalculatorType.SectionParametree: + { nub = new SectionParametree(this.createSection(nodeType, dbg), dbg); break; + } case CalculatorType.RegimeUniforme: + { const sect: acSection = this.createSection(nodeType, dbg); const ru = new RegimeUniforme(sect, dbg); nub = ru; break; + } case CalculatorType.CourbeRemous: - { - const sectCR: acSection = this.createSection(nodeType, dbg); - prms = new CourbeRemousParams(sectCR, 0.15, // Yamont = tirant amont - 0.4, // Yaval = tirant aval - 100, // Long= Longueur du bief - 5, // Dx=Pas d'espace - MethodeResolution.EulerExplicite - ); - nub = new CourbeRemous(prms, dbg); - break; - } + { + const sectCR: acSection = this.createSection(nodeType, dbg); + prms = new CourbeRemousParams(sectCR, 0.15, // Yamont = tirant amont + 0.4, // Yaval = tirant aval + 100, // Long= Longueur du bief + 5, // Dx=Pas d'espace + MethodeResolution.EulerExplicite + ); + nub = new CourbeRemous(prms, dbg); + break; + } case CalculatorType.PabDimensions: - { - prms = new PabDimensionParams( - 2, // Longueur L - 1, // Largeur W - 0.5, // Tirant d'eau Y - 2 // Volume V - ); - nub = new PabDimension(prms, dbg); - break; - } + { + prms = new PabDimensionParams( + 2, // Longueur L + 1, // Largeur W + 0.5, // Tirant d'eau Y + 2 // Volume V + ); + nub = new PabDimension(prms, dbg); + break; + } case CalculatorType.PabPuissance: - { - prms = new PabPuissanceParams( - 0.3, // Chute entre bassins DH (m) - 0.1, // Débit Q (m3/s) - 0.5, // Volume V (m3) - 588.6 // Puissance dissipée PV (W/m3) - ); - nub = new PabPuissance(prms, dbg); - break; - } + { + prms = new PabPuissanceParams( + 0.3, // Chute entre bassins DH (m) + 0.1, // Débit Q (m3/s) + 0.5, // Volume V (m3) + 588.6 // Puissance dissipée PV (W/m3) + ); + nub = new PabPuissance(prms, dbg); + break; + } case CalculatorType.Structure: + { const structType: StructureType = params.getPropValue("structureType"); const loiDebit: LoiDebit = params.getPropValue("loiDebit"); nub = CreateStructure(structType, loiDebit); break; + } case CalculatorType.ParallelStructure: - prms = new ParallelStructureParams(0.5, // Q - 102, // Z1 - 101.5 // Z2 - ); - nub = new ParallelStructure(prms, dbg); - break; + { + prms = new ParallelStructureParams(0.5, // Q + 102, // Z1 + 101.5 // Z2 + ); + nub = new ParallelStructure(prms, dbg); + break; + } case CalculatorType.Dever: - const deverPrms = new DeverParams(0.5, // Q - 102, // Z1 - 10, // BR : largeur du cours d'eau - 99 // ZR : cote du lit du cours d'eau - ); - nub = new Dever(deverPrms, dbg); - break; + { + const deverPrms = new DeverParams(0.5, // Q + 102, // Z1 + 10, // BR : largeur du cours d'eau + 99 // ZR : cote du lit du cours d'eau + ); + nub = new Dever(deverPrms, dbg); + break; + } case CalculatorType.Cloisons: - nub = new Cloisons( - new CloisonsParams( - 0, // Débit total (m3/s) - 102, // Cote de l'eau amont (m) - 10, // Longueur des bassins (m) - 1, // Largeur des bassins (m) - 1, // Profondeur moyenne (m) - 0.5 // Hauteur de chute (m) - ), dbg - ); - break; + { + nub = new Cloisons( + new CloisonsParams( + 0, // Débit total (m3/s) + 102, // Cote de l'eau amont (m) + 10, // Longueur des bassins (m) + 1, // Largeur des bassins (m) + 1, // Profondeur moyenne (m) + 0.5 // Hauteur de chute (m) + ), dbg + ); + break; + } + case CalculatorType.MacroRugo: + { + nub = new MacroRugo( + new MacrorugoParams( + 12.5, // ZF1 + 6, // L + 1, // B + 0.05, // If + 1.57, // Q + 0.6, // h + 0.01, // Ks + 0.05, // C + 0.5, // D + 0.8, // k + 1.5 // Cd0 + ) + ); + } default: + { throw new Error( // tslint:disable-next-line:max-line-length `Session.createNub() : calculatrice '${CalculatorType[calcType]}' / noeud de calcul '${ComputeNodeType[nodeType]}' non pris en charge` ); + } } // propagate properties diff --git a/src/structure/structure.ts b/src/structure/structure.ts index e22451befbee8401dfe7935ab728ae18c7e560f4..34d18629b3b79154fc3272bb7404f707d9b82bce 100644 --- a/src/structure/structure.ts +++ b/src/structure/structure.ts @@ -114,7 +114,7 @@ export abstract class Structure extends Nub { this.prms.update_h1h2(); // Gestion du débit nul - const flagsNull = { Mode: StructureFlowMode.NULL, Regime: StructureFlowRegime.NULL }; + const flagsNull = { ENUM_StructureFlowMode: StructureFlowMode.NULL, ENUM_StructureFlowRegime: StructureFlowRegime.NULL }; if (sVarCalc === "Q") { if (this.prms.h1.v <= 0 || this.prms.Z1.v === this.prms.Z2.v || this.prms.W.v <= 0) { return new Result(0, flagsNull); @@ -170,8 +170,8 @@ export abstract class Structure extends Nub { protected getResultData() { return { - Mode: this.getFlowMode(), - Regime: this.getFlowRegime(), + ENUM_StructureFlowMode: this.getFlowMode(), + ENUM_StructureFlowRegime: this.getFlowRegime(), }; } diff --git a/src/structure/structure_cem88d.ts b/src/structure/structure_cem88d.ts index 0ddac7bdb41187a71f3eed1a07fa3bd55ae6f6a3..2040e13e26f91d4d107a6890a77d8d82767b2b16 100644 --- a/src/structure/structure_cem88d.ts +++ b/src/structure/structure_cem88d.ts @@ -31,9 +31,9 @@ export class StructureCem88d extends RectangularStructure { const b2: number = Math.sqrt(this.prms.h1.v - this.prms.h2.v); const cd1: number = cd * 2.5981; // cd * 3*sqrt(3)/2 this.debug("StructureCem88d.Equation b1=" + b1 + " b2=" + b2 + " cd1=" + cd1); - switch (data.Mode) { + switch (data.ENUM_StructureFlowMode) { case StructureFlowMode.WEIR: - switch (data.Regime) { + switch (data.ENUM_StructureFlowRegime) { case StructureFlowRegime.FREE: v = cd * this.prms.h1.v * b1; break; @@ -45,7 +45,7 @@ export class StructureCem88d extends RectangularStructure { break; case StructureFlowMode.ORIFICE: const b3: number = Math.pow(this.prms.h1.v - this.prms.W.v, 1.5); - switch (data.Regime) { + switch (data.ENUM_StructureFlowRegime) { case StructureFlowRegime.FREE: v = cd * (this.prms.h1.v * b1 - b3); break; diff --git a/src/structure/structure_cem88v.ts b/src/structure/structure_cem88v.ts index c0dc94caaaef79047b945e1fd85517529441c697..a7876538e36369e6fe7e82c87828a9224a861bc6 100644 --- a/src/structure/structure_cem88v.ts +++ b/src/structure/structure_cem88v.ts @@ -21,13 +21,13 @@ export class StructureCem88v extends RectangularStructure { let v: number; const mu0: number = 2 / 3 * this.prms.Cd.v; let KF: number; - if (data.Regime !== StructureFlowRegime.FREE) { + if (data.ENUM_StructureFlowRegime !== StructureFlowRegime.FREE) { KF = this.getKF(Math.sqrt(1 - this.prms.h2.v / this.prms.h1.v), this.getAlfa(this.prms.h2.v)); } - switch (data.Mode) { + switch (data.ENUM_StructureFlowMode) { case StructureFlowMode.WEIR: const muf: number = mu0 - 0.08; - switch (data.Regime) { + switch (data.ENUM_StructureFlowRegime) { case StructureFlowRegime.FREE: v = muf * this.prms.L.v * Structure.R2G * Math.pow(this.prms.h1.v, 1.5); break; @@ -40,12 +40,12 @@ export class StructureCem88v extends RectangularStructure { case StructureFlowMode.ORIFICE: const mu: number = mu0 - 0.08 / (this.prms.h1.v / this.prms.W.v); const mu1: number = mu0 - 0.08 / (this.prms.h1.v / this.prms.W.v - 1); - if (data.Regime === StructureFlowRegime.FREE) { + if (data.ENUM_StructureFlowRegime === StructureFlowRegime.FREE) { v = this.prms.L.v * Structure.R2G * (mu * Math.pow(this.prms.h1.v, 1.5) - mu1 * Math.pow(this.prms.h1.v - this.prms.W.v, 1.5)); } else { - if (data.Regime === StructureFlowRegime.PARTIAL) { + if (data.ENUM_StructureFlowRegime === StructureFlowRegime.PARTIAL) { v = this.prms.L.v * Structure.R2G * ( KF * mu * Math.pow(this.prms.h1.v, 1.5) - mu1 * Math.pow(this.prms.h1.v - this.prms.W.v, 1.5) diff --git a/src/structure/structure_cunge80.ts b/src/structure/structure_cunge80.ts index b278eadca75c9feb578e12a9cfdbcad4162ad789..eee8b87eacfbada092ae550eba1c34095a94979a 100644 --- a/src/structure/structure_cunge80.ts +++ b/src/structure/structure_cunge80.ts @@ -26,9 +26,9 @@ export class StructureCunge80 extends RectangularStructure { const data = this.getResultData(); let v: number; - switch (data.Regime) { + switch (data.ENUM_StructureFlowRegime) { case StructureFlowRegime.FREE: - if (data.Mode === StructureFlowMode.WEIR) { + if (data.ENUM_StructureFlowMode === StructureFlowMode.WEIR) { const R32: number = 3 * Math.sqrt(3) / 2; v = this.prms.Cd.v * this.prms.L.v * Structure.R2G / R32 * Math.pow(this.prms.h1.v, 1.5); this.debug("StructureCunge80.Equation WEIR FREE Q=" + v); @@ -39,7 +39,7 @@ export class StructureCunge80 extends RectangularStructure { } break; case StructureFlowRegime.SUBMERGED: - if (data.Mode === StructureFlowMode.WEIR) { + if (data.ENUM_StructureFlowMode === StructureFlowMode.WEIR) { v = this.prms.Cd.v * this.prms.L.v * Structure.R2G * this.prms.h2.v * Math.sqrt(this.prms.h1.v - this.prms.h2.v); this.debug("StructureCunge80.Equation WEIR SUBMERGED Q=" + v); diff --git a/src/structure/structure_kivi.ts b/src/structure/structure_kivi.ts index 082bdf614540ea028445f620f879ff7663914622..797f184f8bf7de0e3554ba51ea7164501e68681a 100644 --- a/src/structure/structure_kivi.ts +++ b/src/structure/structure_kivi.ts @@ -45,7 +45,7 @@ export class StructureKivi extends Structure { let Q = cd * this.prms.L.v * Structure.R2G * Math.pow(this.prms.h1.v, 1.5); - if (res.extraResults.Regime === StructureFlowRegime.SUBMERGED) { + if (res.extraResults.ENUM_StructureFlowRegime === StructureFlowRegime.SUBMERGED) { Q = Villemonte(this.prms.h1.v, this.prms.h2.v, 1.5) * Q; } diff --git a/src/structure/villemonte.ts b/src/structure/villemonte.ts index c2fdf090eb804fb167e3b52cc2b12e21d2565233..c78a2829f1fe457d8937177d2b63e7496005849f 100644 --- a/src/structure/villemonte.ts +++ b/src/structure/villemonte.ts @@ -2,7 +2,7 @@ * * @param h1 hauteur d'eau amont au dessus de la crête du seuil * @param h2 hauteur d'eau aval au dessus de la crête du seuil - * @param n n est lÂ’exposant dans les relations d’écoulement dénoyé : + * @param n n est l'exposant dans les relations d'écoulement dénoyé : * déversoir proportionnel : n=1 déversoir rectangulaire : n=1,5 * déversoir parabolique : n=2 déversoir triangulaire : n=2,5 */