Calculate Divergence of Velocity Field (3D) in Python

1.2k views Asked by At

I am trying to calculate the divergence of a 3D velocity field in a multi-phase flow setting (with solids immersed in a fluid). If we assume u,v,w to be the three velocity components (each a n x n x n) 3D numpy array, here is the function I have for calculating divergence:

def calc_divergence_velocity(df,h=0.025):
    """        
        @param df: A dataframe with the entire vector field with columns [x,y,z,u,v,w] with 
                   x,y,z indicating the 3D coordinates of each point in the field and u,v,w
                   the velocities in the x,y,z directions respectively.  
        @param h: This is the dimension of a single side of the 3D (uniform) grid. Used 
                  as input to numpy.gradient() function.
    """

    """
        Reshape dataframe columns to get 3D numpy arrays (dim = 80) so each u,v,w is a 
        80x80x80 ndarray.
    """ 
   
    u = df['u'].values.reshape((dim,dim,dim))
    v = df['v'].values.reshape((dim,dim,dim))
    w = df['w'].values.reshape((dim,dim,dim))
    
    #Supply x,y,z coordinates appropriately.
    
    #Note: Only a scalar `h` has been supplied to np.gradient because
    #the type of grid we are dealing with is a uniform grid with each
    #grid cell having the same dimensions in x,y,z directions.

    u_grad = np.gradient(u,h,axis=0)   #central diff. du_dx
    v_grad = np.gradient(v,h,axis=1)   #central diff. dv_dy
    w_grad = np.gradient(w,h,axis=2)   #central diff. dw_dz
    
    """
        The `mask` column in the dataframe is a binary column indicating the locations
        in the field where we are interested in measuring divergence.
        The problem I am looking at is multi-phase flow with solid particles and a fluid
        hence we are only interested in the fluid locations.
    """
    sdf = df['mask'].values.reshape((dim,dim,dim))  
    

    div = (u_grad*sdf) + (v_grad*sdf) + (w_grad*sdf)

    return div

The problem I'm having is that the divergence values that I am seeing are far too high. For example the image below showcases, a distribution with values between [-350,350] whereas most values should technically be close to zero and somewhere between [20,-20] in my case. This tells me I'm calculating the divergence incorrectly and I would like some pointers as to how to correct the above function to calculate the divergence appropriately. As far as I can tell (please correct me if I'm wrong), I think have done something similar to this upvoted SO response. Thanks in advance!

enter image description here

0

There are 0 answers