Calcul d' une courbe de remous avec ressaut

C. Ancey, 2020, EPFL/LHE

Paramètres

q=10;
Ch=50;
p=0.5;
g=9.81;

Valeurs dans le bief 1

i1=0.05;
hn1=(q/Ch/Sqrt[i1])^(2/3)
Frn=q/hn1^1.5/Sqrt[9.81]
hc=(q^2/9.81)^(1/3)
Fr1=q/2^1.5/Sqrt[9.81]

Valeurs dans le bief 2

i2=0.002;
hn2=(q/Ch/Sqrt[i2])^(2/3)
Fr2=q/hn2^1.5/Sqrt[9.81]

Calcul de la courbe de remous dans le bief 1

eqn1=NDSolve[{h'[x]==i1 (1-(hn1/h[x])^3)/(1-(hc/h[x])^3),h[0]==2},h,{x,0,10}]
solution1=h/.eqn1[[1]];
des0=Plot[solution1[x],{x,0,10}]
hA=solution1[10]

Calcul de la courbe de remous dans le bief 2 : prolongement de la branche supercritique

eqn2=NDSolve[{h'[x]==i2 (1-(hn2/h[x])^3)/(1-(hc/h[x])^3),h[10]==hA},h,{x,10,600}];
solution2=h/.eqn2[[1]]
xl=solution2["Domain"][[1,2]]
des1=Plot[solution2[x],{x,10,xl},PlotRange->{0,3}]

Calcul de la courbe de remous conjuguée

conj[h_]:=1/2*h*(Sqrt[8*(q/h^1.5/Sqrt[9.81])^2+1]-1)

des2=Plot[conj@solution2[x],{x,10,xl},PlotRange->All,PlotStyle->Dashing[{0.01,0.01}]];
Show[des0,des1,des2,PlotRange->All]

Hauteur d’eau en B (seuil)

Hb=(q)^(2/3)*3/2/g^(1/3)+p//N
sol=h/.Solve[h+(q/h)^2/2/g==Hb,h]

Show[Plot [h+(q/h)^2/2/g,{h,1,4}],Plot [Hb,{h,1,4},PlotStyle->Dashed],Frame->True,FrameLabel->{"h [m]","H [m] "}]
hb=sol[[3]]

Calcul de la courbe de remous dans le bief 2 : branche subcritique et courbe conjuguée

eqn3=NDSolve[{h'[x]==i2 (1-(hn2/h[x])^3)/(1-(hc/h[x])^3),h[1000]==hb},h,{x,1000,10}]
solution3=h/.eqn3[[1]];
xl2=solution3["Domain"][[1,1]];
des3=Plot[solution3[x],{x,1000,xl2},PlotRange->All,PlotStyle->Red] 
des4=Plot[conj@(solution3[x]),{x,1000,xl2},PlotRange->{0,3},PlotStyle->{Red,Dashed}]

Tracé de tous les courbes

des=Show[des0,des1,des2,des3,des4,Frame->True,Axes->False,FrameLabel->
{StyleForm["  x  ",FontSize->18,FontSlant->"Italic",FontFamily->"Times" ],StyleForm[" h ",FontSize->18,FontFamily->"Times",FontSlant->"Italic" ]},DefaultFont->{"Times",14},ImageSize->500,PlotRange->All]

Position du ressaut

xr=x/.FindRoot[conj@solution2[x]==solution3[x],{x,10,90}][[1]]
FindRoot[conj@solution3[x]==solution2[x],{x,10,90}]

Dessin final

y2=solution3[xr];
y1=solution2[xr];
Print["Le ressaut se produit en ",Round[xr 10]/10., " m et avec les hauteurs (Subscript[h, 1]) = ",Round[y1 10]/10.," et (Subscript[h, 2]) = ",Round[y2 10]/10.]


des5=Plot[solution3[x],{x,1000,xr},PlotRange->All,PlotStyle->Red];
des6=Plot[solution2[x],{x,xr,90},PlotRange->All,PlotStyle->Dashing[{0.02,0.02}]];
des7=ParametricPlot[{xr,y},{y,y1,y2},PlotStyle->{Black,Thickness[0.007]}];
des8= Plot[solution2[x],{x,10,xr},PlotRange->All];
dessin=Show[des0,des2,des8,des5,des6,des7,Frame->True,Axes->False,FrameLabel->
{StyleForm["  x  ",FontSize->18,FontSlant->"Italic",FontFamily->"Times",PrivateFontOptions->{"OperatorSubstitution"->False}],StyleForm[" h ",FontSize->18,FontFamily->"Times",FontSlant->"Italic",PrivateFontOptions->{"OperatorSubstitution"->False}]},DefaultFont->{"Times",14},ImageSize->500,Epilog->{Text["O",{2,2.14}],
Text["A",{10,1.65}],
Text["B",{1000,1.5}],
Text["D'",{xr,y2+0.2}],
Text["D",{xr-2,y1+0.1}],
Text["C",{94,2.2}]
},PlotRange->{{0,100},{1,3}}]