Changeset 5340


Ignore:
Timestamp:
May 16, 2008, 5:12:37 PM (16 years ago)
Author:
ole
Message:

Work on momentum jet

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/Rudy_vandrie_2007/culvert_development.py

    r5314 r5340  
    286286        openings = self.openings
    287287        for i, opening in enumerate(openings):
    288             stage = domain.quantities['stage'].get_values(location='centroids',
     288            # Quantities corresponding to fluid exchange field for this opening           
     289            stage_o = domain.quantities['stage'].get_values(location='centroids',
    289290                                                          indices=opening.indices)
    290             elevation = domain.quantities['elevation'].get_values(location='centroids',
     291            elevation_o = domain.quantities['elevation'].get_values(location='centroids',
    291292                                                                  indices=opening.indices)
     293            xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
     294                                                                  indices=opening.indices)
     295            ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids',
     296                                                                  indices=opening.indices)
     297
     298            # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account)
     299            depth_o = stage_o - elevation_o
     300           
     301            ux_o = xmomentum_o/(depth_o+epsilon)
     302            uy_o = ymomentum_o/(depth_o+epsilon)
     303            v_o = mean(sqrt(ux_o**2+uy_o**2))
     304            d_o = mean(depth_o)
     305            w_o = mean(stage_o)
     306            z_o = mean(elevation_o)
     307            self.mean_xmomentum_o = mean(xmomentum_o)
     308            self.mean_ymomentum_o = mean(ymomentum_o)             
    292309
    293310            # Indices corresponding to energy enquiry field for this opening
     
    317334            w = mean(stage)
    318335            z = mean(elevation)
    319            
     336            self.mean_xmomentum = mean(xmomentum)
     337            self.mean_ymomentum = mean(ymomentum)             
     338
    320339            # Ratio of depth to culvert height
    321340            ratio = d/(2*self.height)
     
    436455                delta_t = time - self.last_time
    437456                if delta_t > 0.0:
    438                     xmomentum_rate = outlet_mom_x - self.xmom_forcing1.value
    439                     xmomentum_rate /= delta_t
     457                    xmomentum_rate = outlet_mom_x - self.mean_xmomentum
     458                    if xmomentum_rate > 0:
     459                        xmomentum_rate /= delta_t
     460                    else:
     461                        xmomentum_rate = 0.0
    440462                   
    441                     ymomentum_rate = outlet_mom_y - self.ymom_forcing1.value
    442                     ymomentum_rate /= delta_t                       
     463                    ymomentum_rate = outlet_mom_y - self.mean_ymomentum
     464                    if ymomentum_rate > 0:
     465                        ymomentum_rate /= delta_t
     466                    else:
     467                        ymomentum_rate = 0.0                   
    443468                else:
    444469                    xmomentum_rate = ymomentum_rate = 0.0
     
    453478
    454479
    455                 if int(domain.time*100) % 100 == 0:
    456                     s = 'T=%.2f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
     480                if int(time*100) % 100 == 0:
     481                    s = 'T=%.3f, Culvert Discharge = %.3f  Culv. Vel. %0.3f'\
    457482                        %(time, flow_rate_control, outlet_vel)
    458                     s += ' Depth= %0.3f  Momentum = (%0.3f, %0.3f)'\
    459                          %(outlet_dep, outlet_mom_x,outlet_mom_y)
    460                     s += ' Momentum rate: (%.4f, %.4f)'\
     483                    s += ' Depth= %0.3f\n    Outlet Momentum = (%0.3f, %0.3f)\n'\
     484                         %(outlet_dep, outlet_mom_x, outlet_mom_y)
     485                    s += '    Avg Momentum at opening = (%0.3f, %0.3f)\n'\
     486                         %(self.mean_xmomentum_o, self.mean_ymomentum_o)
     487                    s += '    Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\
     488                         %(self.mean_xmomentum, self.mean_ymomentum)                   
     489                    s += '    Momentum rate: (%.4f, %.4f)'\
    461490                         %(xmomentum_rate, ymomentum_rate)                   
    462491                    fid.write(s)
     
    488517        # Print out flows every 1 seconds
    489518        if int(time*100) % 100 == 0:
    490             s = 'Time %.2f\n' %time
     519            s = 'Time %.3f\n' %time
    491520            s += '    Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\
    492521                 %(openings[0].depth,
     
    506535           
    507536            print s
     537
    508538            fid.write(s)
    509539
Note: See TracChangeset for help on using the changeset viewer.