source: anuga_work/development/Rudy_vandrie_2007/culvert_development.py @ 5314

Last change on this file since 5314 was 5314, checked in by ole, 16 years ago

Work on culverts

File size: 22.8 KB
Line 
1"""
2
3Ole Check Culvert Routine from Line 258
4
5Although it is Setup as a Culvert with the Opening presented vertically,
6for now the calculation of flow rate is assuming a horizontal hole in the
7ground (Fix this Later)
8
9MOST importantly 2 things...
101. How to use the Create Polygon Routine to enquire Depth ( or later energy)
11infront of the Culvert
12
13Done (Ole)
14
152. How to apply the Culvert velocity and thereby Momentum to the Outlet
16Ject presented at the Outlet
17
18
19
20Testing CULVERT (Changing from Horizontal Abstraction to Vertical Abstraction
21
22This Version CALCULATES the Culvert Velocity and Uses it to establish
23The Culvert Outlet Momentum
24
25The Aim is to define a Flow Transfer function that Simulates a Culvert
26by using the Total Available Energy to Drive the Culvert
27as Derived by determining the Difference in Total Energy
28between 2 Points, Just Up stream and Just Down Stream of the Culvert
29away from the influence of the Flow Abstraction etc..
30
31This example includes a Model Topography that shows a
32TYPICAL Headwall Configuration
33
34The aim is to change the Culvert Routine to Model more precisely the
35abstraction
36from a vertical face.
37
38The inflow must include the impact of Approach velocity.
39Similarly the Outflow has MOMENTUM Not just Up welling as in the
40Horizontal Style
41abstraction
42
43"""
44
45#------------------------------------------------------------------------------
46# Import necessary modules
47#------------------------------------------------------------------------------
48from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
49from anuga.shallow_water import Domain
50from anuga.shallow_water.shallow_water_domain import Reflective_boundary
51from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
52from anuga.shallow_water.shallow_water_domain import Inflow, General_forcing
53from anuga.culvert_flows.culvert_polygons import create_culvert_polygons
54from anuga.utilities.polygon import plot_polygons
55
56from math import pi,sqrt
57
58#------------------------------------------------------------------------------
59# Setup computational domain
60#------------------------------------------------------------------------------
61
62
63# Open file for storing some specific results...
64fid = open('Culvert_Headwall_VarM.txt', 'w')
65
66length = 40.
67width = 5.
68
69#dx = dy = 1           # Resolution: Length of subdivisions on both axes
70#dx = dy = .5           # Resolution: Length of subdivisions on both axes
71dx = 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 ???
76points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
77                                               len1=length, len2=width)
78domain = Domain(points, vertices, boundary)   
79domain.set_name('culvert_HW_Var_Mom')                 # Output name
80domain.set_default_order(2)
81domain.H0 = 0.01
82domain.tight_slope_limiters = True
83domain.set_minimum_storable_height(0.02)
84print 'Size', len(domain)
85
86
87#------------------------------------------------------------------------------
88# Setup initial conditions
89#------------------------------------------------------------------------------
90
91# Define the topography (land scape)
92def 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
139domain.set_quantity('elevation', topography)  # Use function for elevation
140domain.set_quantity('friction', 0.01)         # Constant friction
141domain.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 
172class 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        openings = self.openings
287        for i, opening in enumerate(openings):
288            stage = domain.quantities['stage'].get_values(location='centroids',
289                                                          indices=opening.indices)
290            elevation = domain.quantities['elevation'].get_values(location='centroids',
291                                                                  indices=opening.indices)
292
293            # Indices corresponding to energy enquiry field for this opening
294            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
295            idx = inside_polygon(coordinates, self.enquiry_polygons[i]) 
296           
297            if self.verbose:
298                pass
299                #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
300                #      %(i, len(idx), self.enquiry_polygons[i])
301               
302           
303            # Get average model values for points in
304            # enquiry polygon for this opening
305            dq = domain.quantities
306            stage = dq['stage'].get_values(location='centroids', indices=idx)
307            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx)
308            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx)                       
309            elevation = dq['elevation'].get_values(location='centroids', indices=idx)
310            depth = stage - elevation           
311
312            # Compute mean velocity in the area in front of the culvert (taking zero depths into account)
313            ux = xmomentum/(depth+epsilon)
314            uy = ymomentum/(depth+epsilon)
315            v = mean(sqrt(ux**2+uy**2))
316            d = mean(depth)
317            w = mean(stage)
318            z = mean(elevation)
319           
320            # Ratio of depth to culvert height
321            ratio = d/(2*self.height)
322            if ratio > 1.0:    # Assume culvert is running full &
323                ratio = 1.0    # under pressure. Note this is usually ~ 1.35
324            opening.ratio = ratio
325               
326
327            # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening
328            E = d + 0.5*v**2/# RUDY - please check this
329            opening.energy = E
330           
331           
332
333            # Store current average stage and depth at enquiry field with each opening object
334            opening.depth = d
335            opening.stage = w
336            opening.elevation = z
337           
338
339        # Now we are done calculating energy etc for each opening.
340        # At this point we can work on the transfer functions etc.
341       
342        # Handy values (all calculated at enquiry polygon -
343        # if you need them at exchange polygons we can easily do that.)
344        delta_h = openings[1].stage - openings[0].stage
345        #delta_h = openings[1].depth - openings[0].depth
346        delta_z = openings[1].elevation - openings[0].elevation
347
348        # Ideas.....
349        if openings[0].depth > 0 and openings[1].depth > 0:
350            flow_rate = delta_h
351        else:
352            flow_rate = 0.0
353
354        if openings[0].depth > self.height:
355            # This could be usefull mayby
356            #print 'Inlet has been overflowed'
357            pass
358
359        if openings[1].depth > self.height:
360            # This could be usefull mayby
361            #print 'Outlet has been overflowed'           
362            pass
363
364               
365        if delta_h > 0:
366            # Water Level U/S is higher than DS
367            flow_rate_orifice = 0.6*openings[0].area*(2*g*delta_h)**0.5      # Orifice Eqn Q= cA(2gh)^0.5
368            flow_rate_weir = 1.69*(openings[0].area)*openings[0].depth**1.5  # WEIR Eqn   Q= CLH^1.5
369
370            # Does Weir or Orifice Control Flow Rate ?
371            if flow_rate_weir > flow_rate_orifice:   
372                flow_rate_control = flow_rate_orifice
373            else:
374                flow_rate_control = flow_rate_weir
375           
376        elif delta_h < 0:
377            # Water Level D/S is higher than US
378            # That is reverse flow in culvert
379            flow_rate_orifice = 0.6*openings[0].area*(2*g*-delta_h)**0.5     # Orifice Eqn Q= cA(2gh)^0.5
380            flow_rate_weir = 1.69*openings[0].area*openings[0].depth**1.5    # WEIR Eqn   Q= CLH^1.5
381
382            # Does Weir or Orifice Control Flow Rate ?
383            if flow_rate_weir > flow_rate_orifice: 
384                flow_rate_control = flow_rate_orifice
385            else:
386                flow_rate_control = flow_rate_weir
387
388        if openings[0].depth < epsilon:
389            # Do nothing because No water over Opening  That is Set Flow to Zero!!
390            self.openings[0].rate = 0
391            self.openings[1].rate = 0
392 
393        elif openings[0].depth > 0:
394            # Only calculate fllow if there is some depth over the inlet
395            #if delta_h > 0:
396            if 1:
397                # FIXME (OLE): We never get in here. Why?
398               
399                #print 'We got water at inlet and dh > 0'
400           
401                # Flow will go from opening 0 to opening 1
402                #- that is abstract from [0] and add to [1]
403                self.openings[0].rate = -flow_rate_control
404                self.openings[1].rate = flow_rate_control
405
406                # Add jet in the form of absolute momentum to opening 1
407                #speed = 20.0 # m/s
408                # This Should be Based on the VELOCITY in the Culvert
409                # Based on the Flow Depth in the Culvert for part full
410                # or Flow Jet of the Pressurised Culvert if FUll
411                # If Part Full Flow Calculate Part Full Depth
412                # Based on Depth Calculate Area... the Vel = flow_rate_control / Area
413                # Need to Break Velocity into X & Y Components
414                #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum')
415                # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter
416               
417                #if (ratio*self.area)== 0: # Don't Really need this already established water depth here
418                #   outlet_vel=0.0
419                #else:
420
421
422                # FIXME (Ole): Shouldn't this be openings[1]??
423                outlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area))
424                outlet_dep = 2.0*openings[0].depth*openings[0].ratio #?????
425                outlet_mom = outlet_vel*outlet_dep
426
427               
428                # Eventually will Need Momentum in X & Y Components based on the orientation of
429                # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine
430                # YES - use self.vector which is a unit vector in the direction of the culvert.
431                   
432                outlet_mom_x, outlet_mom_y = self.vector * outlet_mom
433                #print 'Directional momentum', outlet_mom_x, outlet_mom_y
434                   
435                # Update the momentum forcing terms with directional momentum at the outlet
436                delta_t = time - self.last_time
437                if delta_t > 0.0:
438                    xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
439                    xmomentum_rate /= delta_t
440                   
441                    ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
442                    ymomentum_rate /= delta_t                       
443                else:
444                    xmomentum_rate = ymomentum_rate = 0.0
445
446                # Set momentum rates for outlet jet
447                self.xmom_forcing1.rate = xmomentum_rate
448                self.ymom_forcing1.rate = ymomentum_rate
449
450                # Remember this value for next step (IMPORTANT)   
451                self.xmom_forcing1.value = outlet_mom_x
452                self.ymom_forcing1.value = outlet_mom_y                   
453
454
455                if int(domain.time*100) % 100 == 0:
456                    s = 'T=%.2f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
457                        %(time, flow_rate_control, outlet_vel)
458                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
459                         %(outlet_dep, outlet_mom_x,outlet_mom_y)
460                    s += ' Momentum rate: (%.4f, %.4f)'\
461                         %(xmomentum_rate, ymomentum_rate)                   
462                    fid.write(s)
463                    print s
464            else:
465                # Opening 1 has the greatest depth. Therefore Reverse the Flow !!!
466                # Flow will go from opening 1 to opening 0, That is Abstract from [1] and add to [0]
467                self.openings[0].rate = flow_rate_control        # Else it will be Orifice Flow (Going US)
468                self.openings[1].rate = -flow_rate_control
469
470        # Second Else.... if water at outlet before at inlet !!!
471
472
473        self.openings[1].rate = 10
474        # Execute flow term for each opening
475        # This is where Inflow objects are evaluated and update the domain
476        for opening in self.openings:
477            opening(domain)
478
479        # Execute momentum terms
480        # This is where Inflow objects are evaluated and update the domain
481        self.xmom_forcing0(domain)
482        self.ymom_forcing0(domain)               
483        self.xmom_forcing1(domain)
484        self.ymom_forcing1(domain)       
485
486           
487
488        # Print out flows every 1 seconds
489        if int(time*100) % 100 == 0:
490            s = 'Time %.2f\n' %time
491            s += '    Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
492                 %(openings[0].depth,
493                   openings[0].area,
494                   openings[0].energy,
495                   openings[0].ratio)                   
496            s += '    Opening[1]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
497                 %(openings[1].depth,
498                   openings[1].area,
499                   openings[1].energy,
500                   openings[1].ratio)                 
501            s += '    Distance=%.2f, W=%.3f, O=%.3f, C=%.3f\n'\
502                 %(self.distance,
503                   flow_rate_weir,
504                   flow_rate_orifice,
505                   flow_rate_control)
506           
507            print s
508            fid.write(s)
509
510        # Store value of time   
511        self.last_time = time
512
513#------------------------------------------------------------------------------
514# Setup culvert inlets and outlets in current topography
515#------------------------------------------------------------------------------
516
517# Define culvert inlet and outlets
518culvert = Culvert_flow(domain,
519                       end_point0=[9.0, 2.5], 
520                       end_point1=[13.0, 2.5],
521                       width=1.00,                       
522                       verbose=True)
523                       
524domain.forcing_terms.append(culvert)
525
526
527#------------------------------------------------------------------------------
528# Setup boundary conditions
529#------------------------------------------------------------------------------
530#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
531Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
532Br = Reflective_boundary(domain)              # Solid reflective wall
533Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
534
535domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
536
537#------------------------------------------------------------------------------
538# Setup Application of  specialised forcing terms
539#------------------------------------------------------------------------------
540
541# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
542fixed_flow = Inflow(domain,
543                    rate=6.00,
544                    center=(2.1, 2.1),
545                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
546
547#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
548#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         
549
550domain.forcing_terms.append(fixed_flow)
551
552
553#------------------------------------------------------------------------------
554# Evolve system through time
555#------------------------------------------------------------------------------
556
557
558
559for t in domain.evolve(yieldstep = 0.1, finaltime = 10):
560    pass
561    #if int(domain.time*100) % 100 == 0:
562    #    domain.write_time()
563   
564    #if domain.get_time() >= 4 and tap.rate != 0.0:
565    #    print 'Turning tap off'
566    #    tap.rate = 0.0
567       
568    #if domain.get_time() >= 3 and sink.rate < 0.0:
569    #    print 'Turning drain on'
570    #    sink.rate = -0.8       
571    # Close
572
573fid.close()
574
575
576#------------------------------------------------------------------------------
577# Query output
578#------------------------------------------------------------------------------
579
580from anuga.shallow_water.data_manager import get_flow_through_cross_section
581
582swwfilename = domain.get_name()+'.sww'  # Output name from script
583print swwfilename
584
585polyline = [[17., 0.], [17., 5.]]
586
587time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
588
589from pylab import ion, plot
590ion()
591plot(time, Q)
592raw_input('done')
Note: See TracBrowser for help on using the repository browser.