Acoustic Levitation

#!/usr/bin/python3

# Acoustic wave propogation and reflection model used for
# predicting pressure inflection points for acoustic levitationp.
# Determines pressure distribution for a given 2D area assuming a
# circular transducer radiation surface and circular flat reflector.


import pdb
import numpy as np
import matplotlib.pylab as plt
import scipy as sci
######### MODEL PARAMETERS #####################

mm = 1e-3

# Grid sizing (mm) - this can be changed
x_min = -10*mm
x_max = 10*mm
z_min = 0
z_max = 15*mm

disc = 0.5*mm   # Discretization size - this defines the resolution of the data
R = 5*mm        # Transducer radius (mm)

############################################

# Establish THE GRID
x = np.linspace(x_min, x_max, int((x_max-x_min)/disc))
z = np.linspace(z_min, z_max, int((z_max-z_min)/disc))
transducer = np.linspace(-R, R, int((2*R)/disc))


# Array indices so I can keep track of dimensions
N = int((2*R)/disc)
M = int(len(x)*len(z))
L = int(len(x))
Z = int(len(z))

# Distance arrays
r_nm = np.zeros((N, M))
r_im = np.zeros((L, M))
r_in = np.zeros((N, L))

# print(the array sizes to get a feel for how long this shit will take to run
print('N: {}'.format(N))
print('M: {}'.format(M))
print('L: {}'.format(L))
print('Z: {}'.format(Z))

# DISTANCE ARRAYS ####################################

# Calculate r_nm
for i in range(N):
    q = 0
    for j in range(L):
        for k in range(len(z)):
            r_nm[i, q] = np.sqrt((transducer[i]-x[j])**2+(z_max-z[k])**2)
            q = q+1


# Calculate r_im
for i in range(L):
    q = 0

    for j in range(L):
        for k in range(len(z)):
            r_im[i, q] = np.sqrt((x[i]-x[j])**2+(z_min-z[k])**2)
            q = q + 1


# Calculate r_in
for i in range(N):
    for j in range(L):
        r_in[i, j] = np.sqrt((transducer[i]-x[j])**2+(z_max-z_min)**2)
r_ni = r_in.T

##################################################################

# Constants and stuff
# These need to be complex so Python doesn't shit the bed when it sees a complex calculation

A_trans = complex(np.pi*R**2, 0)        # Transducer area, m^2
A_refl = complex(np.pi*(x_max/2)**2, 0)  # Reflector area, m^2
rho = complex(1.214, 0)         # Density of air in Raleigh, kg/m^3
c = complex(340.1, 0)           # Speed of sound in air, m/s
f = complex(28000, 0)           # Frequency, hZ
wavelength = complex(c/f, 0)        # Wavelength, m
omega = complex(2*np.pi*f, 0)       # Angular frequency
cells = complex(100, 0)         # Number of discrete cells being used for calculation
sn = A_trans/cells          # Unit cell area for transducer
# Unit cell area for reflector - factor of 4 to keep area/cell same for both
si = A_refl/(cells*4)
D = (omega*rho*c)/wavelength        # Constant used for primary wave
wavenumber = omega/c            # .....it's the wavenumber - k
# Amplitude (m), estimated from http://www.mtiinstruments.com/mtiuniversity/appnotes/mtiuniversity-appnotes-ultrasonic.aspx
U_0 = 0.0000060
E = complex(0, 1)/wavelength        # Constant used for reflected waves


################# TRANSFER MATRICIES ##################################################

T_TM = np.matrix(np.zeros((N, M)), dtype=complex)
T_TR = np.matrix(np.zeros((N, L)), dtype=complex)
T_RT = np.matrix(np.zeros((L, N)), dtype=complex)
T_RM = np.matrix(np.zeros((L, M)), dtype=complex)


for i in range(N):
    for j in range(M):
        T_TM[i, j] = (sn*np.exp(complex(0, -1) *
                                wavenumber*r_nm[i, j]))/r_nm[i, j]

    for k in range(L):
        T_TR[i, k] = (sn*np.exp(complex(0, -1) *
                                wavenumber*r_in[i, k]))/r_in[i, k]

for i in range(L):
    for j in range(N):
        T_RT[i, j] = (si*np.exp(complex(0, -1) *
                                wavenumber*r_ni[i, j]))/r_ni[i, j]

for i in range(L):
    # Feedback to keep track of where we are in the calculations, this one can take awhile

    for j in range(M):
        if r_im[i, j] == 0.0:
            continue
        T_RM[i, j] = (si*np.exp(complex(0, -1) * wavenumber*r_im[i, j]))/r_im[i, j]
        print('{} {} {}'.format(i,j,T_RM[i,j]))
print('si: {}, wavenumber: {}, r_im[i,j]: {}'.format(
    si, wavenumber, r_im[i, j]))

print('si * exp(0,-1): {}'.format(si*np.exp(complex(0, -1))))
print('wavenumber {} '.format(wavenumber))
print('si * exp(0,-1) * wavenumber {}'.format(si*np.exp(complex(0, -1) * wavenumber)))


# Rotate the matricies so dimensions work out for arithmetic
T_TM = T_TM.T
T_TR = T_TR.T
T_RT = T_RT.T
T_RM = T_RM.T

######################## DISPLACEMENT ARRAY #######################################

# U will be the same for every point N on the transducer
U = np.array(np.zeros((N, 1)), dtype=complex)
for i in range(N):
    U[i] = U_0*np.exp(complex(0, omega))

######################### PRESSURE CALCULATIONS #################################

# Initial wave propogation + 4 reflected waves - higher order terms can be added or removed for accuracy or efficiency
P = (D*E)*T_RM*(T_TR*U)+(D*E**2)*(T_TM*T_RT)*(T_TR*U)+(D*E**3)*T_RM * \
    (T_TR*T_RT)*(T_TR*U)+(D*E**4)*(T_TM*T_RT)*(T_TR*T_RT)*(T_TR*U)

# Reshape the pressure array so Z is up and X is across
P2 = np.array(np.reshape(P, (len(x), len(z))))
P2 = np.transpose(P2)
###################################################################

rho = 1.214
c = 340.1
lamb = 0.0122
omega = 28000
k = omega/c
rad = 2 * np.pi / 180
cell = 100
sn = np.pi*rad**2/cell
u = 0.000006


p = np.zeros(20)
a = np.zeros(20)
for j in range(20):
    p[j] = omega*rho*c*(1/lamb)*sn*(1/(lamb*(j+1)/2) *
                                    u*np.exp(-i*k*(lamb*(j+1)/2)))
    a[j] = rho*c*(1/lamb)*sn*(1/(lamb*(j+1)/2))*u*np.exp(-i*k*(lamb*(j+1)/2))+omega*rho*c * \
        (1/lamb)*sn*(1/(lamb*(j+1)/2))*u*np.exp(-i*k *
                                                (lamb*(j+1)/2))*(-i*(lamb*(j+1)/2)*(1/c))*1000

# print(p)
# print(a)


plt.pcolormesh(x, z, P2.real)
plt.title('Pressure Distribution')
plt.xlabel('X (m)')
plt.ylabel('Z (m)')
# Scale the pressure gradients so the colors make sense
plt.clim(min(P.real)/2, -min(P.real)/2.5)
plt.colorbar().set_label('Relative Pressure')
# savefig('odd_Wavelength.png',dpi=300) # Save the image with given parameters
plt.show()