Python求解矩形通道二维流速方程报MemoryError的解决方法

内存错误的核心诱因不是双曲函数运算开销，是代码逻辑存在严重的重复计算和符号运算滥用：

• 重复计算量级爆炸：solution_func的每轮循环都会完整重跑Q_sheath_func()和Q_total_func()，而这两个函数内部又会反复调用u1()、sumfunc()，重复执行双变量符号积分超过20次。符号积分生成的表达式树本身长度随积分次数指数增长，重复生成直接占满内存。
• 数值/符号混用导致表达式膨胀：代码中把np.pi这类浮点数值直接传入sympy符号运算流程，导致sympy无法做自动化简，生成的中间表达式体积比用符号常量sy.pi大数十倍。
• 无意义保留全符号变量：积分变量y、z的上下限除deld外全部为固定常数，完全不需要在循环中反复对y、z做积分，可提前一次性积分消去两个坐标变量，得到仅含deld的单变量表达式。
• 符号求解器适配错误：最终待解方程含双曲函数，属于超越方程，sy.solve会尝试做符号解析推导，生成巨量无效中间表达式，根本不适合求解这类问题。

核心优化逻辑：仅做1次通用符号积分消去y、z，提前缓存所有固定计算结果，用数值求解器替代符号求解器计算deld，内存占用可从百GB级降至几十MB级。

import numpy as np
import sympy as sy
from scipy.optimize import fsolve

# 通道参数
H = 0.00005  # m 通道高度
w = 0.00002  # m 通道宽度
L = 0.015    # m 通道长度

# 流体参数
eta = 0.001041  # Pas 聚合物溶液黏度

# 符号定义
deld = sy.Symbol('deld', real=True)
y = sy.Symbol('y', real=True)
z = sy.Symbol('z', real=True)

# 提前计算固定的u1系数列表，避免重复计算
def calc_u1_list():
    u1_list = []
    Q_sample = 1*(10**-9)/60
    for i in range(1, 11):
        Q_total = Q_sample * (1 + i)
        delP = (((8*eta*L)*(w+H)**2)/((w**3)*(H**3))) * Q_total
        u1_caculation = (4*(H**2)*delP)/(sy.pi*eta*L)
        u1_list.append(float(u1_caculation))
    return u1_list
u1_list = calc_u1_list()

# 提前计算速度剖面的级数表达式，仅生成1次
def calc_sum_func():
    summation_all = 0
    for k in range(0, 3): # 可自行调整级数项数量
        n = 2*k + 1
        term = (n**(-3)) * (1 - (sy.cosh(n*sy.pi*y/H)/sy.cosh(n*sy.pi*w/(2*H)))) * sy.sin(n*sy.pi*z/H)
        summation_all += term
    return summation_all
u_profile_base = calc_sum_func()

# 提前做1次通用积分，得到鞘层流量关于deld和u系数的表达式，避免重复积分
def calc_sheath_Q_template():
    velocity = u_profile_base
    # 对y积分
    int_y = sy.integrate(velocity, y)
    int_y_sub = int_y.subs(y, w/2 - deld) - int_y.subs(y, -w/2 + deld)
    # 对z积分
    int_z = sy.integrate(int_y_sub, z)
    Q_sheath_template = int_z.subs(z, H) - int_z.subs(z, 0)
    return sy.lambdify((deld,), Q_sheath_template, modules='numpy')
sheath_Q_calc = calc_sheath_Q_template()

# 提前算好总流量（固定值，和deld无关）
def calc_total_Q():
    velocity = u_profile_base
    int_y = sy.integrate(velocity, y)
    int_y_sub = int_y.subs(y, w/2) - int_y.subs(y, -w/2)
    int_z = sy.integrate(int_y_sub, z)
    Q_total_template = int_z.subs(z, H) - int_z.subs(z, 0)
    # 总流量是常数，直接求值
    return float(Q_total_template)
Q_total_base = calc_total_Q()

# 数值求解每个流量比对应的deld
if __name__ == "__main__":
    solution = []
    for idx, u_coeff in enumerate(u1_list):
        a = idx + 1 # 流量比
        Q_total = Q_total_base * u_coeff
        # 构造待解方程 f(deld)=0
        def equation(d_val):
            Q_sheath = sheath_Q_calc(d_val[0]) * u_coeff
            return Q_sheath/(Q_total - Q_sheath) - a
        # 给初始猜测值（取deld=w/4作为初值，可根据物理场景调整）
        d_init = [w/4]
        res = fsolve(equation, d_init)
        solution.append(float(res[0]))
    print("各流量比对应的deld计算结果：", solution)

• 代码中保留了你原有的级数展开逻辑，可通过调整range(0,3)的上限增加级数项提升精度，不会出现内存溢出
• 数值求解的初值d_init可根据实际物理场景调整，避免求解落到非物理解
• 所有固定参数的计算全部提前缓存，没有冗余重复运算，普通8GB内存设备即可秒级完成计算

内容的提问来源于stack exchange，提问作者전현우