source: anuga_work/development/ted_rigby/test_flowhydraulics_in_vee_using_momentum.py @ 7475

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

Added tests from Ted Rigby that cause unexpected flows.

File size: 6.7 KB
Line 
1"""test_flow_in_vee_using_momentum.py
2
3Test the ability of ANUGA to
4a) maintain the flow rate input into the channel at the upstream end
5b) replicate normal depth half way down the channel after half an hour
6c) replicate a typical velocity profile across the channel after half an hour
7
8The vee channel is 20m wide by 1000m long with side slopes of 1:1
9The equations for dn are generalised with side slope as 1 V in Side_Slope H
10The test location is in the invert of (y=10m)and half way down the channel (x=500m)
11This script explores flow behaviour predicted by ANUGA and compares it with known solutions
12
13This test inputs flow using momentum=flow/width .
14
15"""
16
17verbose = True
18import numpy as num
19
20#------------------------------------------------------------------------------
21# Import necessary modules
22#------------------------------------------------------------------------------
23from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
24from anuga.shallow_water import Domain
25from anuga.shallow_water.shallow_water_domain import Reflective_boundary
26from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
27from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary
28from anuga.shallow_water.shallow_water_domain import Transmissive_boundary
29from anuga.shallow_water.shallow_water_domain import Inflow
30from anuga.shallow_water.data_manager import get_flow_through_cross_section
31from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
32
33
34#------------------------------------------------------------------------------
35# Setup computational domain
36#------------------------------------------------------------------------------
37
38finaltime = 1800.0    # 1/2 hour
39
40length     = 1000.
41width      = 20.
42side_slope = 1.0
43dx = dy    = 1          # Resolution: of points (elev) grid on both axes
44ref_flow   = 20.0       # Inflow at head of vee
45slope      = 1.0/1000
46mannings_n = 0.150
47
48points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
49                                               len1=length, len2=width)
50
51
52
53domain = Domain(points, vertices, boundary)   
54domain.set_name('test_flowhydraulics_in_vee_using_momentum') # Output name
55
56
57#------------------------------------------------------------------------------
58# Setup initial conditions
59#------------------------------------------------------------------------------
60
61def topography(x, y):
62    z=-x * slope + abs((y-10)* side_slope)
63    return z
64
65domain.set_quantity('elevation', topography)  # Use function for elevation
66domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
67domain.set_quantity('stage',
68                    expression='elevation')   # zero depth initial condition
69
70
71
72#------------------------------------------------------------------------------
73# Compute normal depth flow characteristics for vee using Mannings equation
74#------------------------------------------------------------------------------
75
76normal_depth  = ( ref_flow*mannings_n* ( (2.0*( (1.0+(1.0/side_slope)**2.0)**0.5) )**0.6667 ) / (side_slope*(slope**0.5))  )**0.375
77area          = (normal_depth**2.0) *side_slope
78wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5)
79hyd_radius    = area/wet_perimeter
80flow_velocity = ref_flow/area
81top_width     = 2* side_slope*normal_depth
82
83print 'Mannings solution for normal depth flow in vee is; '
84print '    Channel sideslopes = ', side_slope
85print '    Channel n          = ', mannings_n
86print '    Channel slope      = ', slope
87print '    Reference flow     = ', ref_flow
88print '    Normal depth       = ', normal_depth
89print '    Flow area          = ', area
90print '    Wet perimeter      = ', wet_perimeter
91print '    Hydraulic radius   = ', hyd_radius
92print '    Avg velocity       = ', flow_velocity
93print '    Flow width         = ', top_width
94
95
96#------------------------------------------------------------------------------
97# Calculate xmomentum and depth to apply to the inflow boundary
98#------------------------------------------------------------------------------
99
100# NOTE: momentum is applied 'uniformly' to the 'full' inflow boundary by ANUGA
101# this is all very unrealistic as only part of the bdry is 'wet'
102# but should reasonably replicate ref_flow in ANUGA nevertheless
103inflow_xmomentum = ref_flow/width 
104inflow_depth     = normal_depth
105
106#------------------------------------------------------------------------------
107# Setup boundary conditions
108#------------------------------------------------------------------------------
109
110
111Br = Reflective_boundary(domain)       # Solid reflective wall on edge of sides
112Bi = Dirichlet_boundary([inflow_depth,inflow_xmomentum,0.0])    # inflow at top
113# set outflow depth to mannings normal depth for flow to avoid drawdown impact upstream
114def outflow_depth(t):
115    return (-slope*length) + normal_depth
116Bo = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=outflow_depth)
117
118
119domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
120
121#------------------------------------------------------------------------------
122# Evolve system through time
123#------------------------------------------------------------------------------
124
125for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
126    if verbose :
127        print domain.timestepping_statistics()
128
129
130#------------------------------------------------------------------------------
131# Compute flow across flowlines -v- ref_flow
132#------------------------------------------------------------------------------
133
134# Check flow across flowline at 2m (immediately DS of inflow bdry)
135q=domain.get_flow_through_cross_section([[2.0,0.0],[2.0,width]])
136print '90 degree flowline just ds of inflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
137
138# Check flow across flowline at 500m (where  it should have normalised)
139q=domain.get_flow_through_cross_section([[500.0,0.0],[500.0,width]])
140print '90 degree flowline half way down channel: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
141
142# Check flow across flowline at 998m (wimmediately US of outflow bdry)
143q=domain.get_flow_through_cross_section([[998.0,0.0],[998.0,width]])
144print '90 degree flowline just US of outflow bdry: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow)
145
146#------------------------------------------------------------------------------
147# Record temporal flow characteristics just ds inflow bdry
148#------------------------------------------------------------------------------
149
150# Store temporal pattern of gauges at head/middle/bottom of channel in csv files
151sww2csv_gauges(sww_file  = 'test_flowhydraulics_in_vee_using_momentum.sww',
152               gauge_file= 'gauge_locations.csv',
153               out_name  = 'momentum_',
154               quantities=['stage','speed','elevation'],
155               verbose=True)
Note: See TracBrowser for help on using the repository browser.