不可压缩流体有限元求解器

项目概述

这是有限元方法项目的一个专门模块，专注于不可压缩流体的数值求解。不可压缩流体在低速流动、传热传质、生物流体力学等领域应用广泛，其求解需要考虑压力-速度耦合、质量守恒等特殊处理。

理论基础

不可压缩Navier-Stokes方程

不可压缩流体的控制方程基于不可压缩Navier-Stokes方程组：

连续性方程： $$\nabla \cdot \mathbf{u} = 0$$

动量方程： $$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}$$

能量方程（传热）： $$\frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \alpha \nabla^2 T + \frac{Q}{\rho c_p}$$

其中： - $\mathbf{u}$ 是速度向量 - $p$ 是压力 - $\rho$ 是流体密度（常数） - $\nu$ 是运动粘度 - $T$ 是温度 - $\alpha$ 是热扩散率

压力-速度耦合

不可压缩流体的核心挑战是压力-速度耦合问题。常用的求解方法包括：

分步法（Fractional Step Method）： 1. 预测步：求解中间速度 2. 压力步：求解压力修正 3. 修正步：更新速度场

SIMPLE算法： - 压力修正方程 - 速度修正方程 - 迭代求解直到收敛

数值方法

有限元离散化

弱形式： 对于动量方程，弱形式为： $$\int_\Omega \frac{\partial \mathbf{u}}{\partial t} \cdot \mathbf{v} d\Omega + \int_\Omega (\mathbf{u} \cdot \nabla)\mathbf{u} \cdot \mathbf{v} d\Omega$$ $$= -\int_\Omega \frac{1}{\rho}\nabla p \cdot \mathbf{v} d\Omega + \int_\Omega \nu \nabla^2 \mathbf{u} \cdot \mathbf{v} d\Omega$$

Galerkin方法： $$\mathbf{u} = \sum_{i=1}^{n} u_i \phi_i, \quad p = \sum_{i=1}^{m} p_i \psi_i$$

其中 $\phi_i$ 是速度形函数，$\psi_i$ 是压力形函数。

时间离散化

向后欧拉方法： $$\frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} + (\mathbf{u}^{n+1} \cdot \nabla)\mathbf{u}^{n+1} = -\frac{1}{\rho}\nabla p^{n+1} + \nu \nabla^2 \mathbf{u}^{n+1}$$

Crank-Nicolson方法： $$\frac{\mathbf{u}^{n+1} - \mathbf{u}^n}{\Delta t} + \frac{1}{2}[(\mathbf{u}^{n+1} \cdot \nabla)\mathbf{u}^{n+1} + (\mathbf{u}^n \cdot \nabla)\mathbf{u}^n]$$ $$= -\frac{1}{\rho}\nabla p^{n+1} + \frac{\nu}{2}[\nabla^2 \mathbf{u}^{n+1} + \nabla^2 \mathbf{u}^n]$$

算法实现

分步法实现

class IncompressibleSolver {
private:
    Mesh mesh;
    Material material;
    LinearSolver linearSolver;

public:
    void solve() {
        initialize();

        while (!converged()) {
            // 预测步：求解中间速度
            VectorXd u_star = predictVelocity();

            // 压力步：求解压力修正
            VectorXd p_correction = solvePressureCorrection(u_star);

            // 修正步：更新速度场
            updateVelocity(u_star, p_correction);

            // 更新压力
            updatePressure(p_correction);

            // 检查收敛性
            checkConvergence();
        }
    }

private:
    VectorXd predictVelocity() {
        // 求解动量方程（忽略压力梯度）
        MatrixXd A = assembleMomentumMatrix();
        VectorXd b = assembleMomentumRHS();

        // 边界条件处理
        applyVelocityBoundaryConditions(A, b);

        return linearSolver.solve(A, b);
    }

    VectorXd solvePressureCorrection(const VectorXd& u_star) {
        // 求解压力修正方程
        MatrixXd A_p = assemblePressureMatrix();
        VectorXd b_p = assemblePressureRHS(u_star);

        // 边界条件处理
        applyPressureBoundaryConditions(A_p, b_p);

        return linearSolver.solve(A_p, b_p);
    }

    void updateVelocity(const VectorXd& u_star, 
                       const VectorXd& p_correction) {
        // 速度修正：u = u* - (Δt/ρ)∇p'
        for (int i = 0; i < mesh.nodes.size(); i++) {
            VectorXd grad_p = calculatePressureGradient(i, p_correction);
            velocity[i] = u_star[i] - (dt / density) * grad_p;
        }
    }
};

