显式求解三维温度场并输出Tecplot文件的C++实现

1. 主程序头文件和参数定义

// thermal_solver_3d.h
#ifndef THERMAL_SOLVER_3D_H
#define THERMAL_SOLVER_3D_H

#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <cmath>

class ThermalSolver3D {
private:
    // 网格参数
    int nx, ny, nz;           // 网格点数
    double dx, dy, dz;        // 网格间距
    double Lx, Ly, Lz;        // 计算域尺寸
    
    // 材料参数
    double k;                 // 热导率 W/(m·K)
    double rho;               // 密度 kg/m³
    double cp;                // 比热容 J/(kg·K)
    double alpha;             // 热扩散系数 m²/s
    
    // 时间参数
    double dt;                // 时间步长
    double total_time;        // 总模拟时间
    double current_time;      // 当前时间
    
    // 稳定性条件
    double Fo;                // 傅里叶数
    
    // 温度场
    std::vector<std::vector<std::vector<double>>> T;      // 当前时间步温度
    std::vector<std::vector<std::vector<double>>> T_old;  // 上一时间步温度
    
    // 源项
    std::vector<std::vector<std::vector<double>>> Q;      // 热源 W/m³
    
    // 边界条件类型
    enum BCType {
        DIRICHLET,    // 第一类边界条件（固定温度）
        NEUMANN,      // 第二类边界条件（固定热流）
        ADIABATIC     // 绝热边界
    };
    
    BCType bc_x_min, bc_x_max;
    BCType bc_y_min, bc_y_max; 
    BCType bc_z_min, bc_z_max;
    
    // 边界值
    double bc_value_x_min, bc_value_x_max;
    double bc_value_y_min, bc_value_y_max;
    double bc_value_z_min, bc_value_z_max;

public:
    ThermalSolver3D(int nx, int ny, int nz, double Lx, double Ly, double Lz);
    ~ThermalSolver3D();
    
    // 设置材料属性
    void setMaterialProperties(double thermal_conductivity, double density, double specific_heat);
    
    // 设置时间参数
    void setTimeParameters(double time_step, double total_simulation_time);
    
    // 设置边界条件
    void setBoundaryConditions(BCType x_min, BCType x_max, BCType y_min, BCType y_max, 
                              BCType z_min, BCType z_max, double value_min = 0, double value_max = 0);
    
    // 设置初始条件
    void setInitialCondition(double initial_temp);
    void setInitialConditionGaussian(double center_x, double center_y, double center_z, 
                                   double amplitude, double spread);
    
    // 设置热源
    void setHeatSource(double power);
    void setLocalHeatSource(int i_min, int i_max, int j_min, int j_max, int k_min, int k_max, double power_density);
    
    // 显式求解
    void solveExplicit();
    
    // 输出结果
    void writeTecplotFile(const std::string& filename);
    void writeTecplotHeader(std::ofstream& file);
    void writeTecplotData(std::ofstream& file);
    
    // 获取温度场
    double getTemperature(int i, int j, int k) const;
    
    // 稳定性检查
    bool checkStability() const;
    
private:
    // 应用边界条件
    void applyBoundaryConditions();
    
    // 计算内部节点
    void updateInternalNodes();
    
    // 初始化内存
    void initializeArrays();
};

#endif

2. 主程序实现

// thermal_solver_3d.cpp
#include "thermal_solver_3d.h"

// 构造函数
ThermalSolver3D::ThermalSolver3D(int nx, int ny, int nz, double Lx, double Ly, double Lz)
    : nx(nx), ny(ny), nz(nz), Lx(Lx), Ly(Ly), Lz(Lz), current_time(0.0) {
    
    // 计算网格间距
    dx = Lx / (nx - 1);
    dy = Ly / (ny - 1);
    dz = Lz / (nz - 1);
    
    // 默认材料属性（铜）
    k = 400.0;      // W/(m·K)
    rho = 8960.0;   // kg/m³
    cp = 385.0;     // J/(kg·K)
    alpha = k / (rho * cp);
    
    // 默认边界条件（绝热）
    bc_x_min = bc_x_max = bc_y_min = bc_y_max = bc_z_min = bc_z_max = ADIABATIC;
    
    // 初始化数组
    initializeArrays();
}

