我正在编写基于格子玻尔兹曼方法 (LBM) 的 CFD 求解器，以分析方形圆柱体周围的流动。但是，当我尝试编译时，遇到以下错误：
“runfile('C:/Users/sayad/Desktop/sanstitre0.py', wdir='C:/Users/sayad/Desktop') 回溯（最近一次调用最后一次）：
compat_exec 中的文件 D:\Spyder\pkgs\spyder_kernels\py3compat.py:356 exec（代码，全局变量，局部变量）
文件 c:\users\sayad\desktop\sanstitre0.py:106 geq,rho,ux,uy=init(M0=0.3)
init 中的文件 c:\users\sayad\desktop\sanstitre0.py:41 eq(geq,rho,ux,uy)
eq中的文件c:\users\sayad\desktop\sanstitre0.py:72 geq[:, :, 1:9] = w[1:] * rho * (1 + (c0**(-2)) * (ca[1:9, 0]*ux + ca[1:9, 1]uy) + 0.5 (c0-4) * (ca[1:9, 0]*ux + ca[1:9, 1]*uy)2 - 0.5 * (c0(-2) ) * (ux2 + uy**2))
ValueError：操作数无法与形状 (8,) (80,40) 一起广播“
这是我的代码：
import numpy as np
import matplotlib.pyplot as plt
import time
# Parametres du modele D2Q9
Re= 20 # Nombre de Reynolds (à faire varier)
ca= np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1], [1, 1], [-1, 1], [-1, -1], [1, -1]]) # vitesses discrètes
op= -ca # symetrique de chaque vitesse discrète
w= np.array([4/9, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36]) # coefficients omega_alpha
c0= 1/np.sqrt(3) # coefficient de vitesse du son
D= 4 # Resolution du cylindre (nombre de point dans le diamètre)
nx,ny= 20*D, 10*D
M0= 0.3
taug= (M0*D)/(c0*Re) + 1 # Temps de relaxation du modèle (défini en fonction du Reynolds)
nt= int((150*D)/(M0*c0)) # Nombre d'itération a définir
mask=np.zeros((nx,ny));
# Definition du cylindre => Choisir les bons indices pour placer
# le cylindre au centre du domaine en y et dans le premier tiers en x
cx1,cx2= int(8*D - 0.5*D), int(8*D + 0.5*D)
cy1,cy2= int(5*D - 0.5*D), int(5*D + 0.5*D)
mask[cx1:cx2,cy1:cy2]=1
def init(M0):
# fonction d'initialisation des variables macroscopiques
# et de la fonction d'équilibre
rho= np.ones((nx, ny)) * c0**2
ux= np.full((nx, ny), M0 * c0)
uy= np.zeros((nx, ny))
geq=np.zeros((nx,ny,9))
eq(geq,rho,ux,uy)
return geq,rho,ux,uy
def collide(gcoll,g,geq,taug):
# Etape de collision
gcoll[:,:,:]= g[:,:,:] - (1/taug)*(g[:,:,:] - geq[:,:,:])
def propagate(g,gcoll):
# Etape de propagation
g[:, :, 0] = gcoll[:, :, 0]
g[:, :, 1] = gcoll[:, :, 1]
g[:, :, 2] = gcoll[:, :, 2]
g[:, :, 3] = gcoll[:, :, 3]
g[:, :, 4] = gcoll[:, :, 4]
g[:, :, 5] = gcoll[:, :, 5]
g[:, :, 6] = gcoll[:, :, 6]
g[:, :, 7] = gcoll[:, :, 7]
g[:, :, 8] = gcoll[:, :, 8]
def macro(g,rho,ux,uy):
# calcul des variables macro
rho[:, :] = np.sum(g, axis=2)
ux[:, :] = (g[:, :, 1] - g[:, :, 3] + g[:, :, 5] - g[:, :, 6] - g[:, :, 7] + g[:, :, 8]) / rho[:, :]
uy[:, :] = (g[:, :, 2] - g[:, :, 4] + g[:, :, 5] + g[:, :, 6] - g[:, :, 7] - g[:, :, 8]) / rho[:, :]
def eq(geq,rho,ux,uy):
# Calcul de la fonction d'équilibre
geq[:, :, 0] = w[0] * rho * (1 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))
geq[:, :, 1:9] = w[1:] * rho * (1 + (c0**(-2)) * (ca[1:9, 0]*ux + ca[1:9, 1]*uy) + 0.5* (c0**-4) * (ca[1:9, 0]*ux + ca[1:9, 1]*uy)**2 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))
def boundary(g,gcoll):
# Conditions aux limites periodiques
# Entree => quelles sont les g inconnues en entree
g[0, :, 0] = gcoll[0, :, 0]
g[0, :, 1:8] = gcoll[0, :, 1:8]
g[0, :, 8] = gcoll[0, :, 8]
# Sortie => quelles sont les g inconnues en sortie
g[-1, :, 0] = gcoll[-1, :, 0]
g[-1, :, 1:8] = gcoll[-1, :, 1:8]
g[-1, :, 8] = gcoll[-1, :, 8]
# Bas => quelles sont les g inconnues en bas
g[:, 0, 0] = gcoll[:, 0, 0]
g[:, 0, 1:8] = gcoll[:, 0, 1:8]
g[:, 0, 8] = gcoll[:, 0, 8]
# Haut => quelles sont les g inconnues en haut
g[:, -1, 0] = gcoll[:, -1, 0]
g[:, -1, 1:8] = gcoll[:, -1, 1:8]
g[:, -1, 8] = gcoll[:, -1, 8]
def wall(g,mask):
# Conditions de parois
g[mask == 1, 1] = g[mask == 1, 3]
g[mask == 1, 2] = g[mask == 1, 4]
g[mask == 1, 3] = g[mask == 1, 1]
g[mask == 1, 4] = g[mask == 1, 2]
g[mask == 1, 5] = g[mask == 1, 7]
g[mask == 1, 6] = g[mask == 1, 8]
g[mask == 1, 7] = g[mask == 1, 5]
g[mask == 1, 8] = g[mask == 1, 6]
#Calcul
geq,rho,ux,uy=init(M0=0.3)
g,gcoll=geq.copy(),geq.copy()
prof=np.zeros(nt) #profil de pression en fonction du temps
start=time.time()
for t in range(nt):
collide(gcoll,g,geq,taug) # Collision
propagate(g,gcoll) # Propagation
boundary(g,gcoll) # Conditions aux limites domaine
wall(g,mask) # Conditions aux limites solides
macro(g,rho,ux,uy) # Calcul des moments
eq(geq,rho,ux,uy) # Calcul de l'equilibre
prof[t]= rho[((cx1 + cx2) // 2) + 5*D, ((cy1 + cy2) // 2) + D] * c0**2 # Ecriture des résultats
tcal=time.time()-start
print("Temps de calcul: "+str(tcal))
我尝试改变 rho、ux 和 uy 的尺寸。但一切都没有真正改变。
编辑：@Marcel 建议在 init 函数中调用 eq(geq, rho, ux, uy) 。但现在我在 eq 函数的第二行遇到错误，如上所示。
假设你的数学是正确的，你有两个错误：
首先，您使用错误的参数调用函数
eq其次，您需要将传入数组的形状扩展为
eqdef eq(geq,rho,ux,uy):
# Calcul de la fonction d'équilibre
uxb = ux[:, :, None]
uyb = uy[:, :, None]
wb = w[None, None, :]
cab = ca[None, None, 1:9, :]
rhob = rho[:, :, None]
geq[:, :, 0] = w[0] * rho * (1 - 0.5 * (c0**(-2)) * (ux**2 + uy**2))
geq[:, :, 1:9] = wb[:, :, 1:] * (rhob * (1 + (c0**(-2)) * (cab[..., 0]*uxb + cab[..., 1]*uyb) + 0.5* (c0**-4) * (cab[..., 0]*uxb + cab[..., 1]*uyb)**2 - 0.5 * (c0**(-2)) * (uxb**2 + uyb**2)))
这里的“b”前缀，我的意思是“用于广播”。