基于Python求解任意偶极子强度方程组以实现NACA0018翼型流场可视化

问题需求 #

刚接触Python，需完成流体力学作业：实现均匀流中偶极子片的流函数可视化，通过求解线性方程组确定任意数量偶极子的强度，使翼型表面流线近似NACA0018翼型（翼型表面总流函数为0）。

核心公式 #

流函数表达式 #

• 均匀流流函数：Psi_uniform = U_inf * r * sin(theta)
• 单个偶极子流函数：Psi_doublet = -(k/(2*pi)) * (sin(theta_ij)/r_ij)，其中r_ij是计算点到偶极子j的距离，theta_ij是对应极角

线性方程组构建 #

已知翼型表面n个点的坐标（已转换为极坐标或保留笛卡尔坐标），每个点满足总流函数为0：

Psi_uniform(i) + Σ(Psi_doublet_j(i)) = 0 （j从1到n）

展开整理为矩阵形式A*K = B：

• K是待求偶极子强度数组[k₁, k₂, ..., kₙ]
• 矩阵A的元素A[i][j] = -sin(theta_ij)/(2*pi*r_ij)，表示偶极子j在点i处的贡献系数
• 向量B的元素B[i] = -U_inf * r_i * sin(theta_i)，对应均匀流在点i处的流函数负值

Python实现步骤 #

1. 导入依赖库 #

import numpy as np
import matplotlib.pyplot as plt

2. 定义基础参数与翼型点数据 #

生成NACA0018翼型的离散点（或读取外部翼型数据）：

U_inf = 1.0  # 来流速度
n = 60  # 偶极子数量（对应翼型离散点数量）

# 生成NACA0018翼型点（标准四位翼型厚度公式）
x = np.linspace(-0.5, 0.5, n)
t = 0.18  # NACA0018厚度系数
y = 5 * t * (0.2969*np.sqrt(np.abs(x)) - 0.1260*x - 0.3516*x**2 + 0.2843*x**3 - 0.1015*x**4)
# 合并上下翼面（可选，示例仅用上翼面，实际可对称生成下翼面）
x = np.hstack([x, x[::-1]])
y = np.hstack([y, -y[::-1]])
n = len(x)  # 更新点数量

3. 构建线性方程组的矩阵A和向量B #

A = np.zeros((n, n))
B = np.zeros(n)

for i in range(n):
    # 计算点i的极坐标参数
    r_i = np.sqrt(x[i]**2 + y[i]**2)
    theta_i = np.arctan2(y[i], x[i])
    B[i] = -U_inf * r_i * np.sin(theta_i)
    
    for j in range(n):
        # 计算点i到偶极子j的距离和极角
        dx = x[i] - x[j]
        dy = y[i] - y[j]
        r_ij = np.sqrt(dx**2 + dy**2)
        theta_ij = np.arctan2(dy, dx)
        # 避免除以0（偶极子自身位置点）
        if r_ij < 1e-6:
            A[i][j] = 0.0
        else:
            A[i][j] = -np.sin(theta_ij) / (2 * np.pi * r_ij)

4. 求解偶极子强度 #

若矩阵非奇异，直接用精确求解；若奇异（如点重合），用最小二乘求解：

try:
    K = np.linalg.solve(A, B)
except np.linalg.LinAlgError:
    # 最小二乘近似求解
    K, _, _, _ = np.linalg.lstsq(A, B, rcond=None)

5. 流函数可视化 #

生成计算域网格，计算总流函数并绘制流线：

# 生成计算域网格
x_grid, y_grid = np.meshgrid(np.linspace(-2, 2, 150), np.linspace(-1.5, 1.5, 150))

# 计算均匀流流函数
r_grid = np.sqrt(x_grid**2 + y_grid**2)
theta_grid = np.arctan2(y_grid, x_grid)
Psi_uniform = U_inf * r_grid * np.sin(theta_grid)

# 计算所有偶极子的流函数叠加
Psi_doublet = np.zeros_like(Psi_uniform)
for j in range(n):
    dx = x_grid - x[j]
    dy = y_grid - y[j]
    r_ij = np.sqrt(dx**2 + dy**2)
    theta_ij = np.arctan2(dy, dx)
    # 避免奇点处数值异常
    mask = r_ij > 1e-6
    Psi_doublet[mask] -= (K[j]/(2*np.pi)) * np.sin(theta_ij[mask]) / r_ij[mask]

# 总流函数
Psi_total = Psi_uniform + Psi_doublet

# 绘制流线与翼型
plt.figure(figsize=(12, 6))
contour = plt.contour(x_grid, y_grid, Psi_total, levels=np.linspace(-2, 2, 60), cmap='coolwarm')
plt.scatter(x, y, color='black', s=10, label='NACA0018翼型')
plt.clabel(contour, inline=True, fontsize=8)
plt.xlabel('x')
plt.ylabel('y')
plt.title('均匀流+偶极子片流函数流线分布')
plt.legend()
plt.axis('equal')
plt.show()

关键注意事项 #

• 偶极子需布置在翼型离散点处，计算流函数时要使用计算点到偶极子的相对距离与角度，而非计算点自身的极坐标参数（这是避免矩阵奇异的关键）
• 若翼型点数量过多，可通过稀疏矩阵优化计算效率，作业级别用常规矩阵即可
• NACA0018翼型点也可通过专业气动库（如pyavl）读取更精确的离散数据

内容的提问来源于stack exchange，提问作者AeroEng43