在这么短的时间内，这可能是一个长镜头，但我一直在努力为我正在从事的项目获取理论数据。该项目涉及使用非常简化的磁粉成像版本来分析流体。
因此，我在一些帮助下编写了一些代码，并且我很确定所有物理特性都是正确的，但是当我以 5E-5 的时间步长 (deltaT) 运行它时，代码会计算数千量级的速度，而它应该是以毫米每秒为单位。但是当代码的时间步长为 5E-6 时，它可以正确运行。问题是在这个时间步完成这些模拟所花费的时间很长，在我的实际实验中，其中一些花费了 20 多分钟。
基本上，我想知道是否有可能在这个时间步显着减少运行时间或能够增加模拟的时间步长。
这是我的代码：
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
#Arrays for storing time
#t1 = [0]
#t2 = [0]
#t3 = [0]
#t4 = [0]
#t5 = [0]
t6 = [0]
#Time Increment
deltaT = 5E-7
#Number of Particles
N = 100
#Viscosity
mu_w = 0.0010049
mu = 3*mu_w
#Geometry
radius = 0.0095
d = 0.023
L = 0.02
#Constants
mu_0 = 1.25663706E-6
#Particle Properties
p_rad = (4.9E-6)/2
m = 1.663223552E-13
moment = 60*m
x, y = np.mgrid[-1.15:1.15:151j, -1:1:151j]*1E-2
#Method for Magnetic Field
def B(x,y):
H11 = 4/((L+2*y)*np.sqrt((4*(d-2*x)**2/(L+2*y)**2+4)))
H12 = 4/((L+2*y)*np.sqrt((4*(d+2*x)**2/(L+2*y)**2+4)))
H13 = 4/((L-2*y)*np.sqrt((4*(d+2*x)**2/(L-2*y)**2+4)))
H14 = 4/((L-2*y)*np.sqrt((4*(d-2*x)**2/(L-2*y)**2+4)))
H21 = 1/((d-2*x)**2*(np.sqrt(4*(d-2*x)**2/(L-2*y)**2+4)))
H22 = 1/((d-2*x)**2*(np.sqrt(4*(d-2*x)**2/(L+2*y)**2+4)))
H23 = 1/((d+2*x)**2*(np.sqrt(4*(d+2*x)**2/(L-2*y)**2+4)))
H24 = 1/((d+2*x)**2*(np.sqrt(4*(d+2*x)**2/(L+2*y)**2+4)))
H1 = (H11 - H12 - H13 + H14)
H2 = (-4*(d-2*x))*(H21 - H22) + (-4*(d+2*x))*(H23 - H24)
B1 = (H1*mu_0)*2666
B2 = (H2*mu_0)*1700
B = np.sqrt(B1**2+B2**2)
return np.dstack((B1,B2,B))
B1 = B(x,y)[:,:,0]
B2 = B(x,y)[:,:,1]
Bt = B(x,y)[:,:,2]
#Contour Plot
levels = np.logspace(-2.5,1.23,10,base=10)
cp = plt.contourf(x, y, Bt, levels=levels, cmap=plt.cm.plasma,extend='max')
cb = plt.colorbar(cp)
plt.clim(0,2)
#X and Y increments
dxy = [0.023/151,0.02/151]
#Gradients of B, Bx & By
dB = np.gradient(Bt,dxy[0],dxy[1])
dBy = np.gradient(B1,dxy[1])[0]
dBx = np.gradient(B2,dxy[0])[1]
#Vector Field Plot
x, y = np.mgrid[-1.15:1.15:20j, -1:1:20j]*1E-2
Bq = B(x,y)[:,:,2]
dxy = [0.023/20, 0.02/20]
dBq = np.gradient(Bq,dxy[0],dxy[1])
quiv = plt.quiver(x, y, dBq[0], dBq[1])
plt.show()
#X Gradient Plot
plt.figure()
x = np.linspace(-0.0115,0.0115,151)
plt.plot(x,B2[75,:])
plt.title('Magnetic Field vs. Position in x-direction')
plt.xlabel('Position (m)')
plt.ylabel('Magnetic Field (T)')
plt.show()
#Y Gradient Plot
plt.figure()
y = np.linspace(-0.01,0.01,151)
plt.plot(y,-B1[:,75])
plt.title('Magnetic Field vs. Position in y-direction')
plt.xlabel('Position (m)')
plt.ylabel('Magnetic Field (T)')
plt.show()
#-----------------------------------------
yslice = 0.00375
xp = np.random.uniform(-radius,radius,N)
yp = np.random.uniform(-radius,radius,N)
v_x = np.zeros(N)
v_y = np.zeros(N)
#Arrays to store counts
#count1 = [0]
#count2 = [0]
#count3 = [0]
#count4 = [0]
#count5 = [0]
count6 = [0]
for i in range(0,N):
while xp[i]**2 + yp[i]**2 > radius**2:
xp[i] = np.random.uniform(-radius,radius,1)
yp[i] = np.random.uniform(-radius,radius,1)
for i in range(0,N):
if -yslice < yp[i] < yslice:
count6[0] = count6[0] + 1
fig1 = plt.figure()
plt.plot(xp,yp,".")
Fmag_x = np.zeros(N)
Fmag_y = np.zeros(N)
Fdrag_x = np.zeros(N)
Fdrag_y = np.zeros(N)
deltaV_x = np.zeros(N)
deltaV_y = np.zeros(N)
timestep = 1
#Method to animate simulation
def update(j,plt,fig1):
global xp, yp, v_x, v_y, t, timestep, count
plt.cla()
if min(xp**2 + yp**2) < radius**2:
count6.append(0)
t6.append(t6[timestep-1] + deltaT)
print(t6[timestep])
#print('fmag: ' + str(Fmag_x[0]))
#print('fdrag: ' + str(Fdrag_x[0]))
print('v(x): ' + str(v_x[0]))
for i in range(0,N):
if xp[i]**2 + yp[i]**2 < radius**2:
#Interpolation Method
Ix = int((xp[i]-(-0.0115))/(0.0115-(-0.0115))*150)
#print(Ix)
Iy = int((yp[i]-(-0.01))/(0.01-(-0.01))*150)
#print(Iy)
dBx1 = dBx[Ix,Iy]
dBx2 = dBx[Ix,Iy+1]
dBx3 = dBx[Ix+1,Iy]
dBx4 = dBx[Ix+1,Iy+1]
dBy1 = dBy[Ix,Iy]
dBy2 = dBy[Ix,Iy+1]
dBy3 = dBy[Ix+1,Iy]
dBy4 = dBy[Ix+1,Iy+1]
T = (xp[i]-x[Ix])/(x[Ix+1]-x[Ix])
#print(T)
U = (yp[i]-y[Iy])/(y[Iy+1]-y[Iy])
#print(U)
#Interpolated Gradients
dBxx = (1-T)*(1-U)*dBx1+(1-T)*U*dBx2+T*(1-U)*dBx3+T*U*dBx4
#print('dBx ' +str(dBxx))
dByy = (1-T)*(1-U)*dBy1+(1-T)*U*dBy2+T*(1-U)*dBy3+T*U*dBy4
#print('dBy ' +str(dByy))
#Magnetic Force in X and Y
Fmag_x[i] = moment*dBxx*(B(xp[i],yp[i])[:,:,0]/np.abs(B(xp[i],yp[i])[:,:,2]))
Fmag_y[i] = moment*dByy*(B(xp[i],yp[i])[:,:,1]/np.abs(B(xp[i],yp[i])[:,:,2]))
#Drag Force in X and Y
Fdrag_x[i] = 6*np.pi*mu*p_rad*v_x[i]
Fdrag_y[i] = 6*np.pi*mu*p_rad*v_y[i]
#Increase in velocity
deltaV_x[i] = (Fmag_x[i] - Fdrag_x[i])*deltaT/m
deltaV_y[i] = (Fmag_y[i] - Fdrag_y[i])*deltaT/m
#Velocity
v_x[i] = v_x[i] + deltaV_x[i]
v_y[i] = v_y[i] + deltaV_y[i]
#v_x = moment*30/(6*np.pi*mu*p_rad)*xp[i]/np.abs(xp[i])
#v_y = moment*40/(6*np.pi*mu*p_rad)*yp[i]/np.abs(yp[i])
#New X and Y coordinates of particles
xp[i] = xp[i] + (v_x[i]*deltaT)
yp[i] = yp[i] + (v_y[i]*deltaT)
#Counter
if -yslice < yp[i] < yslice:
count6[timestep] = count6[timestep] + 1
timestep = timestep + 1
plt.plot(xp,yp,".")
plt.xlim(-radius,radius)
plt.ylim(-radius,radius)
ani = animation.FuncAnimation(fig1, update, frames=1, fargs=(plt, fig1), interval=0.1)
#Distribution of particles
#theta = np.arctan2(yp,xp)
#plt.hist(theta,50)
#plt.title('Distribution of Particles Generated')
#plt.xlabel('Theta (Rads)')
#plt.ylabel('Number of Particles')
#plt.show()