SIMPLE算法实现

class SIMPLEAlgorithm {
private:
    double relaxation_factor_u = 0.7;
    double relaxation_factor_p = 0.3;

public:
    void solve() {
        initialize();

        for (int iteration = 0; iteration < max_iterations; iteration++) {
            // 求解动量方程
            solveMomentumEquations();

            // 计算质量流量
            calculateMassFlux();

            // 求解压力修正方程
            solvePressureCorrection();

            // 修正速度和压力
            correctVelocityAndPressure();

            // 检查收敛性
            if (checkConvergence()) break;
        }
    }

private:
    void solveMomentumEquations() {
        // 求解u动量方程
        MatrixXd A_u = assembleUMomentumMatrix();
        VectorXd b_u = assembleUMomentumRHS();
        VectorXd u_new = linearSolver.solve(A_u, b_u);

        // 松弛更新
        u = relaxation_factor_u * u_new + (1 - relaxation_factor_u) * u;

        // 求解v动量方程
        MatrixXd A_v = assembleVMomentumMatrix();
        VectorXd b_v = assembleVMomentumRHS();
        VectorXd v_new = linearSolver.solve(A_v, b_v);

        v = relaxation_factor_u * v_new + (1 - relaxation_factor_u) * v;
    }

    void solvePressureCorrection() {
        // 组装压力修正方程
        MatrixXd A_p = assemblePressureCorrectionMatrix();
        VectorXd b_p = assemblePressureCorrectionRHS();

        // 求解压力修正
        VectorXd p_correction = linearSolver.solve(A_p, b_p);

        // 松弛更新压力
        p = p + relaxation_factor_p * p_correction;
    }
};

应用案例

方腔流动

问题描述: 二维方腔内的自然对流 - 瑞利数: Ra = 10^6 - 普朗特数: Pr = 0.71 - 几何: 1×1 正方形腔体

数值结果: - 流线分布 - 等温线分布 - 努塞尔数计算

管道流动

问题描述: 圆管内的层流和湍流 - 雷诺数: Re = 100-10000 - 几何: 圆形管道

求解结果: - 速度剖面 - 压力分布 - 摩擦系数

传热问题

问题描述: 强制对流换热 - 雷诺数: Re = 1000 - 普朗特数: Pr = 0.7 - 几何: 圆柱绕流

结果分析: - 温度场分布 - 局部努塞尔数 - 平均换热系数

性能优化

线性求解器优化

class OptimizedLinearSolver {
private:
    // 预处理共轭梯度法
    PreconditionedCG pcg_solver;

    // 多重网格法
    MultigridSolver mg_solver;

public:
    VectorXd solve(const MatrixXd& A, const VectorXd& b) {
        if (A.rows() < 10000) {
            // 小规模问题使用直接法
            return directSolver.solve(A, b);
        } else {
            // 大规模问题使用迭代法
            return pcg_solver.solve(A, b);
        }
    }

private:
    void setupPreconditioner(const MatrixXd& A) {
        // 设置ILU预处理
        ilu_preconditioner.compute(A);
        pcg_solver.setPreconditioner(ilu_preconditioner);
    }
};

并行计算

// OpenMP并行化
void parallelAssemble() {
    #pragma omp parallel for
    for (int i = 0; i < elements.size(); i++) {
        elements[i].assembleElementMatrix();
    }

    // 合并全局矩阵
    assembleGlobalMatrix();
}

// MPI分布式计算
void distributedSolve() {
    // 网格分区
    partitionMesh();

    // 通信边界处理
    setupCommunicationBoundaries();

    // 分布式求解
    solveDistributedSystem();
}

验证与测试

解析解对比

Poiseuille流动: - 圆管内的层流解析解 - 速度剖面验证 - 压力梯度验证

Couette流动: - 两平行板间的剪切流动 - 线性速度分布验证

基准问题

方腔流动基准: - 与Ghia等人结果对比 - 中心线速度分布 - 涡心位置验证

后台阶流动: - 与Arnal等人结果对比 - 再附着长度验证 - 压力系数分布

相关链接

这个模块展示了我在不可压缩流体数值计算和压力-速度耦合算法方面的专业能力。

返回项目列表