# source:anuga_work/development/Rudy_vandrie_2007/culvert_development.py@5406

Last change on this file since 5406 was 5406, checked in by ole, 15 years ago

Added time print out at every timestep (as diagnostics)

File size: 24.5 KB
RevLine
[5305]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
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...
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()
[5406]286
287        print 'Time = %f, delta_t = %f' %(time, time-self.last_time)
288
[5305]289        openings = self.openings
290        for i, opening in enumerate(openings):
[5340]291            # Quantities corresponding to fluid exchange field for this opening
292            stage_o = domain.quantities['stage'].get_values(location='centroids',
[5305]293                                                          indices=opening.indices)
[5340]294            elevation_o = domain.quantities['elevation'].get_values(location='centroids',
[5305]295                                                                  indices=opening.indices)
[5340]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)
[5305]300
[5340]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
[5305]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)
[5340]339            self.mean_xmomentum = mean(xmomentum)
340            self.mean_ymomentum = mean(ymomentum)
341
[5305]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/# 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!!
[5314]412            self.openings[0].rate = 0
413            self.openings[1].rate = 0
[5305]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]
[5314]425                self.openings[0].rate = -flow_rate_control
426                self.openings[1].rate = flow_rate_control
[5305]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:
[5340]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
[5305]465
[5340]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
[5305]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
[5340]483                if int(time*100) % 100 == 0:
484                    s = 'T=%.3f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
[5305]485                        %(time, flow_rate_control, outlet_vel)
[5340]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)'\
[5314]493                         %(xmomentum_rate, ymomentum_rate)
[5305]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]
[5314]499                self.openings[0].rate = flow_rate_control        # Else it will be Orifice Flow (Going US)
500                self.openings[1].rate = -flow_rate_control
[5305]501
502        # Second Else.... if water at outlet before at inlet !!!
503
504
[5314]505        self.openings[1].rate = 10
[5305]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:
[5340]522            s = 'Time %.3f\n' %time
[5305]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
[5340]540
[5305]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
551culvert = 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
557domain.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 !!!
564Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
565Br = Reflective_boundary(domain)              # Solid reflective wall
566Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
567
568domain.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
575fixed_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
583domain.forcing_terms.append(fixed_flow)
584
585
586#------------------------------------------------------------------------------
587# Evolve system through time
588#------------------------------------------------------------------------------
589
590
591
592for t in domain.evolve(yieldstep = 0.1, finaltime = 10):
593    pass
594    #if int(domain.time*100) % 100 == 0:
595    #    domain.write_time()
596
[5314]597    #if domain.get_time() >= 4 and tap.rate != 0.0:
[5305]598    #    print 'Turning tap off'
[5314]599    #    tap.rate = 0.0
[5305]600
[5314]601    #if domain.get_time() >= 3 and sink.rate < 0.0:
[5305]602    #    print 'Turning drain on'
[5314]603    #    sink.rate = -0.8
[5305]604    # Close
605
606fid.close()
607
608
609#------------------------------------------------------------------------------
610# Query output
611#------------------------------------------------------------------------------
612
613from anuga.shallow_water.data_manager import get_flow_through_cross_section
614
615swwfilename = domain.get_name()+'.sww'  # Output name from script
616print swwfilename
617
618polyline = [[17., 0.], [17., 5.]]
619
620time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
621
622from pylab import ion, plot
623ion()
624plot(time, Q)
625raw_input('done')
Note: See TracBrowser for help on using the repository browser.