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

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

Work on momentum jet

File size: 24.4 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            # Quantities corresponding to fluid exchange field for this opening           
289            stage_o = domain.quantities['stage'].get_values(location='centroids',
290                                                          indices=opening.indices)
291            elevation_o = domain.quantities['elevation'].get_values(location='centroids',
292                                                                  indices=opening.indices)
293            xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
294                                                                  indices=opening.indices)
295            ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
296                                                                  indices=opening.indices)
297
298            # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account)
299            depth_o = stage_o - elevation_o
300           
301            ux_o = xmomentum_o/(depth_o+epsilon)
302            uy_o = ymomentum_o/(depth_o+epsilon)
303            v_o = mean(sqrt(ux_o**2+uy_o**2))
304            d_o = mean(depth_o)
305            w_o = mean(stage_o)
306            z_o = mean(elevation_o)
307            self.mean_xmomentum_o = mean(xmomentum_o)
308            self.mean_ymomentum_o = mean(ymomentum_o)             
309
310            # Indices corresponding to energy enquiry field for this opening
311            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
312            idx = inside_polygon(coordinates, self.enquiry_polygons[i]) 
313           
314            if self.verbose:
315                pass
316                #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
317                #      %(i, len(idx), self.enquiry_polygons[i])
318               
319           
320            # Get average model values for points in
321            # enquiry polygon for this opening
322            dq = domain.quantities
323            stage = dq['stage'].get_values(location='centroids', indices=idx)
324            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=idx)
325            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=idx)                       
326            elevation = dq['elevation'].get_values(location='centroids', indices=idx)
327            depth = stage - elevation           
328
329            # Compute mean velocity in the area in front of the culvert (taking zero depths into account)
330            ux = xmomentum/(depth+epsilon)
331            uy = ymomentum/(depth+epsilon)
332            v = mean(sqrt(ux**2+uy**2))
333            d = mean(depth)
334            w = mean(stage)
335            z = mean(elevation)
336            self.mean_xmomentum = mean(xmomentum)
337            self.mean_ymomentum = mean(ymomentum)             
338
339            # Ratio of depth to culvert height
340            ratio = d/(2*self.height)
341            if ratio > 1.0:    # Assume culvert is running full &
342                ratio = 1.0    # under pressure. Note this is usually ~ 1.35
343            opening.ratio = ratio
344               
345
346            # Average measure of total energy (D + V^2/2g) in enquiry field in front of this opening
347            E = d + 0.5*v**2/# RUDY - please check this
348            opening.energy = E
349           
350           
351
352            # Store current average stage and depth at enquiry field with each opening object
353            opening.depth = d
354            opening.stage = w
355            opening.elevation = z
356           
357
358        # Now we are done calculating energy etc for each opening.
359        # At this point we can work on the transfer functions etc.
360       
361        # Handy values (all calculated at enquiry polygon -
362        # if you need them at exchange polygons we can easily do that.)
363        delta_h = openings[1].stage - openings[0].stage
364        #delta_h = openings[1].depth - openings[0].depth
365        delta_z = openings[1].elevation - openings[0].elevation
366
367        # Ideas.....
368        if openings[0].depth > 0 and openings[1].depth > 0:
369            flow_rate = delta_h
370        else:
371            flow_rate = 0.0
372
373        if openings[0].depth > self.height:
374            # This could be usefull mayby
375            #print 'Inlet has been overflowed'
376            pass
377
378        if openings[1].depth > self.height:
379            # This could be usefull mayby
380            #print 'Outlet has been overflowed'           
381            pass
382
383               
384        if delta_h > 0:
385            # Water Level U/S is higher than DS
386            flow_rate_orifice = 0.6*openings[0].area*(2*g*delta_h)**0.5      # Orifice Eqn Q= cA(2gh)^0.5
387            flow_rate_weir = 1.69*(openings[0].area)*openings[0].depth**1.5  # WEIR Eqn   Q= CLH^1.5
388
389            # Does Weir or Orifice Control Flow Rate ?
390            if flow_rate_weir > flow_rate_orifice:   
391                flow_rate_control = flow_rate_orifice
392            else:
393                flow_rate_control = flow_rate_weir
394           
395        elif delta_h < 0:
396            # Water Level D/S is higher than US
397            # That is reverse flow in culvert
398            flow_rate_orifice = 0.6*openings[0].area*(2*g*-delta_h)**0.5     # Orifice Eqn Q= cA(2gh)^0.5
399            flow_rate_weir = 1.69*openings[0].area*openings[0].depth**1.5    # WEIR Eqn   Q= CLH^1.5
400
401            # Does Weir or Orifice Control Flow Rate ?
402            if flow_rate_weir > flow_rate_orifice: 
403                flow_rate_control = flow_rate_orifice
404            else:
405                flow_rate_control = flow_rate_weir
406
407        if openings[0].depth < epsilon:
408            # Do nothing because No water over Opening  That is Set Flow to Zero!!
409            self.openings[0].rate = 0
410            self.openings[1].rate = 0
411 
412        elif openings[0].depth > 0:
413            # Only calculate fllow if there is some depth over the inlet
414            #if delta_h > 0:
415            if 1:
416                # FIXME (OLE): We never get in here. Why?
417               
418                #print 'We got water at inlet and dh > 0'
419           
420                # Flow will go from opening 0 to opening 1
421                #- that is abstract from [0] and add to [1]
422                self.openings[0].rate = -flow_rate_control
423                self.openings[1].rate = flow_rate_control
424
425                # Add jet in the form of absolute momentum to opening 1
426                #speed = 20.0 # m/s
427                # This Should be Based on the VELOCITY in the Culvert
428                # Based on the Flow Depth in the Culvert for part full
429                # or Flow Jet of the Pressurised Culvert if FUll
430                # If Part Full Flow Calculate Part Full Depth
431                # Based on Depth Calculate Area... the Vel = flow_rate_control / Area
432                # Need to Break Velocity into X & Y Components
433                #self.openings[1].set_quantity_values(delta_h*speed, 'xmomentum')
434                # Previous Calculated Depth/ Culvert Heigh Ratio Use it to Determine Velocity !!!! FIX LAter
435               
436                #if (ratio*self.area)== 0: # Don't Really need this already established water depth here
437                #   outlet_vel=0.0
438                #else:
439
440
441                # FIXME (Ole): Shouldn't this be openings[1]??
442                outlet_vel =(flow_rate_control/(openings[0].ratio*openings[0].area))
443                outlet_dep = 2.0*openings[0].depth*openings[0].ratio #?????
444                outlet_mom = outlet_vel*outlet_dep
445
446               
447                # Eventually will Need Momentum in X & Y Components based on the orientation of
448                # the culvert from X0,Y0, X1, Y1 from Create Polygon Routine
449                # YES - use self.vector which is a unit vector in the direction of the culvert.
450                   
451                outlet_mom_x, outlet_mom_y = self.vector * outlet_mom
452                #print 'Directional momentum', outlet_mom_x, outlet_mom_y
453                   
454                # Update the momentum forcing terms with directional momentum at the outlet
455                delta_t = time - self.last_time
456                if delta_t > 0.0:
457                    xmomentum_rate = outlet_mom_x - self.mean_xmomentum
458                    if xmomentum_rate > 0:
459                        xmomentum_rate /= delta_t
460                    else:
461                        xmomentum_rate = 0.0
462                   
463                    ymomentum_rate = outlet_mom_y - self.mean_ymomentum
464                    if ymomentum_rate > 0:
465                        ymomentum_rate /= delta_t
466                    else:
467                        ymomentum_rate = 0.0                   
468                else:
469                    xmomentum_rate = ymomentum_rate = 0.0
470
471                # Set momentum rates for outlet jet
472                self.xmom_forcing1.rate = xmomentum_rate
473                self.ymom_forcing1.rate = ymomentum_rate
474
475                # Remember this value for next step (IMPORTANT)   
476                self.xmom_forcing1.value = outlet_mom_x
477                self.ymom_forcing1.value = outlet_mom_y                   
478
479
480                if int(time*100) % 100 == 0:
481                    s = 'T=%.3f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
482                        %(time, flow_rate_control, outlet_vel)
483                    s += ' Depth= %0.3f\n    Outlet Momentum = (%0.3f, %0.3f)\n'\
484                         %(outlet_dep, outlet_mom_x, outlet_mom_y)
485                    s += '    Avg Momentum at opening = (%0.3f, %0.3f)\n'\
486                         %(self.mean_xmomentum_o, self.mean_ymomentum_o)
487                    s += '    Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\
488                         %(self.mean_xmomentum, self.mean_ymomentum)                   
489                    s += '    Momentum rate: (%.4f, %.4f)'\
490                         %(xmomentum_rate, ymomentum_rate)                   
491                    fid.write(s)
492                    print s
493            else:
494                # Opening 1 has the greatest depth. Therefore Reverse the Flow !!!
495                # Flow will go from opening 1 to opening 0, That is Abstract from [1] and add to [0]
496                self.openings[0].rate = flow_rate_control        # Else it will be Orifice Flow (Going US)
497                self.openings[1].rate = -flow_rate_control
498
499        # Second Else.... if water at outlet before at inlet !!!
500
501
502        self.openings[1].rate = 10
503        # Execute flow term for each opening
504        # This is where Inflow objects are evaluated and update the domain
505        for opening in self.openings:
506            opening(domain)
507
508        # Execute momentum terms
509        # This is where Inflow objects are evaluated and update the domain
510        self.xmom_forcing0(domain)
511        self.ymom_forcing0(domain)               
512        self.xmom_forcing1(domain)
513        self.ymom_forcing1(domain)       
514
515           
516
517        # Print out flows every 1 seconds
518        if int(time*100) % 100 == 0:
519            s = 'Time %.3f\n' %time
520            s += '    Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
521                 %(openings[0].depth,
522                   openings[0].area,
523                   openings[0].energy,
524                   openings[0].ratio)                   
525            s += '    Opening[1]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
526                 %(openings[1].depth,
527                   openings[1].area,
528                   openings[1].energy,
529                   openings[1].ratio)                 
530            s += '    Distance=%.2f, W=%.3f, O=%.3f, C=%.3f\n'\
531                 %(self.distance,
532                   flow_rate_weir,
533                   flow_rate_orifice,
534                   flow_rate_control)
535           
536            print s
537
538            fid.write(s)
539
540        # Store value of time   
541        self.last_time = time
542
543#------------------------------------------------------------------------------
544# Setup culvert inlets and outlets in current topography
545#------------------------------------------------------------------------------
546
547# Define culvert inlet and outlets
548culvert = Culvert_flow(domain,
549                       end_point0=[9.0, 2.5], 
550                       end_point1=[13.0, 2.5],
551                       width=1.00,                       
552                       verbose=True)
553                       
554domain.forcing_terms.append(culvert)
555
556
557#------------------------------------------------------------------------------
558# Setup boundary conditions
559#------------------------------------------------------------------------------
560#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
561Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
562Br = Reflective_boundary(domain)              # Solid reflective wall
563Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
564
565domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
566
567#------------------------------------------------------------------------------
568# Setup Application of  specialised forcing terms
569#------------------------------------------------------------------------------
570
571# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
572fixed_flow = Inflow(domain,
573                    rate=6.00,
574                    center=(2.1, 2.1),
575                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
576
577#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
578#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         
579
580domain.forcing_terms.append(fixed_flow)
581
582
583#------------------------------------------------------------------------------
584# Evolve system through time
585#------------------------------------------------------------------------------
586
587
588
589for t in domain.evolve(yieldstep = 0.1, finaltime = 10):
590    pass
591    #if int(domain.time*100) % 100 == 0:
592    #    domain.write_time()
593   
594    #if domain.get_time() >= 4 and tap.rate != 0.0:
595    #    print 'Turning tap off'
596    #    tap.rate = 0.0
597       
598    #if domain.get_time() >= 3 and sink.rate < 0.0:
599    #    print 'Turning drain on'
600    #    sink.rate = -0.8       
601    # Close
602
603fid.close()
604
605
606#------------------------------------------------------------------------------
607# Query output
608#------------------------------------------------------------------------------
609
610from anuga.shallow_water.data_manager import get_flow_through_cross_section
611
612swwfilename = domain.get_name()+'.sww'  # Output name from script
613print swwfilename
614
615polyline = [[17., 0.], [17., 5.]]
616
617time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
618
619from pylab import ion, plot
620ion()
621plot(time, Q)
622raw_input('done')
Note: See TracBrowser for help on using the repository browser.