Changeset 7403
 Timestamp:
 Aug 24, 2009, 2:17:07 PM (13 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/ted_rigby/test_flow_in_vee_basic.py
r7402 r7403 36 36 width = 20. 37 37 side_slope = 1.0 38 dx = dy = 1# Resolution: of points (elev) grid on both axes38 dx = dy = 2 # Resolution: of points (elev) grid on both axes 39 39 ref_flow = 5.0 # Inflow at head of vee 40 40 slope = 1.0/1000 … … 48 48 domain = Domain(points, vertices, boundary) 49 49 domain.set_name('test_flow_in_vee') # Output name 50 51 52 # 53 # Setup initial conditions 54 # 55 56 def topography(x, y): 57 z=x * slope + abs((y10)* side_slope) 58 return z 59 60 domain.set_quantity('elevation', topography) # Use function for elevation 61 domain.set_quantity('friction', mannings_n) # Constant friction of conc surface 62 domain.set_quantity('stage', 63 expression='elevation') # Dry initial condition 50 print domain.statistics() 64 51 65 52 … … 72 59 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 73 60 area = (normal_depth**2.0) *side_slope 74 wet_perimeter = 2* 61 wet_perimeter = 2*normal_depth*((1+(side_slope)**2.0)**0.5) 75 62 hyd_radius = area/wet_perimeter 76 63 flow_velocity = ref_flow/area 77 top_width = 2* 64 top_width = 2*side_slope*normal_depth 78 65 79 66 # Compute xmomentum to input into the ANUGA inflow boundary 80 67 # assuming momentum is applied uniformly to inflow boundary by ANUGA 81 68 82 xmomentum = ref_flow/top_width # this is all very unrealistic but should replicate ref_flow in ANUGA 69 #xmomentum = ref_flow/top_width # this is all very unrealistic but should replicate ref_flow in ANUGA 70 xmomentum = ref_flow/width # This only applies where channel is flat. 71 72 print 'ref_flow', ref_flow 73 print 'xmomentum', xmomentum 74 print 'normal_depth', normal_depth 75 print 'width', width 76 77 78 79 # 80 # Setup initial conditions 81 # 82 83 def topography(x, y): 84 z = x * slope + abs((y10)*side_slope) 85 86 # Let both ends of domain be square so that the applied xmomentum 87 # represents flow correctly 88 for i in range(len(x)): 89 if x[i] < 20 or x[i] > 380: 90 z[i] = x[i] * slope 91 92 return z 93 94 95 def initial_stage(x, y): 96 """Set constant depth based on calculated normal flow 97 """ 98 99 w = x*slope + normal_depth 100 return w 101 102 103 domain.set_quantity('elevation', topography) # Use function for elevation 104 domain.set_quantity('friction', mannings_n) # Constant friction of conc surface 105 domain.set_quantity('stage', initial_stage) 106 domain.set_quantity('xmomentum', xmomentum) # Initial momentum commensurate with normal flow 107 108 83 109 84 110 # … … 89 115 Br = Reflective_boundary(domain) # Solid reflective wall on edge of sides 90 116 Bi = Dirichlet_boundary([normal_depth,xmomentum,0.0]) # inflow at top 91 Bo = Transmissive_boundary(domain) # momentum as above bdry 117 Bo = Dirichlet_boundary([length*slope+normal_depth,xmomentum,0.0]) # Outflow at bottom 118 92 119 domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br}) 93 120 … … 98 125 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime): 99 126 if verbose : 100 print domain.timestepping_statistics() 127 print domain.timestepping_statistics() 128 129 130 # Check flow across flowline at 200m at every yield step 131 q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,width]]) 132 print '90 degree flowline: ANUGA Q = %f, Ref_flow = %f' %(q, ref_flow) 101 133 102 134
Note: See TracChangeset
for help on using the changeset viewer.