Ignore:
Timestamp:
May 14, 2011, 10:17:08 PM (14 years ago)
Author:
paul
Message:

First implementation of pressurised flow

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/2010-projects/anuga_1d/pipe/pipe_domain.py

    r8153 r8177  
    44This module contains a specialisation of class Generic_domain from module generic_domain.py consisting of methods specific to mixed rectangular Pipe/Channel flow using the Shallow Water Wave Equation.
    55
    6 This particular modification of the Domain class implements the ability to
    7 vary the width of the 1D pipe that the water flows in. As a result the
    8 conserved variables are different than previous implementations and so are the
    9 equations. In particular, the conserved quantites generalise the open channel flow conserved quantities.
     6This particular modification of the Domain class implements the ability to vary the height and width of the 1D pipe that the water flows in. As a result the conserved variables are slightly different than previous implementations and so are the equations. In particular, the conserved quantites generalise the open channel flow conserved quantities.
    107
    118U_t + E_x = S
    … …  
    1411
    1512U = [A, Q] = [b*h, u*b*h]
    16 E = [Q, Q^2/A + g*b*h^2/2]
     13E = [Q, Q^2/A + A]
    1714S represents source terms forcing the system
    18 (e.g. gravity, boundary_stree, friction, wind stress, ...)
     15(e.g. gravity, boundary_stress, friction, wind stress, ...)
    1916gravity = -g*b*h*z_x
    2017boundary_stress = 1/2*g*b_x*h^2
    2118
    22 and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
     19and _t, _x, denote the derivative with respect to t and x respectiely.
    2320
    2421The quantities are
    2522
    2623symbol    variable name    explanation
    27 A         area             equivalent Wetted area = b*h
    28 Q         discharge        flux of water = u*b*h
     24A         area             free surface equivalent wetted area = b*h [m^2]
     25Q         discharge        flux of water = u*b*h [m^3/s]
    2926x         x                horizontal distance from origin [m]
    3027z         elevation        elevation of bottom of pipe [m]
    31 h         height           water height above z [m]
     28h         height           free surface equivalent water height above z [m]
    3229w         stage            equivalent absolute water level, w = z+h [m]
    33 u                          speed in the x direction [m/s]
     30u         xvelocity        speed in the x direction [m/s]
    3431uh        xmomentum        momentum in the x direction [m^2/s]
    35 b         width            width of pipe
     32b         width            width of pipe [m]
     33t         top              height of pipe above z [m]
    3634eta                        mannings friction coefficient [to appear]
    3735
    … …  
    6058
    6159        conserved_quantities = ['area', 'discharge']
    62         evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','stage']
     60        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
    6361        other_quantities = ['friction']
    6462        Generic_domain.__init__(self,
    … …  
    124122        msg = 'Sixth evolved quantity must be "width"'
    125123        assert self.evolved_quantities[5] == 'width', msg
    126         msg = 'Seventh evolved quantity must be "stage"'
    127         assert self.evolved_quantities[6] == 'stage', msg
     124        msg = 'Seventh evolved quantity must be "top"'
     125        assert self.evolved_quantities[6] == 'top', msg
     126        msg = 'Eighth evolved quantity must be "stage"'
     127        assert self.evolved_quantities[7] == 'stage', msg
    128128
    129129        Generic_domain.check_integrity(self)
    … …  
    155155    velocity   = domain.quantities['velocity']
    156156    width      = domain.quantities['width']
     157    top        = domain.quantities['top']
    157158
    158159
    159160    from anuga_1d.pipe.pipe_domain_ext import compute_fluxes_pipe_ext
    160     domain.flux_timestep = compute_fluxes_pipe_ext(timestep,domain,area,discharge,bed,height,velocity,width)
     161    domain.flux_timestep = compute_fluxes_pipe_ext(timestep,domain,area,discharge,bed,height,velocity,width,top)
    161162
    162163#-----------------------------------------------------------------------
    … …  
    183184    velocity  = domain.quantities['velocity']
    184185    width     = domain.quantities['width']
     186    top       = domain.quantities['top']
    185187    stage     = domain.quantities['stage']
    186188
    … …  
    192194    u_C   = velocity.centroid_values
    193195    b_C   = width.centroid_values
     196    t_C   = top.centroid_values
    194197    w_C   = stage.centroid_values
    195198
    … …  
    198201        domain.setstageflag = False
    199202
     203    # Paul: Same for top?
    200204    if domain.discontinousb:
    201205        width.extrapolate_second_order()
    … …  
    204208    h0 = 1.0e-12
    205209
     210    # Paul: How does top (i.e. t_C) affect these?
    206211    h_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C/(b_C+h0/b_C), 0.0 )
    207212    u_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C/(a_C+h0/a_C), 0.0 )
    … …  
    209214    a_C[:] = numpy.where( (a_C>h0) | (b_C>h0), a_C, 0.0 )
    210215    b_C[:] = numpy.where( (a_C>h0) | (b_C>h0), b_C, 0.0 )
     216    t_C[:] = numpy.where( (a_C>h0) | (b_C>h0), t_C, 0.0 )
    211217    d_C[:] = numpy.where( (a_C>h0) | (b_C>h0), d_C, 0.0 )
    212218
    … …  
    214220   
    215221               
    216 
     222    # Paul: same for top?
    217223    bed.extrapolate_second_order()
    218224   
    … …  
    234240    d_V  = discharge.vertex_values
    235241    b_V  = width.vertex_values
     242    t_V  = top.vertex_values
    236243
    237244
    … …  
    274281        self.velocity = domain.quantities['velocity'].vertex_values
    275282        self.width    = domain.quantities['width'].vertex_values
     283        self.top    = domain.quantities['top'].vertex_values
    276284        self.stage    = domain.quantities['stage'].vertex_values
    277285
    278         self.evolved_quantities = numpy.zeros(7, numpy.float)
     286        self.evolved_quantities = numpy.zeros(8, numpy.float)
    279287
    280288    def __repr__(self):
    … …  
    294302        q[4] = -self.velocity[vol_id, edge_id]
    295303        q[5] =  self.width[vol_id,edge_id]
    296         q[6] =  self.stage[vol_id,edge_id]
     304        q[6] =  self.top[vol_id,edge_id]
     305        q[7] =  self.stage[vol_id,edge_id]
    297306
    298307        return q
    … …  
    313322            raise msg
    314323
    315         assert len(evolved_quantities) == 7
     324        assert len(evolved_quantities) == 8
    316325
    317326        self.evolved_quantities=numpy.array(evolved_quantities,numpy.float)
    … …  
    341350    Height    = domain.quantities['height']
    342351    Width     = domain.quantities['width']
     352    Top       = domain.quantities['top']
    343353
    344354    discharge_ud  = Discharge.explicit_update
    … …  
    348358    h = Height.vertex_values
    349359    b = Width.vertex_values
     360    t = Top.vertex_values
    350361    a = Area.vertex_values
    351362    z = Elevation.vertex_values
    … …  
    356367        avg_h = 0.5*(h[k,0] + h[k,1])
    357368        avg_b = 0.5*(b[k,0] + b[k,1])
     369        avg_t = 0.5*(t[k,0] + t[k,1])
    358370       
    359371        #Compute bed slope
    … …  
    376388    Height    = domain.quantities['height']
    377389    Width     = domain.quantities['width']
     390    Top       = domain.quantities['top']
    378391
    379392    discharge_ud  = Discharge.explicit_update
    … …  
    381394    h = Height.vertex_values
    382395    b = Width.vertex_values
     396    t = Top.vertex_values
    383397    a = Area.vertex_values
    384398    z = Elevation.vertex_values
Note: See TracChangeset for help on using the changeset viewer.