Changeset 7403


Ignore:
Timestamp:
Aug 24, 2009, 2:17:07 PM (11 years ago)
Author:
ole
Message:

Played with Ted's V channel

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/ted_rigby/test_flow_in_vee_basic.py

    r7402 r7403  
    3636width      = 20.
    3737side_slope = 1.0
    38 dx = dy    = 1          # Resolution: of points (elev) grid on both axes
     38dx = dy    = 2          # Resolution: of points (elev) grid on both axes
    3939ref_flow   = 5.0        # Inflow at head of vee
    4040slope      = 1.0/1000
     
    4848domain = Domain(points, vertices, boundary)   
    4949domain.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((y-10)* 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
     50print domain.statistics()
    6451
    6552
     
    7259normal_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
    7360area          = (normal_depth**2.0) *side_slope
    74 wet_perimeter = 2* normal_depth*((1+(side_slope)**2.0)**0.5)
     61wet_perimeter = 2*normal_depth*((1+(side_slope)**2.0)**0.5)
    7562hyd_radius    = area/wet_perimeter
    7663flow_velocity = ref_flow/area
    77 top_width     = 2* side_slope*normal_depth
     64top_width     = 2*side_slope*normal_depth
    7865
    7966# Compute xmomentum to input into the ANUGA inflow boundary
    8067# assuming momentum is applied uniformly to inflow boundary by ANUGA
    8168
    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
     70xmomentum = ref_flow/width       # This only applies where channel is flat.
     71
     72print 'ref_flow', ref_flow
     73print 'xmomentum', xmomentum
     74print 'normal_depth', normal_depth
     75print 'width', width
     76
     77
     78
     79#------------------------------------------------------------------------------
     80# Setup initial conditions
     81#------------------------------------------------------------------------------
     82
     83def topography(x, y):
     84    z = -x * slope + abs((y-10)*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   
     95def 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   
     103domain.set_quantity('elevation', topography)  # Use function for elevation
     104domain.set_quantity('friction', mannings_n)   # Constant friction of conc surface   
     105domain.set_quantity('stage', initial_stage)
     106domain.set_quantity('xmomentum', xmomentum) # Initial momentum commensurate with normal flow
     107
     108
    83109
    84110#------------------------------------------------------------------------------
     
    89115Br = Reflective_boundary(domain)                       # Solid reflective wall on edge of sides
    90116Bi = Dirichlet_boundary([normal_depth,xmomentum,0.0])  # inflow at top
    91 Bo = Transmissive_boundary(domain)                     # momentum as above bdry
     117Bo = Dirichlet_boundary([-length*slope+normal_depth,xmomentum,0.0])  # Outflow at bottom
     118
    92119domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})
    93120
     
    98125for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
    99126    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)     
    101133
    102134
Note: See TracChangeset for help on using the changeset viewer.