// 析构函数
ThermalSolver3D::~ThermalSolver3D() {
    // 向量会自动清理
}

// 初始化数组
void ThermalSolver3D::initializeArrays() {
    // 调整温度场大小
    T.resize(nx);
    T_old.resize(nx);
    Q.resize(nx);
    
    for (int i = 0; i < nx; i++) {
        T[i].resize(ny);
        T_old[i].resize(ny);
        Q[i].resize(ny);
        
        for (int j = 0; j < ny; j++) {
            T[i][j].resize(nz, 0.0);
            T_old[i][j].resize(nz, 0.0);
            Q[i][j].resize(nz, 0.0);
        }
    }
}

// 设置材料属性
void ThermalSolver3D::setMaterialProperties(double thermal_conductivity, double density, double specific_heat) {
    k = thermal_conductivity;
    rho = density;
    cp = specific_heat;
    alpha = k / (rho * cp);
}

// 设置时间参数
void ThermalSolver3D::setTimeParameters(double time_step, double total_simulation_time) {
    dt = time_step;
    total_time = total_simulation_time;
    
    // 计算傅里叶数用于稳定性检查
    double dx2 = dx * dx;
    double dy2 = dy * dy;
    double dz2 = dz * dz;
    Fo = alpha * dt * (1.0/dx2 + 1.0/dy2 + 1.0/dz2);
}

// 设置边界条件
void ThermalSolver3D::setBoundaryConditions(BCType x_min, BCType x_max, BCType y_min, BCType y_max,
                                          BCType z_min, BCType z_max, double value_min, double value_max) {
    bc_x_min = x_min; bc_x_max = x_max;
    bc_y_min = y_min; bc_y_max = y_max;
    bc_z_min = z_min; bc_z_max = z_max;
    
    // 对于Dirichlet边界，设置边界值
    if (x_min == DIRICHLET) bc_value_x_min = value_min;
    if (x_max == DIRICHLET) bc_value_x_max = value_max;
    if (y_min == DIRICHLET) bc_value_y_min = value_min;
    if (y_max == DIRICHLET) bc_value_y_max = value_max;
    if (z_min == DIRICHLET) bc_value_z_min = value_min;
    if (z_max == DIRICHLET) bc_value_z_max = value_max;
}

// 设置均匀初始条件
void ThermalSolver3D::setInitialCondition(double initial_temp) {
    for (int i = 0; i < nx; i++) {
        for (int j = 0; j < ny; j++) {
            for (int k = 0; k < nz; k++) {
                T[i][j][k] = initial_temp;
                T_old[i][j][k] = initial_temp;
            }
        }
    }
}

// 设置高斯分布初始条件
void ThermalSolver3D::setInitialConditionGaussian(double center_x, double center_y, double center_z,
                                                 double amplitude, double spread) {
    for (int i = 0; i < nx; i++) {
        double x = i * dx;
        for (int j = 0; j < ny; j++) {
            double y = j * dy;
            for (int k = 0; k < nz; k++) {
                double z = k * dz;
                
                double rx = x - center_x;
                double ry = y - center_y;
                double rz = z - center_z;
                double r2 = rx*rx + ry*ry + rz*rz;
                
                T[i][j][k] = amplitude * exp(-r2 / (2.0 * spread * spread));
                T_old[i][j][k] = T[i][j][k];
            }
        }
    }
}

// 设置均匀热源
void ThermalSolver3D::setHeatSource(double power) {
    double volume = dx * dy * dz;
    double power_density = power / (nx * ny * nz * volume);
    
    for (int i = 0; i < nx; i++) {
        for (int j = 0; j < ny; j++) {
            for (int k = 0; k < nz; k++) {
                Q[i][j][k] = power_density;
            }
        }
    }
}

