| 1 | """ |
|---|
| 2 | |
|---|
| 3 | Ole Check Culvert Routine from Line 258 |
|---|
| 4 | |
|---|
| 5 | Although it is Setup as a Culvert with the Opening presented vertically, |
|---|
| 6 | for now the calculation of flow rate is assuming a horizontal hole in the |
|---|
| 7 | ground (Fix this Later) |
|---|
| 8 | |
|---|
| 9 | MOST importantly 2 things... |
|---|
| 10 | 1. How to use the Create Polygon Routine to enquire Depth ( or later energy) |
|---|
| 11 | infront of the Culvert |
|---|
| 12 | |
|---|
| 13 | Done (Ole) |
|---|
| 14 | |
|---|
| 15 | 2. How to apply the Culvert velocity and thereby Momentum to the Outlet |
|---|
| 16 | Ject presented at the Outlet |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction |
|---|
| 21 | |
|---|
| 22 | This Version CALCULATES the Culvert Velocity and Uses it to establish |
|---|
| 23 | The Culvert Outlet Momentum |
|---|
| 24 | |
|---|
| 25 | The Aim is to define a Flow Transfer function that Simulates a Culvert |
|---|
| 26 | by using the Total Available Energy to Drive the Culvert |
|---|
| 27 | as Derived by determining the Difference in Total Energy |
|---|
| 28 | between 2 Points, Just Up stream and Just Down Stream of the Culvert |
|---|
| 29 | away from the influence of the Flow Abstraction etc.. |
|---|
| 30 | |
|---|
| 31 | This example includes a Model Topography that shows a |
|---|
| 32 | TYPICAL Headwall Configuration |
|---|
| 33 | |
|---|
| 34 | The aim is to change the Culvert Routine to Model more precisely the |
|---|
| 35 | abstraction |
|---|
| 36 | from a vertical face. |
|---|
| 37 | |
|---|
| 38 | The inflow must include the impact of Approach velocity. |
|---|
| 39 | Similarly the Outflow has MOMENTUM Not just Up welling as in the |
|---|
| 40 | Horizontal Style |
|---|
| 41 | abstraction |
|---|
| 42 | |
|---|
| 43 | """ |
|---|
| 44 | |
|---|
| 45 | #------------------------------------------------------------------------------ |
|---|
| 46 | # Import necessary modules |
|---|
| 47 | #------------------------------------------------------------------------------ |
|---|
| 48 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
|---|
| 49 | from anuga.shallow_water import Domain |
|---|
| 50 | from anuga.shallow_water.shallow_water_domain import Reflective_boundary |
|---|
| 51 | from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary |
|---|
| 52 | from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing |
|---|
| 53 | from anuga.culvert_flows.culvert_polygons import create_culvert_polygons |
|---|
| 54 | from anuga.utilities.polygon import plot_polygons |
|---|
| 55 | |
|---|
| 56 | from math import pi,sqrt |
|---|
| 57 | |
|---|
| 58 | #------------------------------------------------------------------------------ |
|---|
| 59 | # Setup computational domain |
|---|
| 60 | #------------------------------------------------------------------------------ |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | # Open file for storing some specific results... |
|---|
| 64 | fid = open('Culvert_Headwall_VarM.txt', 'w') |
|---|
| 65 | |
|---|
| 66 | length = 40. |
|---|
| 67 | width = 5. |
|---|
| 68 | |
|---|
| 69 | #dx = dy = 1 # Resolution: Length of subdivisions on both axes |
|---|
| 70 | #dx = dy = .5 # Resolution: Length of subdivisions on both axes |
|---|
| 71 | dx = dy = .25 # Resolution: Length of subdivisions on both axes |
|---|
| 72 | #dx = dy = .1 # Resolution: Length of subdivisions on both axes |
|---|
| 73 | |
|---|
| 74 | # OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???...... |
|---|
| 75 | # Can I refine in a X & Y Range ??? |
|---|
| 76 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
|---|
| 77 | len1=length, len2=width) |
|---|
| 78 | domain = Domain(points, vertices, boundary) |
|---|
| 79 | domain.set_name('culvert_HW_Var_Mom') # Output name |
|---|
| 80 | domain.set_default_order(2) |
|---|
| 81 | domain.H0 = 0.01 |
|---|
| 82 | domain.tight_slope_limiters = True |
|---|
| 83 | domain.set_minimum_storable_height(0.02) |
|---|
| 84 | print 'Size', len(domain) |
|---|
| 85 | |
|---|
| 86 | |
|---|
| 87 | #------------------------------------------------------------------------------ |
|---|
| 88 | # Setup initial conditions |
|---|
| 89 | #------------------------------------------------------------------------------ |
|---|
| 90 | |
|---|
| 91 | # Define the topography (land scape) |
|---|
| 92 | def topography(x, y): |
|---|
| 93 | """Set up a weir |
|---|
| 94 | |
|---|
| 95 | A culvert will connect either side of an Embankment with a Headwall type structure |
|---|
| 96 | The aim is for the Model to REALISTICALY model flow through the Culvert |
|---|
| 97 | """ |
|---|
| 98 | # General Slope of Topography |
|---|
| 99 | z=-x/100 |
|---|
| 100 | |
|---|
| 101 | # Add bits and Pieces to topography |
|---|
| 102 | N = len(x) |
|---|
| 103 | for i in range(N): |
|---|
| 104 | |
|---|
| 105 | # Sloping Embankment Across Channel |
|---|
| 106 | if 5.0 < x[i] < 10.1: |
|---|
| 107 | if 1.0+(x[i]-5.0)/5.0 < y[i] < 4.0 - (x[i]-5.0)/5.0: # Cut Out Segment for Culvert FACE |
|---|
| 108 | z[i]=z[i] |
|---|
| 109 | else: |
|---|
| 110 | z[i] += 0.5*(x[i] -5.0) # Sloping Segment U/S Face |
|---|
| 111 | if 10.0 < x[i] < 12.1: |
|---|
| 112 | z[i] += 2.5 # Flat Crest of Embankment |
|---|
| 113 | if 12.0 < x[i] < 14.5: |
|---|
| 114 | if 2.0-(x[i]-12.0)/2.5 < y[i] < 3.0 + (x[i]-12.0)/2.5: # Cut Out Segment for Culvert FACE |
|---|
| 115 | z[i]=z[i] |
|---|
| 116 | else: |
|---|
| 117 | z[i] += 2.5-1.0*(x[i] -12.0) # Sloping D/S Face |
|---|
| 118 | |
|---|
| 119 | |
|---|
| 120 | # Constriction |
|---|
| 121 | #if 27 < x[i] < 29 and y[i] > 3: |
|---|
| 122 | # z[i] += 2 |
|---|
| 123 | |
|---|
| 124 | # Pole |
|---|
| 125 | #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2: |
|---|
| 126 | # z[i] += 2 |
|---|
| 127 | |
|---|
| 128 | # HOLE For Pit at Opening[0] |
|---|
| 129 | #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2: |
|---|
| 130 | # z[i]-=1 |
|---|
| 131 | |
|---|
| 132 | # HOLE For Pit at Opening[1] |
|---|
| 133 | #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2: |
|---|
| 134 | # z[i]-=1 |
|---|
| 135 | |
|---|
| 136 | return z |
|---|
| 137 | |
|---|
| 138 | |
|---|
| 139 | domain.set_quantity('elevation', topography) # Use function for elevation |
|---|
| 140 | domain.set_quantity('friction', 0.01) # Constant friction |
|---|
| 141 | domain.set_quantity('stage', |
|---|
| 142 | expression='elevation') # Dry initial condition |
|---|
| 143 | |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | |
|---|
| 147 | #------------------------------------------------------------------------------ |
|---|
| 148 | # Setup specialised forcing terms |
|---|
| 149 | #------------------------------------------------------------------------------ |
|---|
| 150 | |
|---|
| 151 | # NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet) |
|---|
| 152 | # |
|---|
| 153 | # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed |
|---|
| 154 | # Flow Is Removed at a Rate of INFLOW |
|---|
| 155 | # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges |
|---|
| 156 | # |
|---|
| 157 | # This SHould be changed to a Verical Opening Both BOX and Circular |
|---|
| 158 | # There will be several Culvert Routines such as: |
|---|
| 159 | # CULVERT_Simple_FLOOR |
|---|
| 160 | # CULVERT_Simple_WALL |
|---|
| 161 | # CULVERT_Eqn_Floor |
|---|
| 162 | # CULVERT_Eqn_Wall |
|---|
| 163 | # CULVERT_Tab_Floor |
|---|
| 164 | # CULVERT_Tab_Wall |
|---|
| 165 | # BRIDGES..... |
|---|
| 166 | # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!! |
|---|
| 167 | |
|---|
| 168 | # COULD USE EPA SWMM Model !!!! |
|---|
| 169 | |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | class Culvert_flow: |
|---|
| 173 | """Culvert flow - transfer water from one hole to another |
|---|
| 174 | |
|---|
| 175 | Using Momentum as Calculated by Culvert Flow !! |
|---|
| 176 | Could be Several Methods Investigated to do This !!! |
|---|
| 177 | |
|---|
| 178 | 2008_May_08 |
|---|
| 179 | To Ole: |
|---|
| 180 | OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on |
|---|
| 181 | the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon |
|---|
| 182 | |
|---|
| 183 | The two centers are now passed on to create_polygon. |
|---|
| 184 | |
|---|
| 185 | |
|---|
| 186 | Input: Two points, pipe_size (either diameter or width, height), mannings_rougness, |
|---|
| 187 | inlet/outlet energy_loss_coefficients, internal_bend_coefficent, |
|---|
| 188 | top-down_blockage_factor and bottom_up_blockage_factor |
|---|
| 189 | |
|---|
| 190 | |
|---|
| 191 | And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront |
|---|
| 192 | of the culvert |
|---|
| 193 | At the moment this script uses only Depth, later we can change it to Energy... |
|---|
| 194 | |
|---|
| 195 | Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity |
|---|
| 196 | The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet... |
|---|
| 197 | |
|---|
| 198 | """ |
|---|
| 199 | |
|---|
| 200 | def __init__(self, |
|---|
| 201 | domain, |
|---|
| 202 | end_point0=None, |
|---|
| 203 | end_point1=None, |
|---|
| 204 | width=None, |
|---|
| 205 | height=None, |
|---|
| 206 | verbose=False): |
|---|
| 207 | |
|---|
| 208 | from Numeric import sqrt, sum |
|---|
| 209 | |
|---|
| 210 | if height is None: height = width |
|---|
| 211 | |
|---|
| 212 | # Create the fundamental culvert polygons |
|---|
| 213 | P = create_culvert_polygons(end_point0, |
|---|
| 214 | end_point1, |
|---|
| 215 | width=width, |
|---|
| 216 | height=height) |
|---|
| 217 | |
|---|
| 218 | if verbose is True: |
|---|
| 219 | pass |
|---|
| 220 | #plot_polygons([[end_point0, end_point1], |
|---|
| 221 | # P['exchange_polygon0'], |
|---|
| 222 | # P['exchange_polygon1'], |
|---|
| 223 | # P['enquiry_polygon0'], |
|---|
| 224 | # P['enquiry_polygon1']], |
|---|
| 225 | # figname='culvert_polygon_output') |
|---|
| 226 | |
|---|
| 227 | self.openings = [] |
|---|
| 228 | self.openings.append(Inflow(domain, |
|---|
| 229 | polygon=P['exchange_polygon0'])) |
|---|
| 230 | |
|---|
| 231 | self.openings.append(Inflow(domain, |
|---|
| 232 | polygon=P['exchange_polygon1'])) |
|---|
| 233 | |
|---|
| 234 | |
|---|
| 235 | # Assume two openings for now: Referred to as 0 and 1 |
|---|
| 236 | assert len(self.openings) == 2 |
|---|
| 237 | |
|---|
| 238 | # Store basic geometry |
|---|
| 239 | self.end_points = [end_point0, end_point1] |
|---|
| 240 | self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']] |
|---|
| 241 | self.vector = P['vector'] |
|---|
| 242 | self.distance = P['length'] |
|---|
| 243 | self.verbose = verbose |
|---|
| 244 | self.width = width |
|---|
| 245 | self.height = height |
|---|
| 246 | self.last_time = 0.0 |
|---|
| 247 | |
|---|
| 248 | # Create objects to update momentum (a bit crude at this stage) |
|---|
| 249 | self.xmom_forcing0 = General_forcing(domain, |
|---|
| 250 | 'xmomentum', |
|---|
| 251 | polygon=P['exchange_polygon0']) |
|---|
| 252 | |
|---|
| 253 | self.xmom_forcing1 = General_forcing(domain, |
|---|
| 254 | 'xmomentum', |
|---|
| 255 | polygon=P['exchange_polygon1']) |
|---|
| 256 | |
|---|
| 257 | self.ymom_forcing0 = General_forcing(domain, |
|---|
| 258 | 'ymomentum', |
|---|
| 259 | polygon=P['exchange_polygon0']) |
|---|
| 260 | |
|---|
| 261 | self.ymom_forcing1 = General_forcing(domain, |
|---|
| 262 | 'ymomentum', |
|---|
| 263 | polygon=P['exchange_polygon1']) |
|---|
| 264 | |
|---|
| 265 | # Print something |
|---|
| 266 | |
|---|
| 267 | print 'Culvert Effective Length =', self.distance |
|---|
| 268 | print 'Culvert Slope is Delta Z / Dist' |
|---|
| 269 | print 'Culvert Direction is ', self.vector |
|---|
| 270 | print 'Point 1m Up Stream is X,Y =' |
|---|
| 271 | print 'Point 1m Down Stream is X,Y =' |
|---|
| 272 | |
|---|
| 273 | |
|---|
| 274 | |
|---|
| 275 | def __call__(self, domain): |
|---|
| 276 | from anuga.utilities.numerical_tools import mean |
|---|
| 277 | from anuga.utilities.polygon import inside_polygon |
|---|
| 278 | from anuga.config import g, epsilon |
|---|
| 279 | from Numeric import take, sqrt |
|---|
| 280 | |
|---|
| 281 | |
|---|
| 282 | # Get average water depths at each opening |
|---|
| 283 | |
|---|
| 284 | |
|---|
| 285 | time = domain.get_time() |
|---|
| 286 | |
|---|
| 287 | print 'Time = %f, delta_t = %f' %(time, time-self.last_time) |
|---|
| 288 | |
|---|
| 289 | openings = self.openings |
|---|
| 290 | for i, opening in enumerate(openings): |
|---|
| 291 | # Quantities corresponding to fluid exchange field for this opening |
|---|
| 292 | stage_o = domain.quantities['stage'].get_values(location='centroids', |
|---|
| 293 | indices=opening.indices) |
|---|
| 294 | elevation_o = domain.quantities['elevation'].get_values(location='centroids', |
|---|
| 295 | indices=opening.indices) |
|---|
| 296 | xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', |
|---|
| 297 | indices=opening.indices) |
|---|
| 298 | ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', |
|---|
| 299 | indices=opening.indices) |
|---|
| 300 | |
|---|
| 301 | # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account) |
|---|
| 302 | depth_o = stage_o - elevation_o |
|---|
| 303 | |
|---|
| 304 | ux_o = xmomentum_o/(depth_o+epsilon) |
|---|
| 305 | uy_o = ymomentum_o/(depth_o+epsilon) |
|---|
| 306 | v_o = mean(sqrt(ux_o**2+uy_o**2)) |
|---|
| 307 | d_o = mean(depth_o) |
|---|
| 308 | w_o = mean(stage_o) |
|---|
| 309 | z_o = mean(elevation_o) |
|---|
| 310 | self.mean_xmomentum_o = mean(xmomentum_o) |
|---|
| 311 | self.mean_ymomentum_o = mean(ymomentum_o) |
|---|
| 312 | |
|---|
| 313 | # Indices corresponding to energy enquiry field for this opening |
|---|
| 314 | coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y) |
|---|
| 315 | idx = inside_polygon(coordinates, self.enquiry_polygons[i]) |
|---|
| 316 | |
|---|
| 317 | if self.verbose: |
|---|
| 318 | pass |
|---|
| 319 | #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\ |
|---|
| 320 | # %(i, len(idx), self.enquiry_polygons[i]) |
|---|
| 321 | |
|---|
| 322 | |
|---|
| 323 | # Get average model values for points in |
|---|
| 324 | # enquiry polygon for this opening |
|---|
| 325 | dq = domain.quantities |
|---|
| 326 | stage = dq['stage'].get_values(location='centroids', indices=idx) |
|---|
| 327 | xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx) |
|---|
| 328 | ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx) |
|---|
| 329 | elevation = dq['elevation'].get_values(location='centroids', indices=idx) |
|---|
| 330 | depth = stage - elevation |
|---|
| 331 | |
|---|
| 332 | # Compute mean velocity in the area in front of the culvert (taking zero depths into account) |
|---|
| 333 | ux = xmomentum/(depth+epsilon) |
|---|
| 334 | uy = ymomentum/(depth+epsilon) |
|---|
| 335 | v = mean(sqrt(ux**2+uy**2)) |
|---|
| 336 | d = mean(depth) |
|---|
| 337 | w = mean(stage) |
|---|
| 338 | z = mean(elevation) |
|---|
| 339 | self.mean_xmomentum = mean(xmomentum) |
|---|
| 340 | self.mean_ymomentum = mean(ymomentum) |
|---|
| 341 | |
|---|
| 342 | # Ratio of depth to culvert height |
|---|
| 343 | ratio = d/(2*self.height) |
|---|
| 344 | if ratio > 1.0: # Assume culvert is running full & |
|---|
| 345 | ratio = 1.0 # under pressure. Note this is usually ~ 1.35 |
|---|
| 346 | opening.ratio = ratio |
|---|
| 347 | |
|---|
| 348 | |
|---|
| 349 | # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening |
|---|
| 350 | E = d + 0.5*v**2/g # RUDY - please check this |
|---|
| 351 | opening.energy = E |
|---|
| 352 | |
|---|
| 353 | |
|---|
| 354 | |
|---|
| 355 | # Store current average stage and depth at enquiry field with each opening object |
|---|
| 356 | opening.depth = d |
|---|
| 357 | opening.stage = w |
|---|
| 358 | opening.elevation = z |
|---|
| 359 | |
|---|
| 360 | |
|---|
| 361 | # Now we are done calculating energy etc for each opening. |
|---|
| 362 | # At this point we can work on the transfer functions etc. |
|---|
| 363 | |
|---|
| 364 | # Handy values (all calculated at enquiry polygon - |
|---|
| 365 | # if you need them at exchange polygons we can easily do that.) |
|---|
| 366 | delta_h = openings[1].stage - openings[0].stage |
|---|
| 367 | #delta_h = openings[1].depth - openings[0].depth |
|---|
| 368 | delta_z = openings[1].elevation - openings[0].elevation |
|---|
| 369 | |
|---|
| 370 | # Ideas..... |
|---|
| 371 | if openings[0].depth > 0 and openings[1].depth > 0: |
|---|
| 372 | flow_rate = delta_h |
|---|
| 373 | else: |
|---|
| 374 | flow_rate = 0.0 |
|---|
| 375 | |
|---|
| 376 | if openings[0].depth > self.height: |
|---|
| 377 | # This could be usefull mayby |
|---|
| 378 | #print 'Inlet has been overflowed' |
|---|
| 379 | pass |
|---|
| 380 | |
|---|
| 381 | if openings[1].depth > self.height: |
|---|
| 382 | # This could be usefull mayby |
|---|
| 383 | #print 'Outlet has been overflowed' |
|---|
| 384 | pass |
|---|
| 385 | |
|---|
| 386 | |
|---|
| 387 | if delta_h > 0: |
|---|
| 388 | # Water Level U/S is higher than DS |
|---|
| 389 | flow_rate_orifice = 0.6*openings[0].area*(2*g*delta_h)**0.5 # Orifice Eqn Q= cA(2gh)^0.5 |
|---|
| 390 | flow_rate_weir = 1.69*(openings[0].area)*openings[0].depth**1.5 # WEIR Eqn Q= CLH^1.5 |
|---|
| 391 | |
|---|
| 392 | # Does Weir or Orifice Control Flow Rate ? |
|---|
| 393 | if flow_rate_weir > flow_rate_orifice: |
|---|
| 394 | flow_rate_control = flow_rate_orifice |
|---|
| 395 | else: |
|---|
| 396 | flow_rate_control = flow_rate_weir |
|---|
| 397 | |
|---|
| 398 | elif delta_h < 0: |
|---|
| 399 | # Water Level D/S is higher than US |
|---|
| 400 | # That is reverse flow in culvert |
|---|
| 401 | flow_rate_orifice = 0.6*openings[0].area*(2*g*-delta_h)**0.5 # Orifice Eqn Q= cA(2gh)^0.5 |
|---|
| 402 | flow_rate_weir = 1.69*openings[0].area*openings[0].depth**1.5 # WEIR Eqn Q= CLH^1.5 |
|---|
| 403 | |
|---|
| 404 | # Does Weir or Orifice Control Flow Rate ? |
|---|
| 405 | if flow_rate_weir > flow_rate_orifice: |
|---|
| 406 | flow_rate_control = flow_rate_orifice |
|---|
| 407 | else: |
|---|
| 408 | flow_rate_control = flow_rate_weir |
|---|
| 409 | |
|---|
| 410 | if openings[0].depth < epsilon: |
|---|
| 411 | # Do nothing because No water over Opening That is Set Flow to Zero!! |
|---|
| 412 | self.openings[0].rate = 0 |
|---|
| 413 | self.openings[1].rate = 0 |
|---|
| 414 | |
|---|
| 415 | elif openings[0].depth > 0: |
|---|
| 416 | # Only calculate fllow if there is some depth over the inlet |
|---|
| 417 | #if delta_h > 0: |
|---|
| 418 | if 1: |
|---|
| 419 | # FIXME (OLE): We never get in here. Why? |
|---|
| 420 | |
|---|
| 421 | #print 'We got water at inlet and dh > 0' |
|---|
| 422 | |
|---|
| 423 | # Flow will go from opening 0 to opening 1 |
|---|
| 424 | #- that is abstract from [0] and add to [1] |
|---|
| 425 | self.openings[0].rate = -flow_rate_control |
|---|
| 426 | self.openings[1].rate = flow_rate_control |
|---|
| 427 | |
|---|
| 428 | # Add jet in the form of absolute momentum to opening 1 |
|---|
| 429 | #speed = 20.0 # m/s |
|---|
| 430 | # This Should be Based on the VELOCITY in the Culvert |
|---|
| 431 | # Based on the Flow Depth in the Culvert for part full |
|---|
| 432 | # or Flow Jet of the Pressurised Culvert if FUll |
|---|
| 433 | # If Part Full Flow Calculate Part Full Depth |
|---|
| 434 | # Based on Depth Calculate Area... the Vel = flow_rate_control / Area |
|---|
| 435 | # Need to Break Velocity into X & Y Components |
|---|
| 436 | #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum') |
|---|
| 437 | # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter |
|---|
| 438 | |
|---|
| 439 | #if (ratio*self.area)== 0: # Don't Really need this already established water depth here |
|---|
| 440 | # outlet_vel=0.0 |
|---|
| 441 | #else: |
|---|
| 442 | |
|---|
| 443 | |
|---|
| 444 | # FIXME (Ole): Shouldn't this be openings[1]?? |
|---|
| 445 | outlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area)) |
|---|
| 446 | outlet_dep = 2.0*openings[0].depth*openings[0].ratio #????? |
|---|
| 447 | outlet_mom = outlet_vel*outlet_dep |
|---|
| 448 | |
|---|
| 449 | |
|---|
| 450 | # Eventually will Need Momentum in X & Y Components based on the orientation of |
|---|
| 451 | # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine |
|---|
| 452 | # YES - use self.vector which is a unit vector in the direction of the culvert. |
|---|
| 453 | |
|---|
| 454 | outlet_mom_x, outlet_mom_y = self.vector * outlet_mom |
|---|
| 455 | #print 'Directional momentum', outlet_mom_x, outlet_mom_y |
|---|
| 456 | |
|---|
| 457 | # Update the momentum forcing terms with directional momentum at the outlet |
|---|
| 458 | delta_t = time - self.last_time |
|---|
| 459 | if delta_t > 0.0: |
|---|
| 460 | xmomentum_rate = outlet_mom_x - self.mean_xmomentum |
|---|
| 461 | if xmomentum_rate > 0: |
|---|
| 462 | xmomentum_rate /= delta_t |
|---|
| 463 | else: |
|---|
| 464 | xmomentum_rate = 0.0 |
|---|
| 465 | |
|---|
| 466 | ymomentum_rate = outlet_mom_y - self.mean_ymomentum |
|---|
| 467 | if ymomentum_rate > 0: |
|---|
| 468 | ymomentum_rate /= delta_t |
|---|
| 469 | else: |
|---|
| 470 | ymomentum_rate = 0.0 |
|---|
| 471 | else: |
|---|
| 472 | xmomentum_rate = ymomentum_rate = 0.0 |
|---|
| 473 | |
|---|
| 474 | # Set momentum rates for outlet jet |
|---|
| 475 | self.xmom_forcing1.rate = xmomentum_rate |
|---|
| 476 | self.ymom_forcing1.rate = ymomentum_rate |
|---|
| 477 | |
|---|
| 478 | # Remember this value for next step (IMPORTANT) |
|---|
| 479 | self.xmom_forcing1.value = outlet_mom_x |
|---|
| 480 | self.ymom_forcing1.value = outlet_mom_y |
|---|
| 481 | |
|---|
| 482 | |
|---|
| 483 | if int(time*100) % 100 == 0: |
|---|
| 484 | s = 'T=%.3f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\ |
|---|
| 485 | %(time, flow_rate_control, outlet_vel) |
|---|
| 486 | s += ' Depth= %0.3f\n Outlet Momentum = (%0.3f, %0.3f)\n'\ |
|---|
| 487 | %(outlet_dep, outlet_mom_x, outlet_mom_y) |
|---|
| 488 | s += ' Avg Momentum at opening = (%0.3f, %0.3f)\n'\ |
|---|
| 489 | %(self.mean_xmomentum_o, self.mean_ymomentum_o) |
|---|
| 490 | s += ' Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\ |
|---|
| 491 | %(self.mean_xmomentum, self.mean_ymomentum) |
|---|
| 492 | s += ' Momentum rate: (%.4f, %.4f)'\ |
|---|
| 493 | %(xmomentum_rate, ymomentum_rate) |
|---|
| 494 | fid.write(s) |
|---|
| 495 | print s |
|---|
| 496 | else: |
|---|
| 497 | # Opening 1 has the greatest depth. Therefore Reverse the Flow !!! |
|---|
| 498 | # Flow will go from opening 1 to opening 0, That is Abstract from [1] and add to [0] |
|---|
| 499 | self.openings[0].rate = flow_rate_control # Else it will be Orifice Flow (Going US) |
|---|
| 500 | self.openings[1].rate = -flow_rate_control |
|---|
| 501 | |
|---|
| 502 | # Second Else.... if water at outlet before at inlet !!! |
|---|
| 503 | |
|---|
| 504 | |
|---|
| 505 | self.openings[1].rate = 10 |
|---|
| 506 | # Execute flow term for each opening |
|---|
| 507 | # This is where Inflow objects are evaluated and update the domain |
|---|
| 508 | for opening in self.openings: |
|---|
| 509 | opening(domain) |
|---|
| 510 | |
|---|
| 511 | # Execute momentum terms |
|---|
| 512 | # This is where Inflow objects are evaluated and update the domain |
|---|
| 513 | self.xmom_forcing0(domain) |
|---|
| 514 | self.ymom_forcing0(domain) |
|---|
| 515 | self.xmom_forcing1(domain) |
|---|
| 516 | self.ymom_forcing1(domain) |
|---|
| 517 | |
|---|
| 518 | |
|---|
| 519 | |
|---|
| 520 | # Print out flows every 1 seconds |
|---|
| 521 | if int(time*100) % 100 == 0: |
|---|
| 522 | s = 'Time %.3f\n' %time |
|---|
| 523 | s += ' Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ |
|---|
| 524 | %(openings[0].depth, |
|---|
| 525 | openings[0].area, |
|---|
| 526 | openings[0].energy, |
|---|
| 527 | openings[0].ratio) |
|---|
| 528 | s += ' Opening[1]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ |
|---|
| 529 | %(openings[1].depth, |
|---|
| 530 | openings[1].area, |
|---|
| 531 | openings[1].energy, |
|---|
| 532 | openings[1].ratio) |
|---|
| 533 | s += ' Distance=%.2f, W=%.3f, O=%.3f, C=%.3f\n'\ |
|---|
| 534 | %(self.distance, |
|---|
| 535 | flow_rate_weir, |
|---|
| 536 | flow_rate_orifice, |
|---|
| 537 | flow_rate_control) |
|---|
| 538 | |
|---|
| 539 | print s |
|---|
| 540 | |
|---|
| 541 | fid.write(s) |
|---|
| 542 | |
|---|
| 543 | # Store value of time |
|---|
| 544 | self.last_time = time |
|---|
| 545 | |
|---|
| 546 | #------------------------------------------------------------------------------ |
|---|
| 547 | # Setup culvert inlets and outlets in current topography |
|---|
| 548 | #------------------------------------------------------------------------------ |
|---|
| 549 | |
|---|
| 550 | # Define culvert inlet and outlets |
|---|
| 551 | culvert = Culvert_flow(domain, |
|---|
| 552 | end_point0=[9.0, 2.5], |
|---|
| 553 | end_point1=[13.0, 2.5], |
|---|
| 554 | width=1.00, |
|---|
| 555 | verbose=True) |
|---|
| 556 | |
|---|
| 557 | domain.forcing_terms.append(culvert) |
|---|
| 558 | |
|---|
| 559 | |
|---|
| 560 | #------------------------------------------------------------------------------ |
|---|
| 561 | # Setup boundary conditions |
|---|
| 562 | #------------------------------------------------------------------------------ |
|---|
| 563 | #Bi = Dirichlet_boundary([0.5, 0.0, 0.0]) # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!! |
|---|
| 564 | Bi = Dirichlet_boundary([0.0, 0.0, 0.0]) # Inflow based on Flow Depth and Approaching Momentum !!! |
|---|
| 565 | Br = Reflective_boundary(domain) # Solid reflective wall |
|---|
| 566 | Bo = Dirichlet_boundary([-5, 0, 0]) # Outflow |
|---|
| 567 | |
|---|
| 568 | domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br}) |
|---|
| 569 | |
|---|
| 570 | #------------------------------------------------------------------------------ |
|---|
| 571 | # Setup Application of specialised forcing terms |
|---|
| 572 | #------------------------------------------------------------------------------ |
|---|
| 573 | |
|---|
| 574 | # This is the new element implemented by Ole to allow direct input of Inflow in m^3/s |
|---|
| 575 | fixed_flow = Inflow(domain, |
|---|
| 576 | rate=6.00, |
|---|
| 577 | center=(2.1, 2.1), |
|---|
| 578 | radius=1.261566) # Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s |
|---|
| 579 | |
|---|
| 580 | # flow=file_function('Q/QPMF_Rot_Sub13.tms')) # Read Time Series in from File |
|---|
| 581 | # flow=lambda t: min(0.01*t, 0.01942)) # Time Varying Function Tap turning up |
|---|
| 582 | |
|---|
| 583 | domain.forcing_terms.append(fixed_flow) |
|---|
| 584 | |
|---|
| 585 | |
|---|
| 586 | #------------------------------------------------------------------------------ |
|---|
| 587 | # Evolve system through time |
|---|
| 588 | #------------------------------------------------------------------------------ |
|---|
| 589 | |
|---|
| 590 | |
|---|
| 591 | |
|---|
| 592 | for t in domain.evolve(yieldstep = 0.1, finaltime = 10): |
|---|
| 593 | pass |
|---|
| 594 | #if int(domain.time*100) % 100 == 0: |
|---|
| 595 | # domain.write_time() |
|---|
| 596 | |
|---|
| 597 | #if domain.get_time() >= 4 and tap.rate != 0.0: |
|---|
| 598 | # print 'Turning tap off' |
|---|
| 599 | # tap.rate = 0.0 |
|---|
| 600 | |
|---|
| 601 | #if domain.get_time() >= 3 and sink.rate < 0.0: |
|---|
| 602 | # print 'Turning drain on' |
|---|
| 603 | # sink.rate = -0.8 |
|---|
| 604 | # Close |
|---|
| 605 | |
|---|
| 606 | fid.close() |
|---|
| 607 | |
|---|
| 608 | |
|---|
| 609 | #------------------------------------------------------------------------------ |
|---|
| 610 | # Query output |
|---|
| 611 | #------------------------------------------------------------------------------ |
|---|
| 612 | |
|---|
| 613 | from anuga.shallow_water.data_manager import get_flow_through_cross_section |
|---|
| 614 | |
|---|
| 615 | swwfilename = domain.get_name()+'.sww' # Output name from script |
|---|
| 616 | print swwfilename |
|---|
| 617 | |
|---|
| 618 | polyline = [[17., 0.], [17., 5.]] |
|---|
| 619 | |
|---|
| 620 | time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True) |
|---|
| 621 | |
|---|
| 622 | from pylab import ion, plot |
|---|
| 623 | ion() |
|---|
| 624 | plot(time, Q) |
|---|
| 625 | raw_input('done') |
|---|
