source: anuga_work/development/culvert_flow/culvert_boyd_channel.py @ 5428

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

First workable version of culvert flows based on Boyd's generalisation of the US department of transportation's culvert model. A momentum 'jet' has been added.
The ANUGA version has the potential improvement of working with the total energy difference and the specific energy at the inlet rather than just the head water.

File size: 27.9 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
63def log(fid, s):
64    print s
65    fid.write(s + '\n')
66   
67
68# Open file for storing some specific results...
69fid = open('Culvert_Headwall_VarM.txt', 'w')
70
71length = 40.
72width = 5.
73
74#dx = dy = 1           # Resolution: Length of subdivisions on both axes
75#dx = dy = .5           # Resolution: Length of subdivisions on both axes
76dx = dy = .25           # Resolution: Length of subdivisions on both axes
77#dx = dy = .1           # Resolution: Length of subdivisions on both axes
78
79# OLE.... How do I refine the resolution... in the area where I have the Culvert Opening ???......
80#  Can I refine in a X & Y Range ???
81points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
82                                               len1=length, len2=width)
83domain = Domain(points, vertices, boundary)   
84domain.set_name('culv_dev_HW_Var_Mom')                 # Output name
85domain.set_default_order(2)
86domain.H0 = 0.01
87domain.tight_slope_limiters = True
88domain.set_minimum_storable_height(0.001)
89
90s='Size'+str(len(domain))
91log(fid, s)
92
93velocity_protection = 1.0e-4
94
95
96
97#------------------------------------------------------------------------------
98# Setup initial conditions
99#------------------------------------------------------------------------------
100
101# Define the topography (land scape)
102def topography(x, y):
103    """Set up a weir
104   
105    A culvert will connect either side of an Embankment with a Headwall type structure
106    The aim is for the Model to REALISTICALY model flow through the Culvert
107    """
108    # General Slope of Topography
109    z=-x/100
110    floorhgt = 5
111    embank_hgt=10
112    embank_upslope=embank_hgt/5
113    embank_dnslope=embank_hgt/2.5
114    # Add bits and Pieces to topography
115    N = len(x)
116    for i in range(N):
117
118        # Sloping Embankment Across Channel
119       
120        if 0.0 < x[i] < 7.51:
121            z[i]=z[i]+5.0
122        if 7.5 < x[i] < 10.1:
123            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
124               z[i]=z[i]+5.0
125            else:
126               z[i] +=  embank_upslope*(x[i] -5.0)    # Sloping Segment  U/S Face
127        if 10.0 < x[i] < 12.1:
128            if  2.0 < y[i]  < 3.0: # Cut Out Segment for Culvert (open Channel)
129                #z[i] +=  z[i]+5-(x[i]-10)*2              # Sloping Channel in Embankment
130                z[i] +=  embank_hgt                 # Flat Crest of Embankment
131            else:
132                z[i] +=  embank_hgt                 # Flat Crest of Embankment
133        if 12.0 < x[i] < 14.5:
134            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
135               z[i]=z[i]
136            else:
137               z[i] += embank_hgt-embank_dnslope*(x[i] -12.0)       # Sloping D/S Face
138     
139       
140        # Constriction
141        #if 27 < x[i] < 29 and y[i] > 3:
142        #    z[i] += 2       
143       
144        # Pole
145        #if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
146        #    z[i] += 2
147
148        # HOLE For Pit at Opening[0]
149        #if (x[i]-4)**2 + (y[i]-1.5)**2<0.75**2:
150        #  z[i]-=1
151
152        # HOLE For Pit at Opening[1]
153        #if (x[i]-20)**2 + (y[i]-3.5)**2<0.5**2:
154        #  z[i]-=1
155
156    return z
157
158
159domain.set_quantity('elevation', topography)  # Use function for elevation
160domain.set_quantity('friction', 0.01)         # Constant friction
161domain.set_quantity('stage',
162                    expression='elevation')   # Dry initial condition
163
164
165
166
167#------------------------------------------------------------------------------
168# Setup specialised forcing terms
169#------------------------------------------------------------------------------
170
171 #   NEW DEFINED CULVERT FLOW---- Flow from INLET 1 ------> INLET 2 (Outlet)
172 #
173 # The First Attempt has a Simple Horizontal Circle as a Hole on the Bed
174 # Flow Is Removed at a Rate of INFLOW
175 # Downstream there is a similar Circular Hole on the Bed where INFLOW effectively Surcharges
176 #
177 # This SHould be changed to a Vertical Opening Both BOX and Circular
178 # There will be several Culvert Routines such as:
179 # CULVERT_Boyd_Channel
180 # CULVERT_Orifice_and_Weir
181 # CULVERT_Simple_FLOOR
182 # CULVERT_Simple_WALL
183 # CULVERT_Eqn_Floor
184 # CULVERT_Eqn_Wall
185 # CULVERT_Tab_Floor
186 # CULVERT_Tab_Wall
187 # BRIDGES.....
188 # NOTE NEED TO DEVELOP CONCEPT 1D Model for Linked Pipe System !!!!
189 
190 #  COULD USE EPA SWMM Model !!!!
191 
192 
193 
194class Culvert_flow:
195    """Culvert flow - transfer water from one hole to another
196   
197    Using Momentum as Calculated by Culvert Flow !!
198    Could be Several Methods Investigated to do This !!!
199
200    2008_May_08
201    To Ole:
202    OK so here we need to get the Polygon Creating code to create polygons for the culvert Based on
203    the two input Points (X0,Y0) and (X1,Y1) - need to be passed to create polygon
204
205    The two centers are now passed on to create_polygon.
206   
207
208    Input: Two points, pipe_size (either diameter or width, height), mannings_rougness,
209    inlet/outlet energy_loss_coefficients, internal_bend_coefficent,
210    top-down_blockage_factor and bottom_up_blockage_factor
211   
212   
213    And the Delta H enquiry should be change from Openings in line 412 to the enquiry Polygons infront
214    of the culvert
215    At the moment this script uses only Depth, later we can change it to Energy...
216
217    Once we have Delta H can calculate a Flow Rate and from Flow Rate an Outlet Velocity
218    The Outlet Velocity x Outlet Depth = Momentum to be applied at the Outlet...
219       
220    """ 
221
222    def __init__(self,
223                 domain,
224                 label=None, 
225                 description=None,
226                 end_point0=None, 
227                 end_point1=None,
228                 width=None,
229                 height=None,
230                 manning=None,          # Mannings Roughness for Culvert
231                 invert_level0=None,    # Invert level if not the same as the Elevation on the Domain
232                 invert_level1=None,    # Invert level if not the same as the Elevation on the Domain
233                 loss_exit=None,
234                 loss_entry=None,
235                 loss_bend=None,
236                 loss_special=None,
237                 blockage_topdwn=None,
238                 blockage_bottup=None,
239                 verbose=False):
240       
241        from Numeric import sqrt, sum
242
243        # Input check
244        if height is None: height = width
245
246        # Set defaults
247        if manning is None: manning = 0.012   # Set a Default Mannings Roughness for Pipe
248        if loss_exit is None: loss_exit = 1.00
249        if loss_entry is None: loss_entry = 0.50
250        if loss_bend is None: loss_bend=0.00
251        if loss_special is None: loss_special=0.00
252        if blockage_topdwn is None: blockage_topdwn=0.00
253        if blockage_bottup is None: blockage_bottup=0.00
254       
255
256        # Create the fundamental culvert polygons from POLYGON
257        P = create_culvert_polygons(end_point0,
258                                    end_point1,
259                                    width=width,   
260                                    height=height)
261       
262        if verbose is True:
263            pass
264            #plot_polygons([[end_point0, end_point1],
265            #               P['exchange_polygon0'],
266            #               P['exchange_polygon1'],
267            #               P['enquiry_polygon0'],
268            #               P['enquiry_polygon1']],
269            #              figname='culvert_polygon_output')
270       
271        self.openings = []
272        self.openings.append(Inflow(domain,
273                                    polygon=P['exchange_polygon0']))
274
275        self.openings.append(Inflow(domain,
276                                    polygon=P['exchange_polygon1']))                                   
277
278
279        # Assume two openings for now: Referred to as 0 and 1
280        assert len(self.openings) == 2
281       
282        # Store basic geometry
283        self.end_points = [end_point0, end_point1]
284        self.invert_levels = [invert_level0, invert_level1]               
285        self.enquiry_polygons = [P['enquiry_polygon0'], P['enquiry_polygon1']]
286        self.vector = P['vector']
287        self.distance = P['length']; assert self.distance > 0.0
288        self.verbose = verbose
289        self.width = width
290        self.height = height
291        self.last_time = 0.0       
292        self.temp_keep_delta_t = 0.0
293       
294
295        # Store hydraulic parameters
296        self.manning = manning
297        self.loss_exit = loss_exit
298        self.loss_entry = loss_entry
299        self.loss_bend = loss_bend
300        self.loss_special = loss_special
301        self.sum_loss = loss_exit + loss_entry + loss_bend + loss_special
302        self.blockage_topdwn = blockage_topdwn
303        self.blockage_bottup = blockage_bottup
304
305       
306        # Create objects to update momentum (a bit crude at this stage)
307
308       
309        xmom0 = General_forcing(domain, 'xmomentum',
310                                polygon=P['exchange_polygon0'])
311
312        xmom1 = General_forcing(domain, 'xmomentum',
313                                polygon=P['exchange_polygon1'])
314
315        ymom0 = General_forcing(domain, 'ymomentum',
316                                polygon=P['exchange_polygon0'])
317
318        ymom1 = General_forcing(domain, 'ymomentum',
319                                polygon=P['exchange_polygon1'])
320
321        self.opening_momentum = [ [xmom0, ymom0], [xmom1, ymom1] ]
322       
323
324        # Print something
325        s = 'Culvert Effective Length = %.2f m' %(self.distance)
326        log(fid, s)
327   
328        s = 'Culvert Direction is %s\n' %str(self.vector)
329        log(fid, s)       
330       
331    def __call__(self, domain):
332        from anuga.utilities.numerical_tools import mean
333        from anuga.utilities.polygon import inside_polygon
334        from anuga.config import g, epsilon
335        from Numeric import take, sqrt
336
337         
338        # Time stuff
339        time = domain.get_time()
340        delta_t = time-self.last_time
341        s = '\nTime = %.2f, delta_t = %f' %(time, delta_t)
342        log(fid, s)
343       
344        msg = 'Time did not advance'
345        if time > 0.0: assert delta_t > 0.0, msg
346
347
348        # Get average water depths at each opening       
349        openings = self.openings   # There are two Opening [0] and [1]
350        for i, opening in enumerate(openings):
351            stage = domain.quantities['stage'].get_values(location='centroids',
352                                                          indices=opening.exchange_indices)
353            elevation = domain.quantities['elevation'].get_values(location='centroids',
354                                                                  indices=opening.exchange_indices)
355
356            # Indices corresponding to energy enquiry field for this opening
357            coordinates = domain.get_centroid_coordinates() # Get all centroid points (x,y)
358            enquiry_indices = inside_polygon(coordinates, self.enquiry_polygons[i]) 
359           
360            if self.verbose:
361                pass
362                #print 'Opening %d: Got %d points in enquiry polygon:\n%s'\
363                #      %(i, len(idx), self.enquiry_polygons[i])
364               
365           
366            # Get model values for points in enquiry polygon for this opening
367            dq = domain.quantities
368            stage = dq['stage'].get_values(location='centroids', indices=enquiry_indices)
369            xmomentum = dq['xmomentum'].get_values(location='centroids', indices=enquiry_indices)
370            ymomentum = dq['ymomentum'].get_values(location='centroids', indices=enquiry_indices)                       
371            elevation = dq['elevation'].get_values(location='centroids', indices=enquiry_indices)
372            depth = stage - elevation
373
374            # Compute mean values of selected quantitites in the enquiry area in front of the culvert
375            # Epsilon handles a dry cell case
376            ux = xmomentum/(depth+velocity_protection/depth)   # Velocity (x-direction)
377            uy = ymomentum/(depth+velocity_protection/depth)   # Velocity (y-direction)
378            v = mean(sqrt(ux**2+uy**2))      # Average velocity
379            w = mean(stage)                  # Average stage
380
381            # Store values at enquiry field
382            opening.velocity = v
383
384
385            # Compute mean values of selected quantitites in the exchange area in front of the culvert
386            # Stage and velocity comes from enquiry area and elevation from exchange area
387           
388            # Use invert level instead of elevation if specified
389            invert_level = self.invert_levels[i]
390            if invert_level is not None:
391                z = invert_level
392            else:
393                elevation = dq['elevation'].get_values(location='centroids', indices=opening.exchange_indices)
394                z = mean(elevation)                   # Average elevation
395               
396            # Estimated depth above the culvert inlet
397            d = w - z
398
399            if d < 0.0:
400                # This is possible since w and z are taken at different locations
401                #msg = 'D < 0.0: %f' %d
402                #raise msg
403                d = 0.0
404           
405            # Ratio of depth to culvert height.
406            # If ratio > 1 then culvert is running full
407            ratio = d/self.height 
408            opening.ratio = ratio
409               
410            # Average measures of energy in front of this opening
411            Es = d + 0.5*v**2/#  Specific energy in exchange area
412            Et = w + 0.5*v**2/#  Total energy in the enquiry area
413            opening.total_energy = Et
414            opening.specific_energy = Es           
415           
416            # Store current average stage and depth with each opening object
417            opening.depth = d
418            opening.stage = w
419            opening.elevation = z
420           
421
422        #################  End of the FOR loop ################################################
423
424           
425        # We now need to deal with each opening individually
426
427        # Determine flow direction based on total energy difference
428        delta_Et = openings[0].total_energy - openings[1].total_energy
429
430        if delta_Et > 0:
431            #print 'Flow U/S ---> D/S'
432            inlet=openings[0]
433            outlet=openings[1]
434
435            inlet.momentum = self.opening_momentum[0]
436            outlet.momentum = self.opening_momentum[1]           
437        else:
438            #print 'Flow D/S ---> U/S'
439            inlet=openings[1]
440            outlet=openings[0]
441
442            inlet.momentum = self.opening_momentum[1]
443            outlet.momentum = self.opening_momentum[0]
444           
445            delta_Et = -delta_Et
446
447        msg = 'Total energy difference is negative'
448        assert delta_Et > 0.0, msg
449
450        delta_h = inlet.stage - outlet.stage
451        delta_z = inlet.elevation - outlet.elevation
452        culvert_slope = (delta_z/self.distance)
453
454        if culvert_slope < 0.0:
455            # Adverse gradient - flow is running uphill
456            # Flow will be purely controlled by uphill outlet face
457            print 'WARNING: Flow is running uphill. Watch Out!', inlet.elevation, outlet.elevation
458
459
460        s = 'Delta total energy = %.3f' %(delta_Et)
461        log(fid, s)
462       
463
464        # ====================== NOW ENTER INTO THE CULVERT EQUATIONS AS DEFINED BY BOYD GENERALISED METHOD
465        # == The quantity of flow passing through a culvert is controlled by many factors
466        # == It could be that the culvert is controled by the inlet only, with it being Un submerged this is effectively equivalent to the WEIR Equation
467        # == Else the culvert could be controlled by the inlet, with it being Submerged, this is effectively the Orifice Equation
468        # == Else it may be controlled by Down stream conditions where depending on the down stream depth, the momentum in the culvert etc. flow is controlled
469       
470       
471        Q_outlet_tailwater = 0.0
472        inlet.rate = 0.0
473        outlet.rate = 0.0
474        Q_inlet_unsubmerged = 0.0
475        Q_inlet_submerged = 0.0
476        Q_outlet_critical_depth = 0.0
477       
478        if inlet.depth >= 0.01: 
479            # Water has risen above inlet
480            if self.width==self.height:  # Need something better to Flag Diamater !!!!!!! Can't we just have DIAMETER as well as Width & Height ???
481                pass
482                #Q1[t]= 0.421*g**0.5*Diam**0.87*HW**1.63    # Inlet Ctrl Inlet Unsubmerged
483                #Q2[t]= 0.530*g**0.5*Diam**1.87*HW**0.63    # Inlet Ctrl Inlet Submerged
484                #PipeDcrit=
485                #Delta_E=HW-TW
486            else:
487                # Box culvert
488               
489                sum_loss=self.loss_exit+self.loss_entry+self.loss_bend+self.loss_special
490                manning=self.manning
491
492                # Calculate flows for inlet control
493                E = inlet.specific_energy
494                #E = min(inlet.specific_energy, delta_Et)
495
496                s = 'Driving energy  = %f m' %E
497                log(fid, s)
498
499                msg = 'Driving energy is negative'
500                assert E >= 0.0, msg
501                         
502                Q_inlet_unsubmerged = 0.540*g**0.5*self.width*E**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
503                Q_inlet_submerged = 0.702*g**0.5*self.width*self.height**0.89*E**0.61  # Flow based on Inlet Ctrl Inlet Submerged
504
505                s = 'Q_inlet_unsubmerged = %.6f, Q_inlet_submerged = %.6f' %(Q_inlet_unsubmerged, Q_inlet_submerged)
506                log(fid, s)
507
508                case = ''
509                if Q_inlet_unsubmerged < Q_inlet_submerged:
510                    Q = Q_inlet_unsubmerged
511                    flow_area = self.width*inlet.depth
512                    outlet_culv_depth = inlet.depth
513                    case = 'Inlet unsubmerged'
514                else:   
515                    Q = Q_inlet_submerged
516                    flow_area = self.width*self.height
517                    outlet_culv_depth = self.height
518                    case = 'Inlet submerged'                   
519
520                if delta_Et < Es:
521                    # Calculate flows for outlet control
522                    # Determine the depth at the outlet relative to the depth of flow in the Culvert
523
524                    sum_loss = self.sum_loss
525               
526                    if outlet.depth > self.height:  # The Outlet is Submerged
527                        outlet_culv_depth=self.height
528                        flow_area=self.width*self.height # Cross sectional area of flow in the culvert
529                        perimeter=2.0*(self.width+self.height)
530                        case = 'Outlet submerged'
531                    elif outlet.depth==0.0: 
532                        outlet_culv_depth=inlet.depth   # For purpose of calculation assume the outlet depth = the inlet depth
533                        flow_area=self.width*inlet.depth
534                        perimeter=(self.width+2.0*inlet.depth)
535                        case = 'Outlet depth is zero'                       
536                    else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
537                        outlet_culv_depth=outlet.depth
538                        flow_area=self.width*outlet.depth
539                        perimeter=(self.width+2.0*outlet.depth)
540                        case = 'Outlet is open channel flow'
541                       
542                    hyd_rad = flow_area/perimeter
543                    s = 'hydraulic radius at outlet = %f' %hyd_rad
544                    log(fid, s)
545                   
546
547                    # Outlet control velocity using tail water
548                    culvert_velocity = sqrt(delta_Et/((sum_loss/2*g)+(manning**2*self.distance)/hyd_rad**1.33333)) 
549                    Q_outlet_tailwater = flow_area * culvert_velocity
550
551                    s = 'Q_outlet_tailwater = %.6f' %Q_outlet_tailwater
552                    log(fid, s)
553
554                    Q = min(Q, Q_outlet_tailwater)
555
556
557                   
558                log(fid, 'Case: "%s"' %case)
559               
560                flow_rate_control=Q
561
562                s = 'Flow Rate Control = %f' %flow_rate_control
563                log(fid, s)
564
565                inlet.rate = -flow_rate_control
566                outlet.rate = flow_rate_control               
567               
568                culv_froude=sqrt(flow_rate_control**2*self.width/(g*flow_area**3))
569                s = 'Froude in Culvert = %f' %culv_froude
570                log(fid, s)
571
572               
573
574                # Determine momentum at the outlet
575                barrel_velocity = Q/(flow_area + velocity_protection/flow_area)
576                barrel_momentum = barrel_velocity*outlet_culv_depth
577               
578                outlet_E_Loss= 0.5*0.5*barrel_velocity**2/g   # Ke v^2/2g
579                s = 'Barrel velocity = %f' %barrel_velocity
580                log(fid, s)
581
582                # Compute momentum vector at outlet
583                outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum
584               
585                s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)
586                log(fid, s)
587
588                delta_t = time - self.last_time
589                if delta_t > 0.0:
590                    xmomentum_rate = outlet_mom_x - outlet.momentum[0].value
591                    xmomentum_rate /= delta_t
592                   
593                    ymomentum_rate = outlet_mom_y - outlet.momentum[1].value
594                    ymomentum_rate /= delta_t
595                   
596                    s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)
597                    log(fid, s)                   
598                else:
599                    xmomentum_rate = ymomentum_rate = 0.0
600
601
602                # Set momentum rates for outlet jet
603                outlet.momentum[0].rate = xmomentum_rate
604                outlet.momentum[1].rate = ymomentum_rate
605
606                # Remember this value for next step (IMPORTANT)
607                outlet.momentum[0].value = outlet_mom_x
608                outlet.momentum[1].value = outlet_mom_y                   
609
610                if int(domain.time*100) % 100 == 0:
611                    s = 'T=%.5f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
612                        %(time, flow_rate_control, barrel_velocity)
613                    s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
614                         %(outlet_culv_depth, outlet_mom_x,outlet_mom_y)
615                    s += ' Momentum rate: (%.4f, %.4f)'\
616                         %(xmomentum_rate, ymomentum_rate)                   
617                    s+='Outlet Vel= %.3f'\
618                         %(barrel_velocity)
619                    log(fid, s)
620         
621       
622        # Execute flow term for each opening
623        # This is where Inflow objects are evaluated and update the domain
624        for opening in self.openings:
625            opening(domain)
626
627        # Execute momentum terms
628        # This is where Inflow objects are evaluated and update the domain
629        outlet.momentum[0](domain)
630        outlet.momentum[1](domain)       
631           
632        # Store value of time   
633        self.last_time = time
634
635
636
637
638#------------------------------------------------------------------------------
639# Setup culvert inlets and outlets in current topography
640#------------------------------------------------------------------------------
641
642# Define culvert inlet and outlets
643# NEED TO ADD Mannings N as Fixed Value or Function
644# Energy Loss Coefficients as Fixed or Function
645# Also Set the Shape & Gap Factors fo rthe Enquiry PolyGons
646# ALSO Allow the Invert Level to be provided by the USER
647culvert = Culvert_flow(domain,
648                       label='Culvert No. 1',
649                       description=' This culvert is a test unit 1.2m Wide by 0.75m High',   
650                       end_point0=[9.0, 2.5], 
651                       end_point1=[13.0, 2.5],
652                       width=1.20,height=0.75,                       
653                       verbose=True)
654                       
655domain.forcing_terms.append(culvert)
656
657
658#------------------------------------------------------------------------------
659# Setup boundary conditions
660#------------------------------------------------------------------------------
661#Bi = Dirichlet_boundary([0.5, 0.0, 0.0])          # Inflow based on Flow Depth (0.5m) and Approaching Momentum !!!
662Bi = Dirichlet_boundary([0.0, 0.0, 0.0])          # Inflow based on Flow Depth and Approaching Momentum !!!
663Br = Reflective_boundary(domain)              # Solid reflective wall
664Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
665
666domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
667
668#------------------------------------------------------------------------------
669# Setup Application of  specialised forcing terms
670#------------------------------------------------------------------------------
671
672# This is the new element implemented by Ole to allow direct input of Inflow in m^3/s
673fixed_flow = Inflow(domain,
674                    rate=20.00,
675                    center=(2.1, 2.1),
676                    radius=1.261566)   #   Fixed Flow Value Over Area of 5m2 at 1m/s = 5m^3/s
677
678#            flow=file_function('Q/QPMF_Rot_Sub13.tms'))        # Read Time Series in  from File
679#             flow=lambda t: min(0.01*t, 0.01942))                             # Time Varying Function   Tap turning up         
680
681domain.forcing_terms.append(fixed_flow)
682
683
684#------------------------------------------------------------------------------
685# Evolve system through time
686#------------------------------------------------------------------------------
687
688temp_keep_delta_t=0.0
689
690for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
691    pass
692    #if int(domain.time*100) % 100 == 0:
693    #    domain.write_time()
694   
695    #if domain.get_time() >= 4 and tap.flow != 0.0:
696    #    print 'Turning tap off'
697    #    tap.flow = 0.0
698       
699    #if domain.get_time() >= 3 and sink.flow < 0.0:
700    #    print 'Turning drain on'
701    #    sink.flow = -0.8       
702    # Close
703
704fid.close()
705
706
707#------------------------------------------------------------------------------
708# Query output
709#------------------------------------------------------------------------------
710
711from anuga.shallow_water.data_manager import get_flow_through_cross_section
712
713swwfilename = domain.get_name()+'.sww'  # Output name from script
714print swwfilename
715
716polyline = [[17., 0.], [17., 5.]]
717
718time, Q = get_flow_through_cross_section(swwfilename, polyline, verbose=True)
719
720from pylab import ion, plot
721ion()
722plot(time, Q)
723raw_input('done')
Note: See TracBrowser for help on using the repository browser.