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

news/2025/10/30 15:47:56/文章来源:https://www.cnblogs.com/m877087643/p/19177264

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°Csolver.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°Csolver.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等软件打开进行三维可视化,显示温度场的时空演化。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/950893.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

哈希重温

字符串哈希(进制哈希)把字符串通过哈希函数映射成一个数,类似进制的方法处理即可。 下文我们约束字符串为 $str$,长度分别为 $n$,以 $1$ 为起始下标。 类比成 $base$ 进制的理解,有: $$hash(str)=\sum_{i=1}^{n…

.NET6 Web程序部署在IIS上

1、应用程序在IDE中进行发布 2、需要在Window 服务器上开启IIS 3、ASP.NET Core 应用针对IIS部署依赖于一个IIS针对ASP.NET Core的扩展模块。所以需要下载ASP.NET Core 运行时的Hosting Bundle。ASP.NET Core ModuleAS…

20251030模拟赛

T1 小模拟题,场切了。

HarmonyOS自动化测试与持续集成实战指南

1. 自动化测试概述与工具链介绍 自动化测试是HarmonyOS应用开发流程中保证质量的关键环节。随着HarmonyOS生态的快速发展,应用功能日益复杂,分布式特性、多设备适配需求以及快速迭代的开发模式,使得传统手动测试无法…

HarmonyOS应用性能调优与内存管理实战

1. 性能优化概述与指标体系 性能优化是HarmonyOS应用开发中不可或缺的一环,它直接影响用户体验和应用稳定性。一个高性能应用应具备快速启动、流畅交互和低资源消耗等特点。在HarmonyOS应用性能评估中,我们需要关注几…

GEO优化源头厂家怎么选?这篇干货帮你摸透门道!

最近不少老板都在问:“GEO优化到底是个啥?源头厂家哪家强?” 别急,咱们今天就用大白话,掰开揉碎了聊一聊,顺便给大家推荐一个实力派选手——讯灵AI(GEO+Agent)双引擎系统。 一、GEO是啥?为啥它成了香饽饽? 简单…

HarmonyOS大型项目架构与模块化开发指南

1. HarmonyOS大型项目架构设计 HarmonyOS大型项目开发需要采用分层架构和模块化设计,确保代码的可维护性、可扩展性和团队协作效率。合理的架构设计是项目成功的基础。 1.1 分层架构设计原则 HarmonyOS推荐采用四层架…

鸿蒙NDK开发实战指南:从ArkTS到C/C++的高性能桥梁

1. NDK概述与核心价值 HarmonyOS NDK(Native Development Kit)是HarmonyOS SDK中提供的Native API、编译脚本和编译工具链的集合,它让开发者能够使用C或C++语言实现应用的关键功能模块。NDK主要覆盖了HarmonyOS的基…

HarmonyOS后台任务管理:短时任务与长驻任务实战

一、HarmonyOS后台任务概述 HarmonyOS的后台任务管理旨在平衡任务执行需求与系统资源消耗,提供多种后台任务类型以满足不同场景的需求。后台任务主要分为短时任务和长驻任务两大类,每种类型都有特定的使用场景和限制…

GEO 源头厂家独家王炸:南方网通讯灵 AI 业内首创“3+4+3” 智能生态营销体系,领爆AI搜索新浪潮

近日,深圳彻底沸腾!一场聚焦GEOAI 搜索时代的“流量破局盛会”高能集结 —— 这不是普通行业会,而是专门为 “苦流量枯竭久矣” 的企业主量身定制的增长急救场! 现场汇聚了200+企业掌舵人,他们带着各自行业的实战…

HarmonyOS分布式硬件共享:调用手机摄像头的手表应用

一、分布式硬件共享概述 分布式硬件共享是HarmonyOS的核心能力之一,它基于分布式硬件池理念,将网络中多个物理设备的硬件资源进行统一虚拟化管理。这意味着应用程序可以按需调用任意可信设备的硬件能力,打破传统单设…

HarmonyOS应用配置文件与资源组织深度解析

一、应用程序包结构概述 HarmonyOS应用以APP Pack形式发布,它由一个或多个HAP以及描述每个HAP属性的pack.info组成。HAP是Ability的部署包,可分为entry和feature两种类型。 Entry类型的HAP是应用的主模块,一个应用中…

OpenHarmony内核基础:LiteOS-M内核与POSIX/CMSIS接口

1. OpenHarmony内核架构概述 OpenHarmony采用多内核设计理念,根据设备资源能力匹配不同的内核形态,为各种物联网设备提供精准化的系统支持。这种设计使得OpenHarmony能够灵活适应从低端资源受限设备到高端智能设备的…

GEO源头厂家怎么选?看这3点:研发实力、产品核心、交付标准,精准避坑不选错!

当AI搜索的浪潮席卷而来,无数企业将“GEO优化”视为破局的关键。然而,面对市场上纷繁复杂的技术服务商,一个核心问题浮出水面:GEO源头厂家怎么选?是选择仅有概念包装的中间商,还是拥有核心技术与深厚底蕴的研发源…

2025年防腐蚀地坪生产厂家权威推荐榜单:聚脲防腐地坪/化工厂防腐工程/三布六油防腐蚀地坪源头厂家精选

防腐蚀地坪作为工业设施中的关键组成部分,广泛应用于化工、食品、医药、环保等腐蚀性环境。根据行业数据,全球防腐蚀地坪市场规模预计在2025年将达到120亿美元,年复合增长率稳定在6.5%。在中国市场,防腐蚀地坪需求…

2025年可靠的水电镀表面处理厂家推荐及选购参考榜

2025年可靠的水电镀表面处理厂家推荐及选购参考榜 行业概述 水电镀表面处理作为现代制造业中不可或缺的工艺环节,广泛应用于金属、塑料、玻璃、陶瓷等材料的表面处理,以提高产品的耐腐蚀性、耐磨性、美观度及功能性…

机器学习中,验证阶段为什么还要返回损失?

为什么验证阶段还要返回损失? 在验证阶段返回损失(val_loss)是模型训练中评估性能、指导训练的核心逻辑,主要有以下几方面原因:评估模型泛化能力 训练阶段的损失(train_loss)只能反映模型对训练数据的拟合程度,…

JYU-ACM算法协会每日一题题解(每日刷新)

P8754 [蓝桥杯 2021 省 AB2] 完全平方数 点击跳转 P2818 天使的起誓 点击跳转 P5707 【深基2.例12】上学迟到 点击跳转

revit api previewcontrol wpf预览窗口

revit api previewcontrol wpf预览窗口族库管理插件 不直接打开文件,预览族文件的指定view, 进一步可以使用using 进行资源释放操作,这里没有写出来using Autodesk.Revit.Attributes; using Autodesk.Revit.DB; u…

2025年VOC废气治理RTO蓄热焚烧炉供应商权威推荐榜单:TO废气焚烧炉/氟化工废液废气焚烧炉 /含氯含氟废气处理厂家精选

随着环保政策持续收紧,高效可靠的VOC废气治理设备正成为工业企业实现达标排放和绿色发展的关键。 在“十四五”期间VOCs排放总量需较2020年下降10%的政策硬约束下,蓄热式热力氧化技术(RTO)凭借其高热回收效率和卓越…