// 设置局部热源
void ThermalSolver3D::setLocalHeatSource(int i_min, int i_max, int j_min, int j_max, int k_min, int k_max, double power_density) {
    for (int i = i_min; i <= i_max && i < nx; i++) {
        for (int j = j_min; j <= j_max && j < ny; j++) {
            for (int k = k_min; k <= k_max && k < nz; k++) {
                Q[i][j][k] = power_density;
            }
        }
    }
}

// 检查稳定性条件
bool ThermalSolver3D::checkStability() const {
    double stability_limit = 1.0 / 6.0;  // 三维显式格式稳定性限制
    
    if (Fo > stability_limit) {
        std::cerr << "警告：稳定性条件不满足！Fo = " << Fo 
                  << " > " << stability_limit << std::endl;
        return false;
    }
    return true;
}

// 应用边界条件
void ThermalSolver3D::applyBoundaryConditions() {
    // X方向边界
    for (int j = 0; j < ny; j++) {
        for (int k = 0; k < nz; k++) {
            // X最小边界
            switch (bc_x_min) {
                case DIRICHLET:
                    T[0][j][k] = bc_value_x_min;
                    break;
                case NEUMANN:
                    // 使用前向差分
                    T[0][j][k] = T[1][j][k] - bc_value_x_min * dx / k;
                    break;
                case ADIABATIC:
                    T[0][j][k] = T[1][j][k];
                    break;
            }
            
            // X最大边界
            switch (bc_x_max) {
                case DIRICHLET:
                    T[nx-1][j][k] = bc_value_x_max;
                    break;
                case NEUMANN:
                    // 使用后向差分
                    T[nx-1][j][k] = T[nx-2][j][k] + bc_value_x_max * dx / k;
                    break;
                case ADIABATIC:
                    T[nx-1][j][k] = T[nx-2][j][k];
                    break;
            }
        }
    }
    
    // Y方向边界
    for (int i = 0; i < nx; i++) {
        for (int k = 0; k < nz; k++) {
            // Y最小边界
            switch (bc_y_min) {
                case DIRICHLET:
                    T[i][0][k] = bc_value_y_min;
                    break;
                case NEUMANN:
                    T[i][0][k] = T[i][1][k] - bc_value_y_min * dy / k;
                    break;
                case ADIABATIC:
                    T[i][0][k] = T[i][1][k];
                    break;
            }
            
            // Y最大边界
            switch (bc_y_max) {
                case DIRICHLET:
                    T[i][ny-1][k] = bc_value_y_max;
                    break;
                case NEUMANN:
                    T[i][ny-1][k] = T[i][ny-2][k] + bc_value_y_max * dy / k;
                    break;
                case ADIABATIC:
                    T[i][ny-1][k] = T[i][ny-2][k];
                    break;
            }
        }
    }
    
    // Z方向边界
    for (int i = 0; i < nx; i++) {
        for (int j = 0; j < ny; j++) {
            // Z最小边界
            switch (bc_z_min) {
                case DIRICHLET:
                    T[i][j][0] = bc_value_z_min;
                    break;
                case NEUMANN:
                    T[i][j][0] = T[i][j][1] - bc_value_z_min * dz / k;
                    break;
                case ADIABATIC:
                    T[i][j][0] = T[i][j][1];
                    break;
            }
            
            // Z最大边界
            switch (bc_z_max) {
                case DIRICHLET:
                    T[i][j][nz-1] = bc_value_z_max;
                    break;
                case NEUMANN:
                    T[i][j][nz-1] = T[i][j][nz-2] + bc_value_z_max * dz / k;
                    break;
                case ADIABATIC:
                    T[i][j][nz-1] = T[i][j][nz-2];
                    break;
            }
        }
    }
}

