Ignore:
Timestamp:
Nov 8, 2004, 6:32:02 PM (20 years ago)
Author:
ole
Message:

First play with wind stress

File:
1 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r462 r501  
    804804
    805805    w = domain.quantities['level'].centroid_values
     806    z = domain.quantities['elevation'].centroid_values
     807    h = w-z
     808   
    806809    uh = domain.quantities['xmomentum'].centroid_values
    807810    vh = domain.quantities['ymomentum'].centroid_values
     
    816819   
    817820    for k in range(N):
    818         if w[k] >= eps:
     821        if h[k] >= eps:
    819822            S = -g * eta[k]**2 * sqrt((uh[k]**2 + vh[k]**2))
    820             S /= w[k]**(7.0/3)
     823            S /= h[k]**(7.0/3)
    821824
    822825            #Update momentum
     
    834837       
    835838    w = domain.quantities['level'].centroid_values
     839    z = domain.quantities['elevation'].centroid_values
     840    h = w-z
    836841    uh = xmom.centroid_values
    837842    vh = ymom.centroid_values
     
    846851
    847852    from shallow_water_ext import manning_friction
    848     manning_friction(g, eps, w, uh, vh, eta, xmom_update, ymom_update)
    849        
     853    manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
     854       
     855
     856
     857       
     858def wind_stress(domain):
     859    """Apply wind stress to water momentum
     860    """
     861
     862    #FIXME: Under construction
     863
     864    #Hardwired velocity field - should be initialised as a
     865    #callable object
     866
     867    u = 1.0
     868    v = 0.0
     869
     870    Cw = 3.0e-3 #Wind stress coeffficient
     871    rho_a = 1.2e-3 #Atmospheric density
     872    rho = 1023 #Density of water
     873   
     874    from math import sqrt
     875
     876    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
     877    ymom_update = domain.quantities['ymomentum'].semi_implicit_update   
     878   
     879    N = domain.number_of_elements
     880
     881    c = Cw*rho_a/rho
     882    S = c * sqrt(u**2 + v**2)   
     883    for k in range(N):
     884        #Update momentum
     885        xmom_update[k] += S*u
     886        ymom_update[k] += S*v
     887           
     888       
     889
    850890
    851891###########################
Note: See TracChangeset for help on using the changeset viewer.