1 | """test_flow_in_vee_using_momentum.py |
---|
2 | |
---|
3 | Test the ability of ANUGA to |
---|
4 | a) maintain the flow rate input into the channel at the upstream end |
---|
5 | b) replicate normal depth half way down the channel after half an hour |
---|
6 | c) replicate a typical velocity profile across the channel after half an hour |
---|
7 | |
---|
8 | The vee channel is 20m wide by 1000m long with side slopes of 1:1 |
---|
9 | The equations for dn are generalised with side slope as 1 V in Side_Slope H |
---|
10 | The test location is in the invert of (y=10m)and half way down the channel (x=500m) |
---|
11 | This script explores flow behaviour predicted by ANUGA and compares it with known solutions |
---|
12 | |
---|
13 | This test inputs flow using momentum=flow/width . |
---|
14 | |
---|
15 | """ |
---|
16 | |
---|
17 | verbose = True |
---|
18 | import numpy as num |
---|
19 | |
---|
20 | #------------------------------------------------------------------------------ |
---|
21 | # Import necessary modules |
---|
22 | #------------------------------------------------------------------------------ |
---|
23 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
24 | from anuga.shallow_water import Domain |
---|
25 | from anuga.shallow_water.shallow_water_domain import Reflective_boundary |
---|
26 | from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary |
---|
27 | from anuga.shallow_water.shallow_water_domain import Transmissive_Momentum_Set_Stage_boundary |
---|
28 | from anuga.shallow_water.shallow_water_domain import Transmissive_boundary |
---|
29 | from anuga.shallow_water.shallow_water_domain import Inflow |
---|
30 | from anuga.shallow_water.data_manager import get_flow_through_cross_section |
---|
31 | from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs |
---|
32 | |
---|
33 | |
---|
34 | #------------------------------------------------------------------------------ |
---|
35 | # Setup computational domain |
---|
36 | #------------------------------------------------------------------------------ |
---|
37 | |
---|
38 | finaltime = 1800.0 # 1/2 hour |
---|
39 | |
---|
40 | length = 1000. |
---|
41 | width = 20. |
---|
42 | side_slope = 1.0 |
---|
43 | dx = dy = 1 # Resolution: of points (elev) grid on both axes |
---|
44 | ref_flow = 20.0 # Inflow at head of vee |
---|
45 | slope = 1.0/1000 |
---|
46 | mannings_n = 0.150 |
---|
47 | |
---|
48 | points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy), |
---|
49 | len1=length, len2=width) |
---|
50 | |
---|
51 | |
---|
52 | |
---|
53 | domain = Domain(points, vertices, boundary) |
---|
54 | domain.set_name('test_flowhydraulics_in_vee_using_momentum') # Output name |
---|
55 | |
---|
56 | |
---|
57 | #------------------------------------------------------------------------------ |
---|
58 | # Setup initial conditions |
---|
59 | #------------------------------------------------------------------------------ |
---|
60 | |
---|
61 | def topography(x, y): |
---|
62 | z=-x * slope + abs((y-10)* side_slope) |
---|
63 | return z |
---|
64 | |
---|
65 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
66 | domain.set_quantity('friction', mannings_n) # Constant friction of conc surface |
---|
67 | domain.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 | |
---|
76 | normal_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 |
---|
77 | area = (normal_depth**2.0) *side_slope |
---|
78 | wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5) |
---|
79 | hyd_radius = area/wet_perimeter |
---|
80 | flow_velocity = ref_flow/area |
---|
81 | top_width = 2* side_slope*normal_depth |
---|
82 | |
---|
83 | print 'Mannings solution for normal depth flow in vee is; ' |
---|
84 | print ' Channel sideslopes = ', side_slope |
---|
85 | print ' Channel n = ', mannings_n |
---|
86 | print ' Channel slope = ', slope |
---|
87 | print ' Reference flow = ', ref_flow |
---|
88 | print ' Normal depth = ', normal_depth |
---|
89 | print ' Flow area = ', area |
---|
90 | print ' Wet perimeter = ', wet_perimeter |
---|
91 | print ' Hydraulic radius = ', hyd_radius |
---|
92 | print ' Avg velocity = ', flow_velocity |
---|
93 | print ' 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 |
---|
103 | inflow_xmomentum = ref_flow/width |
---|
104 | inflow_depth = normal_depth |
---|
105 | |
---|
106 | #------------------------------------------------------------------------------ |
---|
107 | # Setup boundary conditions |
---|
108 | #------------------------------------------------------------------------------ |
---|
109 | |
---|
110 | |
---|
111 | Br = Reflective_boundary(domain) # Solid reflective wall on edge of sides |
---|
112 | Bi = 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 |
---|
114 | def outflow_depth(t): |
---|
115 | return (-slope*length) + normal_depth |
---|
116 | Bo = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=outflow_depth) |
---|
117 | |
---|
118 | |
---|
119 | domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) |
---|
120 | |
---|
121 | #------------------------------------------------------------------------------ |
---|
122 | # Evolve system through time |
---|
123 | #------------------------------------------------------------------------------ |
---|
124 | |
---|
125 | for 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) |
---|
135 | q=domain.get_flow_through_cross_section([[2.0,0.0],[2.0,width]]) |
---|
136 | print '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) |
---|
139 | q=domain.get_flow_through_cross_section([[500.0,0.0],[500.0,width]]) |
---|
140 | print '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) |
---|
143 | q=domain.get_flow_through_cross_section([[998.0,0.0],[998.0,width]]) |
---|
144 | print '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 |
---|
151 | sww2csv_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) |
---|