// 更新内部节点温度
void ThermalSolver3D::updateInternalNodes() {
    double dx2 = dx * dx;
    double dy2 = dy * dy;
    double dz2 = dz * dz;
    
    double factor_x = alpha * dt / dx2;
    double factor_y = alpha * dt / dy2;
    double factor_z = alpha * dt / dz2;
    double factor_source = dt / (rho * cp);
    
    // 更新内部节点 (1 to n-2 in each direction)
    for (int i = 1; i < nx - 1; i++) {
        for (int j = 1; j < ny - 1; j++) {
            for (int k = 1; k < nz - 1; k++) {
                double laplacian = 
                    (T_old[i+1][j][k] - 2.0 * T_old[i][j][k] + T_old[i-1][j][k]) / dx2 +
                    (T_old[i][j+1][k] - 2.0 * T_old[i][j][k] + T_old[i][j-1][k]) / dy2 +
                    (T_old[i][j][k+1] - 2.0 * T_old[i][j][k] + T_old[i][j][k-1]) / dz2;
                
                T[i][j][k] = T_old[i][j][k] + alpha * dt * laplacian + factor_source * Q[i][j][k];
            }
        }
    }
}

// 显式求解主函数
void ThermalSolver3D::solveExplicit() {
    if (!checkStability()) {
        std::cerr << "稳定性检查失败，计算可能发散！" << std::endl;
        return;
    }
    
    int time_step = 0;
    double output_interval = total_time / 10.0;  // 输出10个时间点的结果
    
    std::cout << "开始显式求解三维温度场..." << std::endl;
    std::cout << "时间步长: " << dt << " s" << std::endl;
    std::cout << "总模拟时间: " << total_time << " s" << std::endl;
    std::cout << "傅里叶数 Fo = " << Fo << std::endl;
    
    while (current_time < total_time) {
        // 保存上一时间步的温度
        for (int i = 0; i < nx; i++) {
            for (int j = 0; j < ny; j++) {
                for (int k = 0; k < nz; k++) {
                    T_old[i][j][k] = T[i][j][k];
                }
            }
        }
        
        // 更新内部节点
        updateInternalNodes();
        
        // 应用边界条件
        applyBoundaryConditions();
        
        // 更新时间
        current_time += dt;
        time_step++;
        
        // 输出进度
        if (time_step % 100 == 0 || current_time >= total_time) {
            std::cout << "时间步: " << time_step << ", 当前时间: " << current_time << " s" << std::endl;
        }
        
        // 输出结果文件
        if (fmod(current_time, output_interval) < dt || current_time >= total_time) {
            std::string filename = "temperature_3d_t" + std::to_string(int(current_time*1000)) + ".dat";
            writeTecplotFile(filename);
            std::cout << "输出文件: " << filename << std::endl;
        }
    }
    
    std::cout << "计算完成！总时间步数: " << time_step << std::endl;
}

// 写入Tecplot文件
void ThermalSolver3D::writeTecplotFile(const std::string& filename) {
    std::ofstream file(filename);
    if (!file.is_open()) {
        std::cerr << "无法创建文件: " << filename << std::endl;
        return;
    }
    
    writeTecplotHeader(file);
    writeTecplotData(file);
    file.close();
}

// 写入Tecplot文件头
void ThermalSolver3D::writeTecplotHeader(std::ofstream& file) {
    file << "TITLE = \"3D Temperature Field\"" << std::endl;
    file << "VARIABLES = \"X\", \"Y\", \"Z\", \"Temperature\"" << std::endl;
    file << "ZONE I=" << nx << ", J=" << ny << ", K=" << nz 
         << ", F=POINT, SOLUTIONTIME=" << current_time << std::endl;
}

// 写入Tecplot数据
void ThermalSolver3D::writeTecplotData(std::ofstream& file) {
    file.precision(6);
    file << std::scientific;
    
    for (int k = 0; k < nz; k++) {
        double z = k * dz;
        for (int j = 0; j < ny; j++) {
            double y = j * dy;
            for (int i = 0; i < nx; i++) {
                double x = i * dx;
                file << x << " " << y << " " << z << " " << T[i][j][k] << std::endl;
            }
        }
    }
}

// 获取温度值
double ThermalSolver3D::getTemperature(int i, int j, int k) const {
    if (i >= 0 && i < nx && j >= 0 && j < ny && k >= 0 && k < nz) {
        return T[i][j][k];
    }
    return 0.0;
}

