Ticket #2376: myShallowRiver.mo

File myShallowRiver.mo, 7.2 KB (added by Adeel Asghar, 12 years ago)
Line 
1package myShallowRiver
2 // Package for simulating a Shallow River
3 // author: Bernt Lie
4 // Telemark University College
5 // September 10, 2013
6 model compareShallowRiver
7 simShallowRiver sr5(NVd = 5),sr10(NVd = 10),sr20(NVd = 20); //,sr50(NVd = 50);
8 end compareShallowRiver;
9 //
10 model simShallowRiver
11 // Main Shallow River model
12 // author: Bernt Lie
13 // Telemark University College
14 // September 10, 2013
15 //
16 // Constants
17 constant Real g = 9.81 "Gravity, m/s2";
18 // Parameters
19 parameter Real w = 100 "River width, m";
20 parameter Real L = 5e3 "River length, m";
21 parameter Real H = 57 "River bed drop, m";
22 parameter Real th = asin(H/L) "River bed angle, rad";
23 parameter Real rho = 1e3 "Density of water, kg/m3";
24 parameter Real kS = 20 "Strickler's friction factor, m0.33/s";
25 parameter Integer NVd = 18 "Number of momentum balance segments, -";
26 parameter Real dx = L/NVd "discretization step, m";
27 parameter Real oneh[NVd+1] = ones(NVd+1);
28 // Initial state parameters
29 parameter Real h0 = 4;// 0.5275 "Initial guess river level, m";
30 parameter Real Vd0 = 120 "Initial guess volumetric flow rate, m3/s";
31 // Declaring variables
32 // -- states
33 Real h[NVd+1](each start = h0); //, each fixed = true) "Level states, m";
34 Real Vd[NVd](each start = Vd0); //, each fixed = true) "Volumetric flow rates, m3/s";
35 // -- auxiliary variables at mass knots
36 Real A[NVd+1] "Area perpendicular to flow direction, m2";
37 Real per[NVd+1] "Wetting perimeter, m";
38 // -- auxiliary variables at momentum knots
39 Real h_[NVd] "Staggered levels, m";
40 Real A_[NVd] "Staggered areas, m2";
41 Real Md_rho_i[NVd] "Input momentum flow/density, m4/s2";
42 Real Md_rho_o[NVd] "Output momentum flow/density, m4/s2";
43 Real per_[NVd] "Staggered wetted perimeters, m";
44 Real R_[NVd] "Staggered hydraulic radii, m";
45 Real C_[NVd] "Staggered Chezy coefficient, ...";
46 // Real vMd[NVd] "stupid vector";
47 // Real vp[NVd] "stupid vector";
48 // Real vg[NVd] "stupid vector";
49 // Real vf[NVd] "stupid vector";
50 Real vMdip[NVd-2] "stupid vector";
51 Real vMdim[NVd-2] "stupid vector";
52 // -- input variables
53 Real Vdin(start = 120);
54 Real Vdout(start = 120);
55 // -- output variables, scaled
56 output Real _Vd, _h;
57 // Initializating in steady state
58
59 initial equation
60 der(h[1:end-1]) = zeros(NVd);
61 der(Vd[:]) = zeros(NVd);
62
63 // Equations constituting the model
64 equation
65 // Input values
66 //Vdin = 120;
67 Vdout = 120;
68
69 if time < 1e3 then
70 Vdin = 120;
71 elseif time < 1.2e3 then
72 Vdin = 145;
73 else
74 Vdin = 120;
75 end if;
76
77 // Defining auxiliary variables
78 A[:] = h[:]*w;
79 per[:] = w*oneh[:] + 2*h[:];
80 //
81 h_[:] = (h[1:end-1] + h[2:end])/2;
82 A_[:] = h_[:]*w;
83 per_[:] = (per[1:end-1] + per[2:end])/2;
84 R_[:] = A_[:]./per_[:];
85 C_[:] = kS*((R_[:]).^(1/6));
86 //
87 Md_rho_i[1] = Vdin^2/A[1];
88 for n in 2:NVd-1 loop
89 vMdip[n-1] = abs(Vd[n-1])*max(Vd[n-1],0)/A_[n-1];
90 vMdim[n-1] = abs(Vd[n+1])*max(-Vd[n+1],0)/A_[n+1];
91 Md_rho_i[n] = vMdip[n-1] + vMdim[n-1];
92 end for;
93 Md_rho_i[end] = abs(Vd[end-1])*max(Vd[end-1],0)/A_[end-1];
94 //
95 for n in 1:NVd-1 loop
96 Md_rho_o[n] = Vd[n]^2/A_[n];
97 end for;
98 Md_rho_o[NVd] = Vdout^2/A[NVd+1];
99 //
100 // vMd[:] = (Md_rho_i[:] - Md_rho_o[:])/dx;
101 // vp[:] = g*w*cos(th)*(h[1:end-1].^2 - h[2:end].^2)/(2*dx);
102 // vg[:] = A_[:]*g*sin(th);
103 // vf[:] = - g*per_[:].*abs(Vd[:]).*Vd[:]./(C_[:].*A_[:]).^2;
104 // Differential equations
105 der(h[1]) = (Vdin-Vd[1])/(w*dx/2);
106 der(h[2:NVd]) = (Vd[1:end-1] - Vd[2:end])/(w*dx);
107 der(h[NVd+1]) = (Vd[NVd] - Vdout)/(w*dx/2);
108 //
109 der(Vd[:]) = (Md_rho_i[:] - Md_rho_o[:])/dx
110 + g*w*cos(th)*(h[1:end-1].^2 - h[2:end].^2)/(2*dx)
111 + A_[:]*g*sin(th)
112 - g*per_[:].*abs(Vd[:]).*Vd[:]./(C_[:].*A_[:]).^2;
113 // Outputs
114 _Vd = Vd[end];
115 _h = h[end];
116 end simShallowRiver;
117 //
118 model simShallowRiverWorks
119 // Main Shallow River model
120 // author: Bernt Lie
121 // Telemark University College
122 // September 10, 2013
123 //
124 // Constants
125 constant Real g = 9.81 "Gravity, m/s2";
126 // Parameters
127 parameter Real w = 100 "River width, m";
128 parameter Real L = 5e3 "River length, m";
129 parameter Real H = 57 "River bed drop, m";
130 parameter Real th = asin(H/L) "River bed angle, rad";
131 parameter Real rho = 1e3 "Density of water, kg/m3";
132 parameter Real kS = 20 "Strickler's friction factor, m0.33/s";
133 parameter Integer NVd = 18 "Number of momentum balance segments, -";
134 parameter Real dx = L/NVd "discretization step, m";
135 parameter Real oneh[NVd+1] = ones(NVd+1);
136 // Initial state parameters
137 parameter Real h0 = 4;// 0.5275 "Initial guess river level, m";
138 parameter Real Vd0 = 120 "Initial guess volumetric flow rate, m3/s";
139 // Declaring variables
140 // -- states
141 Real h[NVd+1](each start = h0); //, each fixed = true) "Level states, m";
142 Real Vd[NVd](each start = Vd0); //, each fixed = true) "Volumetric flow rates, m3/s";
143 // -- auxiliary variables at mass knots
144 Real A[NVd+1] "Area perpendicular to flow direction, m2";
145 Real per[NVd+1] "Wetting perimeter, m";
146 // -- auxiliary variables at momentum knots
147 Real h_[NVd] "Staggered levels, m";
148 Real A_[NVd] "Staggered areas, m2";
149 Real Md_rho_i[NVd] "Input momentum flow/density, m4/s2";
150 Real Md_rho_o[NVd] "Output momentum flow/density, m4/s2";
151 Real per_[NVd] "Staggered wetted perimeters, m";
152 Real R_[NVd] "Staggered hydraulic radii, m";
153 Real C_[NVd] "Staggered Chezy coefficient, ...";
154 Real vMd[NVd] "stupid vector";
155 Real vp[NVd] "stupid vector";
156 Real vg[NVd] "stupid vector";
157 Real vf[NVd] "stupid vector";
158 Real vMdip[NVd-2] "stupid vector";
159 Real vMdim[NVd-2] "stupid vector";
160 // -- input variables
161 Real Vdin(start = 120);
162 Real Vdout(start = 120);
163 // -- output variables, scaled
164 output Real _Vd, _h;
165 // Initializating in steady state
166
167 initial equation
168 der(h[1:end-1]) = zeros(NVd);
169 der(Vd[:]) = zeros(NVd);
170
171 // Equations constituting the model
172 equation
173 // Input values
174 //Vdin = 120;
175 Vdout = 120;
176
177 if time < 1e3 then
178 Vdin = 120;
179 elseif time < 1.2e3 then
180 Vdin = 145;
181 else
182 Vdin = 120;
183 end if;
184
185 // Defining auxiliary variables
186 A[:] = h[:]*w;
187 per[:] = w*oneh[:] + 2*h[:];
188 //
189 h_[:] = (h[1:end-1] + h[2:end])/2;
190 A_[:] = h_[:]*w;
191 per_[:] = (per[1:end-1] + per[2:end])/2;
192 R_[:] = A_[:]./per_[:];
193 C_[:] = kS*((R_[:]).^(1/6));
194 //
195 Md_rho_i[1] = Vdin^2/A[1];
196 for n in 2:NVd-1 loop
197 vMdip[n-1] = abs(Vd[n-1])*max(Vd[n-1],0)/A_[n-1];
198 vMdim[n-1] = abs(Vd[n+1])*max(-Vd[n+1],0)/A_[n+1];
199 Md_rho_i[n] = vMdip[n-1] + vMdim[n-1];
200 end for;
201 Md_rho_i[end] = abs(Vd[end-1])*max(Vd[end-1],0)/A_[end-1];
202 //
203 for n in 1:NVd-1 loop
204 Md_rho_o[n] = Vd[n]^2/A_[n];
205 end for;
206 Md_rho_o[NVd] = Vdout^2/A[NVd+1];
207 //
208 vMd[:] = (Md_rho_i[:] - Md_rho_o[:])/dx;
209 vp[:] = g*w*cos(th)*(h[1:end-1].^2 - h[2:end].^2)/(2*dx);
210 vg[:] = A_[:]*g*sin(th);
211 vf[:] = - g*per_[:].*abs(Vd[:]).*Vd[:]./(C_[:].*A_[:]).^2;
212 // Differential equations
213 der(h[1]) = (Vdin-Vd[1])/(w*dx/2);
214 der(h[2:NVd]) = (Vd[1:end-1] - Vd[2:end])/(w*dx);
215 der(h[NVd+1]) = (Vd[NVd] - Vdout)/(w*dx/2);
216 //
217 der(Vd[:]) = vMd[:] + vp[:] + vg[:] + vf[:];
218 // Outputs
219 _Vd = Vd[end];
220 _h = h[end];
221 end simShallowRiverWorks;
222 // End package
223end myShallowRiver;
224