3. 主函数和测试用例

// main.cpp
#include "thermal_solver_3d.h"
#include <iostream>

void testCase1() {
    std::cout << "=== 测试用例1：立方体瞬态导热 ===" << std::endl;
    
    // 创建求解器：21x21x21网格，1m x 1m x 1m计算域
    ThermalSolver3D solver(21, 21, 21, 1.0, 1.0, 1.0);
    
    // 设置材料属性（铝）
    solver.setMaterialProperties(237.0, 2700.0, 900.0);
    
    // 设置时间参数
    solver.setTimeParameters(0.1, 100.0);  // 时间步长0.1s，总时间100s
    
    // 设置边界条件：底面固定温度100°C，其他面绝热
    solver.setBoundaryConditions(
        ThermalSolver3D::ADIABATIC, ThermalSolver3D::ADIABATIC,  // X方向
        ThermalSolver3D::ADIABATIC, ThermalSolver3D::ADIABATIC,  // Y方向  
        ThermalSolver3D::DIRICHLET, ThermalSolver3D::ADIABATIC,  // Z方向
        100.0, 0.0  // 底面100°C，顶面绝热
    );
    
    // 设置初始条件：均匀温度20°C
    solver.setInitialCondition(20.0);
    
    // 设置局部热源
    solver.setLocalHeatSource(8, 12, 8, 12, 15, 18, 100000.0);  // 100 kW/m³
    
    // 求解
    solver.solveExplicit();
}

void testCase2() {
    std::cout << "\n=== 测试用例2：高斯初始温度分布 ===" << std::endl;
    
    // 创建求解器
    ThermalSolver3D solver(31, 31, 31, 2.0, 2.0, 2.0);
    
    // 设置材料属性
    solver.setMaterialProperties(50.0, 2000.0, 1000.0);
    
    // 设置时间参数
    solver.setTimeParameters(0.05, 50.0);
    
    // 设置边界条件：所有边界绝热
    solver.setBoundaryConditions(
        ThermalSolver3D::ADIABATIC, ThermalSolver3D::ADIABATIC,
        ThermalSolver3D::ADIABATIC, ThermalSolver3D::ADIABATIC, 
        ThermalSolver3D::ADIABATIC, ThermalSolver3D::ADIABATIC
    );
    
    // 设置高斯初始条件：中心温度500°C
    solver.setInitialConditionGaussian(1.0, 1.0, 1.0, 500.0, 0.2);
    
    // 求解
    solver.solveExplicit();
}

int main() {
    std::cout << "三维温度场显式求解器 - Tecplot输出" << std::endl;
    std::cout << "=====================================" << std::endl;
    
    // 运行测试用例
    testCase1();
    testCase2();
    
    std::cout << "\n所有计算完成！结果已输出为Tecplot格式文件。" << std::endl;
    std::cout << "可以使用Tecplot或Paraview进行可视化。" << std::endl;
    
    return 0;
}

4. 编译和运行

CMakeLists.txt

cmake_minimum_required(VERSION 3.10)
project(ThermalSolver3D)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

add_executable(thermal_solver_3d 
    main.cpp 
    thermal_solver_3d.cpp
)

# 设置编译器优化
if(CMAKE_BUILD_TYPE STREQUAL "Release")
    target_compile_options(thermal_solver_3d PRIVATE -O3)
endif()

编译运行

mkdir build
cd build
cmake ..
make
./thermal_solver_3d

参考代码 显式求解三维温度场并输出Tecplot文件 www.youwenfan.com/contentcnk/70480.html

5. 主要特性

1. 显式时间积分：简单易实现，适合瞬态问题
2. 多种边界条件：Dirichlet、Neumann、绝热边界
3. 热源设置：支持均匀热源和局部热源
4. 稳定性检查：自动检查傅里叶数稳定性条件
5. Tecplot输出：标准格式，支持多种后处理软件

6. 输出文件示例

生成的Tecplot文件可以直接用Tecplot、Paraview等软件打开进行三维可视化，显示温度场的时空演化。