WO2024041642A1 - 时域热传导仿真方法及存储介质 - Google Patents

时域热传导仿真方法及存储介质 Download PDF

Info

Publication number
WO2024041642A1
WO2024041642A1 PCT/CN2023/114983 CN2023114983W WO2024041642A1 WO 2024041642 A1 WO2024041642 A1 WO 2024041642A1 CN 2023114983 W CN2023114983 W CN 2023114983W WO 2024041642 A1 WO2024041642 A1 WO 2024041642A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
time domain
heat conduction
time
finite element
Prior art date
Application number
PCT/CN2023/114983
Other languages
English (en)
French (fr)
Inventor
蒲菠
何秋森
范峻
Original Assignee
宁波德图科技有限公司
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 宁波德图科技有限公司 filed Critical 宁波德图科技有限公司
Publication of WO2024041642A1 publication Critical patent/WO2024041642A1/zh

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Definitions

  • the present invention relates to the technical field of thermal simulation, and specifically relates to a time domain heat conduction simulation method based on high-order tetrahedral units, and a computer-readable storage medium for executing steps in the time domain heat conduction simulation method.
  • Kane Yee first proposed an electromagnetic time domain simulation method based on the Finite Difference Time Domain (FDTD) method.
  • FDTD Finite Difference Time Domain
  • Jin-fa Lee et al. proposed and introduced in detail the electromagnetic time domain numerical simulation method based on the Finite Element Time Domain (FETD) method.
  • FETD Finite Element Time Domain
  • the present invention provides a time-domain heat conduction simulation method based on high-order tetrahedral units.
  • the time domain heat conduction simulation method provided by the present invention includes the following steps:
  • S1 use finite element meshing tools for geometric modeling and meshing to read geometric information
  • S3 through the number of each unit, node, and edge, correspond the elements in the unit matrix to the elements of the global matrix, and use the unit matrix through traversal and superposition to obtain the global matrix;
  • the time domain heat conduction simulation method also includes:
  • step S1 geometric modeling and meshing are performed through a finite element meshing tool (such as Trelis), and read through a third-party simulation software (such as Matlab) platform Get geometric information.
  • a finite element meshing tool such as Trelis
  • a third-party simulation software such as Matlab
  • the materials in each region Relevant parameters include: thermal conductivity, mass density, specific heat capacity, boundary conditions, and domain point probes required for observation data.
  • T is the temperature
  • C is the material specific heat capacity
  • t is time
  • ⁇ t is the simulation time interval
  • [C T ] is the damping matrix
  • [K T ] is the stiffness matrix
  • ⁇ p T ⁇ is the excitation source
  • ⁇ T ⁇ is the discrete form
  • ⁇ T ⁇ t is the temperature field scalar in discrete form at time t
  • ⁇ T ⁇ t- ⁇ t is the temperature field scalar in discrete form at time t- ⁇ t.
  • the unit matrix of each unit includes a mass matrix, a damping matrix and a stiffness matrix.
  • step S4 of the time domain heat conduction simulation method the LUPQ matrix decomposition method is used for iterative operations.
  • step S5 of the time domain heat conduction simulation method third-party simulation software (such as COMSOL Multiphysics) is used as a verification tool by constructing a calculation example.
  • third-party simulation software such as COMSOL Multiphysics
  • the present invention also provides a computer-readable storage medium for executing steps in any one of the above-mentioned time domain heat conduction simulation methods.
  • This invention uses finite element meshing tools for geometric modeling and meshing to read geometric information; constructs basic equations of thermal time domain finite elements based on high-order tetrahedral units to calculate and obtain the global matrix; and then calculates the global matrix through thermal time domain finite elements.
  • the basic equation of the thermal time domain finite element is iteratively calculated for each time step; the verification tool is used to verify the results of the basic equation of the thermal time domain finite element to obtain accurate time domain heat conduction simulation results.
  • Reliable and efficient time domain heat conduction equation simulation method is suitable for high The comprehensive design of performance electronic products and the research and verification of scientific theoretical issues are of great significance.
  • Figure 1 is a schematic diagram of meshing the object to be analyzed
  • Figure 2 is a schematic diagram comparing the FETD thermal solution and the COMSOL Multiphysics results
  • Figure 3 is a schematic diagram of the relative errors comparing FETD thermal solution and COMSOL Multiphysics
  • Figure 4 Schematic diagram of a body heat source equipped with a heat sink.
  • the physical field involved in the present invention is a thermal field (temperature distribution), which is described accordingly.
  • the equation described above is the heat conduction equation.
  • the heat conduction equation describes the phenomenon of heat transfer in solid matter.
  • T is the temperature
  • t is the time
  • ⁇ is the mass density
  • c is the material specific heat capacity
  • ⁇ is the thermal conductivity coefficient
  • Q is the thermal power density
  • Laplacian operator is the temperature gradient.
  • T is the temperature
  • T f is the fluid environment temperature
  • ⁇ is the thermal conductivity coefficient
  • ⁇ is the boundary surface normal vector
  • Laplacian operator is the Laplacian operator
  • h is the convection heat transfer coefficient.
  • Numerical methods use computers to solve approximate solutions to mathematical problems. Many related numerical algorithms have been proposed for different problems. Algorithms commonly used to solve systems of partial differential equations include the method of moments, the finite difference method and the finite element method used in the present invention.
  • the present invention provides a time domain heat conduction simulation method, which includes the following steps:
  • S1 use finite element meshing tools for geometric modeling and meshing to read geometric information
  • Geometric modeling is performed in Cubit Trelis by writing .jou scripts in combination with the GUI.
  • Trelis is a dedicated finite element meshing tool that has built-in geometric modeling tools and allows users to set label information (such as material number, geometry number, etc.) for the constructed geometric model.
  • Cubit Trelis's geometric modeling supports geometric structures from one to three dimensions. Following a professional tool-based approach, Cubit Trelis offers a variety of different meshing techniques, collaboration infrastructure and algorithms. Cubit Trelis supports .jou scripts or python scripts for geometric modeling. When generating calculation examples in batches, you can use script modeling to increase speed.
  • Meshing is implemented in Cubit Trelis by writing .jou scripts combined with GUI. Export as .g file.
  • Cubit Trelis' meshing supports hexahedral elements, tetrahedral elements and other common element types.
  • Basis Functions In the finite element method, basis functions are local. The finite element method divides large geometric objects into smaller and simpler units through meshing. Based on these units, basis functions are constructed.
  • Import geometric information Convert the .g file to a text file, and then write a script in the Matlab platform to read the geometric information contained in the file.
  • post-processing can be performed as needed.
  • the general approach is to use a separate .m pair of parameter information, that is, the thermal conductivity and mass density of the material in each area in the solid heat transfer analysis. , specific heat capacity and boundary conditions and other parameters are set.
  • the domain point probes required for observation data need to be set in this step.
  • Finite element elements can be classified according to the geometry of the base element.
  • common shapes include: cuboid, tetrahedron, pentahedron, curved hexahedron, etc.
  • the cuboid element is a relatively simple element.
  • the subdivision algorithm is relatively simple, and relatively simple objects can even be divided manually, making it more direct and convenient to obtain geometric information.
  • the cuboid unit has low performance when simulating geometric structures. Therefore, simple cuboids are rarely used directly as the basic unit of finite elements in industry.
  • the tetrahedral unit unit is more flexible in simulating geometric shapes and produces smaller errors. Therefore, it is widely used in thermal, electromagnetic, and mechanical finite element analysis.
  • the pentahedral unit Compared with the tetrahedron, the pentahedral unit has no obvious improvement in accuracy, but each unit has more vertices and edges, which increases the amount of calculation to a certain extent.
  • the scope of application of pentahedral units is not as wide as that of tetrahedral units, and it is difficult to find suitable references, which to a certain extent increases the difficulty of writing programs to implement finite element guidelines.
  • the geometric simulation capability of the curved hexahedron is stronger than that of the tetrahedral unit, but the curved hexahedral unit cannot be simply described by vertices, which brings great trouble to the processing of geometric information.
  • additional auxiliary points are needed when defining their shape to ensure that the element accurately fits the local geometry of the simulation object.
  • the process of meshing and spatial discretization is complicated and difficult to implement. The degree is higher.
  • tetrahedral units with high performance and feasibility are selected as basic units.
  • the basis functions of finite elements are distinguished by function variables and can be divided into two categories: vector basis functions and scalar basis functions.
  • the field quantity of thermal analysis is mainly temperature, which is a scalar quantity, so when implementing thermal finite elements, scalar basis functions are used for analysis.
  • the basis functions can be defined according to the vertices or edges of the finite element unit. For vector basis functions, it is more convenient to use edge definitions; for scalar basis functions, it is more convenient to use vertex definitions. Therefore, thermal finite elements use point-based scalar basis functions.
  • a first-order tetrahedron has four nodes, one at each vertex.
  • the second-order tetrahedron unit adds one node to each of the six edges of the tetrahedron.
  • the third-order tetrahedron unit adds six new nodes to the six edges of the second-order tetrahedron, and adds one node to each of the six edges of the tetrahedron.
  • a node is also added to the four faces of .
  • Higher-order elements will contain more node information, so by using higher-order elements to divide the solution area, the finite element calculation accuracy will be greatly improved. Due to the consideration of computational difficulty and complexity, the present invention selects a second-order tetrahedral structure to mesh the solution area.
  • the next step is to select a set of trial functions with certain properties as interpolation functions for the unknown solutions in the element.
  • the selection of the interpolation function will be directly related to the accuracy of the finite element calculation results. For a unit e containing n nodes, no matter what order the unit is, each node will correspond to an interpolation function. Although the interpolation function of high-order units will be more complex than the interpolation formula of low-order units, its calculation accuracy will also be greatly improved.
  • the interpolation function is divided into two types: corner nodes and inner-edge nodes:
  • N 5 4L 1 L 2
  • N 6 4L 1 L 3
  • Ni is the node basis function
  • Li is the coordinate representation in the tetrahedral element coordinate system.
  • the spatial discretization of geometry and physical fields in the finite element method is based on many tiny units, so it is called the finite element method.
  • the finite element method In the early days, the theoretical development of finite element was based on the variational method. When the finite element method was first used for calculation and simulation, the simulations used mostly used the variational method to derive its basic equations. Therefore, its application scenarios in the field of physics are mostly scenarios described by the Laplace equation and the Poisson equation. The finite element method is therefore often used for problems closely related to functional extreme value problems. In subsequent research, it was discovered that Galerkin’s Method (Galerkin’s Method) in the Weighted Residual Method (WRM) can also be used to obtain the same form of basic finite element equations.
  • Galerkin’s Method Galerkin’s Method
  • WRM Weighted Residual Method
  • the finite element method can be used to simulate a variety of physical fields including thermal, optics, electromagnetics, fluid mechanics, and solid mechanics. Therefore, it has been widely used in various fields of industry and academia. Because it is suitable for most physical scenarios, finite element is also a more suitable algorithm for simulating multi-physics coupling problems.
  • Finite elements realize numerical calculations by performing specific spatial discretizations in spatial dimensions.
  • the specific process includes three steps: geometric meshing, spatial discretization based on meshing and selected basis functions, and final calculation (iterative calculation for the time-domain finite element method).
  • the geometric gridding provides basic geometric information for space discretization, and space discretization converts the original complex problem into a linear equation system for final calculation.
  • Time-domain finite element is the time-domain version of the finite element method. It performs iterative calculations through a certain time-stepping method to obtain the change of field quantities over time in the entire observation time domain.
  • Common time-stepping methods include before/after Difference method and Runge-Kutta method.
  • T is the temperature
  • [C T ] is the damping matrix
  • [K T ] is the stiffness matrix
  • ⁇ T ⁇ is the temperature field scalar in discrete form
  • ⁇ p T ⁇ is the excitation source.
  • N is the basis function used in the finite element
  • ⁇ is the mass density
  • c is the material specific heat capacity
  • ⁇ is the thermal conductivity coefficient
  • Q is the thermal power density
  • h is the convection heat transfer coefficient.
  • the thermal time discretization method uses the backward difference method (Backward Difference Method, BDM) [19].
  • BDM Backward Difference Method
  • T is temperature
  • t time
  • T(t) is the temperature at time t
  • T(t- ⁇ t) is the temperature at time t- ⁇ t
  • O( ⁇ t) is the error term.
  • T is the temperature
  • C is the material specific heat capacity
  • t time
  • ⁇ t is the simulation time interval
  • [C T ] is the damping matrix
  • [K T ] is the stiffness matrix
  • ⁇ p T ⁇ is the excitation source
  • ⁇ T ⁇ is the discrete
  • the temperature field scalar in the form of, ⁇ T ⁇ t is the temperature field scalar in the discrete form at time t, ⁇ T ⁇ t- ⁇ t is the temperature field scalar in the discrete form at time t- ⁇ t.
  • each element in the element matrix is calculated.
  • Thermal uses linear point-based scalar basis functions based on tetrahedral elements. It is crucial to obtain several matrices in the equation (mass matrix, damping matrix, stiffness matrix). Since the degrees of freedom of the limited basis functions are assigned to the edges or points of each unit, the physical field in a unit is determined only by the physical field values and basis functions of the nodes and edges of this unit, and is different from the basis functions of other units. Regardless, the basis function in this method is generally called the local basis function. Therefore, when constructing the matrix in the finite element equation, the element matrix (Element matrix) of each unit is generally calculated according to the formula for calculating matrix elements, and finally the matrices of each unit are assembled to obtain the final required global matrix (Global matrix). ).
  • S3 through the number of each unit, node, and edge, correspond the elements in the unit matrix to the elements of the global matrix, and use the unit matrix through traversal and superposition to obtain the global matrix;
  • the basis function used by the finite element method is a local basis function, this basis function only has a non-zero field value in units that contain points or edges that define its degrees of freedom, so the basis functions in adjacent units are not It is generally orthogonal.
  • the matrix obtained by using this finite element basis function is generally a sparse matrix.
  • matrix decomposition methods such as LU decomposition can be used to improve the computing performance of the program (the efficiency of using matrix decomposition to solve large linear equations is much higher than directly inverting large matrices).
  • LUPQ matrix decomposition method is used.
  • a s is in the final formula of the thermal time domain finite element term, x is ⁇ T ⁇ t item, b is item.
  • the second embodiment of the present invention is formed based on the above-mentioned first embodiment, and the same parts will not be described again. The following steps are added on the basis of the above-mentioned first embodiment;
  • COMSOL Multiphysics As an example, the reference software used for verification and final verification is COMSOL Multiphysics.
  • COMSOL Multiphysics was originally developed as an application on the Matlab platform. It is widely used in industry and academia to numerically solve physical problems described by partial differential equations through certain mathematical methods.
  • COMSOL Multiphysics supports a wide range of algorithms, and there are many optimization algorithms for different scenarios, but these algorithms are essentially finite element methods.
  • solver and multi-physics simulation software As a cross-platform finite element analysis, solver and multi-physics simulation software.
  • COMSOL Multiphysics allows modeling and constructing coupled systems of partial differential equations through a traditional physics-based user interface. It is a software with a high degree of freedom in use.
  • COMSOL provides IDE and unified workflow for simulation research in different fields such as mechanics, fluids, electrical, radio frequency, acoustics, electrochemical reactions, plasma, and optics.
  • the present invention constructed several calculation examples and compared them with the results of COMSOL Multiphysics.
  • the first calculation example is copper block cooling: a 0.1m ⁇ 0.1m ⁇ 0.1m copper block, its thermal parameters: thermal conductivity is 377W/(m ⁇ K), density is 8900kg/m 3 , constant pressure heat capacity is 386J /(kg ⁇ K), the boundary condition of the solid is the convection boundary condition, the temperature of the convection fluid is 300K, and the heat transfer coefficient is 500.
  • the results obtained by the self-written FETD thermal solver are basically consistent with the results in COMSOL Multiphysics, and the error is within an acceptable range.
  • the second calculation example is a body heat source equipped with a heat sink, as shown in Figure 4.
  • the entire object is in convective boundary conditions. Since there are many relevant parameters for this model design, the specific parameters will not be described in detail. The specific parameters will be given in the appendix.
  • the present invention provides a computer-readable storage medium for executing steps in the time-domain heat conduction simulation method described in any one of the first embodiment or the second embodiment.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种时域热传导仿真方法，包括：利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；设置固体传热学分析中各个区域材料相关参数；构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；通过热学时域有限元的基本方程对每一个时间步进行迭代运算。本发明能获得准确的时域热传导仿真结果。

Description

时域热传导仿真方法及存储介质
本申请要求于2022年08月26日提交中国专利局、申请号为202211032489.1、发明名称为“时域热传导仿真方法及存储介质”的中国专利申请的优先权，其全部内容通过引用结合在本申请中。
技术领域
本发明涉及热学仿真技术领域，具体涉及一种基于高阶四面体单元的时域热传导仿真方法，以及一种用于执行所述时域热传导仿真方法中步骤的计算机可读存储介质。
背景技术
20世纪40年代早期，A.Hrennikoff和R.Courant针对建筑结构和飞行器结构的设计的力学分析提出发展了现在称为有限元法(Finite Element Method)的数值方法。
1965年，Winslow将FEM应用于电气工程分析中。1966年，Kane Yee首先提出了基于时域有限差分法(Finite Difference Time Domain，FDTD)的电磁学时域仿真方法。
1969年，Silvester首先将FEM运用到高频领域，应用于微波波导设计。
1997年，Jin-fa Lee等人提出并详细介绍了基于时域有限元法(Finite Element Time Domain，FETD)的电磁学时域数值仿真方法。
2016年，Jian-Ming Jin等人基于FEM，FETD实现了电热耦合仿真，利用FETI提升计算效率，应用于射频器件的热分析。
2019年，Zhang Hao-Xuan等人采用并行计算技术优化射频元器件多物理仿真的效率。
2020年，Yan Su，Jian-Ming Jin等人发表“An enhanced transient solver with dynamic p‐adaptation and multirate time integration for electromagnetic and multiphysics simulations”。基于时域离散伽辽金法(Discontinuous Galerkin Time Domain，DGTD)开发了多物理(电磁学与等离子体)的仿真解算器，针对多尺度问题(Multi-scale problem)利用动态p优化适应和多时间步的方法提升计算效率。
关于有限元的算法研究已经比较完善，但是，对于基于高阶四面体单元的热学时域有限元求解方法还比较少。部分研究在进行电热多物理仿真时，只采用了单项耦合的模型，计算了电磁场的热效应而忽略了温度对于电磁场的影响。目前，对各个物理场都采用时域分析(即全时域)的电热双向耦合仿真方法的研究是很少的。因此，基于高阶四面体单元的热学时域有限元求解方法对于多物理问题有很大的作用。
发明内容
在发明内容部分中引入了一系列简化形式的概念，该简化形式的概念均为本领域现有技术简化，这将在具体实施方式部分中进一步详细说明。本发明的发明内容部分并不意味着要试图限定出所要求保护的技术方案的关键特征和必要技术特征，更不意味着试图确定所要求保护的技术方案的保护范围。
为了克服上述现有技术的不足，本发明提供了一种基于高阶四面体单元的时域热传导仿真方法。
为解决上述技术问题，本发明提供的时域热传导仿真方法，包括以下步骤：
S1，利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；
S2，设置固体传热学分析中各个区域材料相关参数；
S3，构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；
S3，通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；
S4，通过热学时域有限元的基本方程对每一个时间步进行迭代运算。
可选择的，所述的时域热传导仿真方法，还包括：
S5，通过验证工具对热学时域有限元的基本方程结果进行验证。
可选择的，所述的时域热传导仿真方法，步骤S1中，通过有限元网格剖分工具(例如Trelis)进行几何建模和网格剖分，通过第三方仿真软件(例如Matlab)平台读取几何信息。
可选择的，所述的时域热传导仿真方法，步骤S2中，各个区域材料 相关参数包括：热导率，质量密度，比热容、边界条件以及观测数据所需的域点探针。
可选择的，所述的时域热传导仿真方法，热学时域有限元的基本方程表达如下；
T为温度，C为材料比热容，t为时间，Δt为仿真时间间隔，[CT]为阻尼矩阵，[KT]为刚度矩阵，{pT}为激励源，{T}为离散形式的温度场标量，{T}t为t时刻离散形式的温度场标量，{T}t-Δt为t-Δt时刻离散形式的温度场标量。
可选择的，所述的时域热传导仿真方法，步骤S3中，每个单元的单元矩阵包括质量矩阵、阻尼矩阵和刚度矩阵。
可选择的，所述的时域热传导仿真方法，步骤S4中，进行迭代运算采用LUPQ矩阵分解法。
可选择的，所述的时域热传导仿真方法，步骤S5中，通过构建算例采用第三方仿真软件(例如COMSOL Multiphysics)作为验证工具。
本发明还提供了一种用于执行上述任意一项所述的时域热传导仿真方法中步骤的计算机可读存储介质。
传统的电子器件设计环节中，仅对与电子器件工作相关的电磁场进行分析，经常忽略了力学热学等其他物理场对于器件性能的影响。21世纪以来，电子产品向着大功率，高性能，小型化的方向发展，电子元件、设备的发热问题越来越不可忽视。因此温度场的分析于设计在工业中愈发重要。任何让工程师更好地分析不同物理场之间的联系，成为了高性能产品工业设计中一个十分重要的问题。科学研究上，随着跨领域交叉学科的发展，对于热学时域有限元研究的需求也逐渐增大。
本发明利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；基于高阶四面体单元构建热学时域有限元的基本方程计算获得全局矩阵；再通过热学时域有限元的基本方程对每一个时间步进行迭代运算；通过验证工具对热学时域有限元的基本方程结果进行验证进而能获得准确的时域热传导仿真结果。可靠高效的时域热传导方程仿真方法对于高 性能电子产品综合设计与科学理论问题的研究验证有着十分重要的意义。
说明书附图
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1是对待分析对象进行网格剖分示意图；
图2是FETD热学求解与COMSOL Multiphysics结果对比示意图；
图3是FETD热学求解与COMSOL Multiphysics对比相对误差示意图；
图4搭载散热片的体热源示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明的描述中，需要说明的是，术语“上”、“下”、“水平”、“内”、“外”、“顶”、“底”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“形成”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
本发明所涉及的物理场为热学场(温度分布)，对应的对其进行描 述的方程为热传导方程。
热传导方程描述了固体物质中的热传递现象。
其中T为温度，t为时间，ρ为质量密度，c为材料比热容，λ为热传导系数，Q为热功率体密度，为拉普拉斯算符，为温度梯度。
在热分析中，采用对流边界条件(ConvictionBoundary)：
其中T为温度，Tf为流体环境温度，λ为热传导系数，为边界面法向量，为拉普拉斯算符，为温度梯度，h为对流换热系数。
数值方法是采用计算机求解数学问题近似解的问题。对于不同的问题，已经提出了很多相关的数值算法。常用于求解偏微分方程组的算法有矩量法，有限差分法以及本发明使用的有限元法。
本发明提供一种时域热传导仿真方法，包括以下步骤：
S1，利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；
示例性的，几何建模：通过编写.jou脚本结合GUI在Cubit Trelis中进行几何建模。Trelis是一专用的有限元网格剖分工具，其中内置了几何建模工具，并允许用户对构建的几何模型设置标签信息(如材料编号，几何体编号等)。Cubit Trelis的几何建模支持从一维到三维的几何结构。遵循基于专业工具的方法，Cubit Trelis提供了各种不同的网格划分技术，协作基础架构和算法。Cubit Trelis支持.jou脚本或python脚本进行几何建模，在批量生成算例时，可以使用脚本建模提升速度。
网格剖分：通过编写.jou脚本结合GUI在Cubit Trelis中实现网格剖分。以.g file文件导出。处理三维模型时Cubit Trelis的网格剖分支持六面体单元，四面体单元等常用单元类型，在批量生成算例时，可以使用脚本建模提升速度。基函数在有限元法中，基函数是局部的。有限元法通过网格剖分将大型的几何对象剖分为更小，更简单的单元，基于这些单元，基函数得以构建。
导入几何信息：将.g file转换为文本文件，再在Matlab平台中编写脚本读取文件中包含的几何信息。
S2，设置固体传热学分析中各个区域材料相关参数；
示例性的，将读取的信息导入Matlab workspace后，可以根据需要进行后处理，一般的做法采用单独的.m对参数信息，即固体传热学分析中各个区域材料的热导率，质量密度，比热容以及边界条件等参数进行设置。此外，对于观测数据所需的域点探针，需要在此步骤设置。
S3，构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；
有限元单元可根据基础单元的几何形状分类。针对三维有限元仿真，常见的形状有：长方体，四面体，五面体，曲面六面体等。
长方体单元是比较简单的单元。选用长方体单元时，剖分算法较为简单，对于比较简单的对象甚至可以人工进行划分，获取几何信息时较为直接方便。但是受到阶梯效应的影响，长方体单元在模拟几何结构时性能较低，因此在工业中比较少直接使用简单长方体作为有限元的基本单元。
四面体单元单元单元相对于长方体单元而言，在模拟几何外形时更加灵活，产生的误差较小，因此广泛的应用于热学，电磁学，力学有限元分析中。
五面体单元相较于四面体而言没有明显的精度上的提升，而每个单元的顶点数和棱边数更多，一定程度上提高了计算量。此外，五面体单元的适用范围不如四面体单元广泛，较难找到合适的参考文献，这一定程度上增大了编写程序实现有限元方针的难度。
曲面六面体的几何模拟能力比四面体单元更强，但是曲面六面体单元不能通过顶点简单的描述，这给几何信息的处理带来很大的麻烦。对于具有弯曲边界面的单元，在定义其形状时，需要额外的辅助点来保证单元足够准确地拟合了仿真对象的局部几何形状。选用曲面六面体单元实现有限元，在进行网格剖分以及空间离散化时过程比较复杂，实现难 度较高。
本发明选取性能与可行性皆较高的四面体单元作为基本单元。
有限元的基函数，以函数变量进行区分，可分为矢量基函数与标量基函数两大类。热学分析的场量主要为温度，属于标量，所以实现热学有限元时采用标量基函数进行分析。基函数可以根据有限元单元的顶点或棱边进行定义，对于矢量基函数，采用棱边定义是比较方便的；而对于标量基函数，采用顶点定义比较方便。因此，热学有限元采用点基标量基函数。
为了求解得到更加精确的结果，可以应用更高次的单元对求解区域进行剖分。图中显示了各阶四面体单元。一阶四面体有四个结点，分别在各个顶点处。二阶四面体单元在四面体的六条棱边内各增加了一个结点，三阶四面体在二阶四面体的基础上又在六条棱边内新增加了六个结点，并在四面体的四个面内也增加了一个结点。高阶单元将会包含更多的结点信息，因此采用高阶单元剖分求解区域，有限元计算精度会极大改善。出于计算难度和复杂度的考虑，本发明选用二阶四面体结构来对求解区域进行网格剖分。
使用高阶四面体单元对求解域进行网格剖分后，接下来要选择一组具有某些性质的试探函数作为单元中未知解的插值函数。插值函数的选取将直接关系到有限元计算结果的准确性。对于含有n个结点的单元e中，无论该单元是多少阶的，每一个结点将对应一个插值函数。虽然高阶单元的插值函数会比低阶单元的插值公式更复杂，但它的计算精度也会大大提高。
对于二阶四面体单元，插值函数分为角结点和棱内结点两种：
角结点：
Ni＝(2Li-1)Li (i＝1，2，3，4）
棱内结点：
N5＝4L1L2
N6＝4L1L3
其中Ni为结点基函数，Li为四面体单元坐标系中的坐标表示。
有限元法中几何和物理场的空间离散化是基于许多微小单元的，故被称作有限元法。在早期，有限元的理论发展是基于变分方法的，在最开始使用有限元法计算仿真时，使用的仿真多用变分法推导其基本方程。因此它在物理领域的应用场景多为由Laplace方程和Poisson方程描述的场景。有限元法也因此常被用于与功能极值问题密切相关的问题。在之后的研究中，发现可以使用加权余量法(Weighted Residual Method，WRM)中的伽辽金法(Galerkin’s Method)也可以得出相同形式的有限元基本方程。这一理论上的发展很大程度拓宽了有限元法的应用场景。理论上，任何偏微分方程组描述的物理场都可以用有限元法进行仿真计算。因此，近年来越来越多关于有限元法的研究都改使用伽辽金法作为推导有限元法方程的基础，并把有限元法归为加权余量法的一个分支。在本毕业论文中的理论推导部分，也采用伽辽金法推导。
有限元法能用于仿真包括热学、光学、电磁学、流体力学、固体力学在内的等多种物理场，因此在工业界和学术界的各个领域得到了广泛的应用。因为适用于大多数物理场景，有限元也是比较适合用于仿真多物理场耦合问题的算法。
有限元是通过在空间尺寸上进行特定的空间离散化来实现数值计算的。具体的过程包括几何网格剖分，基于剖分的与选定基函数的空间离散化与最终的计算(对于时域有限元法而言是迭代计算)三个步骤。其中几何网格剖分为空间离散化提供了基础的几何信息，而空间离散化将原本复杂的问题转换为一个线性方程组以便最终的计算。
时域有限元则是有限元法的时域版本，通过某种时间步进方法进行迭代计算，从而得到整个观测时间域上场量随时间变化情况的方法，常见的时间步进方法有前/后向差分法和龙格-库塔法。
之前已经给出，热传导分析的域内控制方程为：
之前已经给出，热传导分析的边界内控制方程为：
根据之前给出的加权余量法，可以推导出热传导分析的矩阵形式有限元方程，这里直接给出：
上式中，T为温度，[CT]为阻尼矩阵，[KT]为刚度矩阵，{T}为离散形式的温度场标量，为温度对时间的微分，为离散形式的温度场标量对时间的微分，{pT}为激励源。上式中的矩阵以及列向量中的元素可由以下公式给出：
[CT]ij＝∫VρcNiNjdV

{pT}i＝∫VQNidV+∫shNiTadS
式中N是有限元中使用的基函数，ρ为质量密度，c为材料比热容，λ为热传导系数，Q为热功率体密度，为拉普拉斯算符，h为对流换热系数。
由于温度是标量，在热学分析中这些基函数采用点基(Nodal based)标量基函数。热学的时间离散方法采用后向差分法(Backward Difference Method，BDM)[19]。针对温度-时间的变化关系，后向差分法方程如下:
其中，T为温度，t为时间，T(t)为t时刻的温度，T(t-Δt)为t-Δt时刻的温度，为温度对时间的微分，Δt为时间间隔，O(Δt)为误差项。
无论dt的值如何选取，采用后向差分法构建的系统是无条件稳定的。将后向差分法带入原式，可以得到热学时域有限元的最终公式：
其中，T为温度，C为材料比热容，t为时间，Δt为仿真时间间隔，[CT]为阻尼矩阵，[KT]为刚度矩阵，{pT}为激励源，{T}为离散形式的温度场标量，{T}t为t时刻离散形式的温度场标量，{T}t-Δt为t-Δt时刻离散形式的温度场标量。
计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数。获取方程中的几个矩阵(质量矩阵，阻尼矩阵，刚度矩阵)至关重要。由于有限的基函数的自由度是分配给每个单元的边或点上的，所以一个单元中的物理场仅有此单元节点和边的物理场值与基函数决定，与其他单元的基函数无关，一般称这种方法中的基函数为局部基函数(Local basis function)。所以构建有限元方程中的矩阵时，一般是根据计算矩阵元素的公式先计算每个单元的单元矩阵(Element matrix)，最后再将各个单元的矩阵组装起来，得到最终需要的全局矩阵(Global matrix)。
S3，通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；
由于有限元法使用的基函数为局部基函数(Local basis function)，此基函数只在包含定义其自由度的点或边的单元中场值不为零，所以不相邻单元中的基函数一般是正交的，根据热学时域有限元的基本方程，利用这种有限元基函数得到的矩阵一般是稀疏矩阵。在Matlab中实现有限元矩阵构造时，可以采用Matlab中的稀疏矩阵(Sparse matrix)的数据格式即相关的数学方法，这样可以减少程序的内存占用并加快程序运行速度。
S4，通过热学时域有限元的基本方程对每一个时间步进行迭代运算；
在进行矩阵线性运算时，可以通过LU分解等矩阵分解方法来提高程序的运算性能(采用矩阵分解法求解大型线性方程组的效率远高于直接对大矩阵求逆)。在具体毕业设计中，采用LUPQ矩阵分解法。利用Matlab中的lu函数，将稀疏矩阵S分解为一个单位下三角矩阵L、一个上三角矩阵U、一个行置换矩阵P以及一个列置换矩阵Q，并满足：
PSQ＝LU
热学时域有限元的基本方程都可以转换线性方程组求解的一般形式：
Asx＝b
其中As为热学时域有限元的最终公式中的项，x为{T}t 项，b为项。
利用LUPQ矩阵分解方程：
PAsQ＝LU
可以得到用LUPQ表示的线性方程组：
P-1LUQ-1x＝b
从而得到x的解：
x＝PL-1U-1Qb
经过大量的实验，发现在利用Matlab进行线性方程组进行求解时，利用LUPQ方法的求解效率与不做矩阵分解的直接求解的效率相比有明显的优势，且在绝大多数情况下利用LUPQ求解方法进行线性问题求解的效率也比其他矩阵分解的效率高。
第二实施例；
本发明第二实施例是基于上述第一实施例形成，相同部分不再赘述，在上述第一实施例基础上增加下述步骤；
S5，通过验证工具对热学时域有限元的基本方程结果进行验证。
示例性的，验证和最终验证采用的参考软件为COMSOL Multiphysics。COMSOL Multiphysics最初也是于Matlab平台作为应用开发的，在工业界和学术界广泛用于对偏微分方程描述的物理问题通过某些数学方法进行数值求解。COMSOL Multiphysics支持使用的算法繁多，且有许多优化性算法以用于不同场景，但是这些算法本质上都属于有限元方法。作为一款跨平台的有限元分析，求解器和多物理场仿真软件。COMSOL Multiphysics允许通过传统的基于物理学的用户界面进行建模并构造偏微分方程的耦合系统，是一款使用自由度较高的软件。COMSOL为机械、流体、电气、射频、声学、电化学反应、等离子体和光学等不同领域的仿真研究工作提供IDE和统一的工作流程。
为验证热学时域有限元求解器的准确性，本发明构建了几个算例，并与COMSOL Multiphysics的结果进行比对。
第一个算例为铜块冷却：0.1m×0.1m×0.1m的铜块，其热学参数：导热系数为377W/(m·K)，密度为8900kg/m3，恒压热容为386J/(kg·K)，固体的边界条件为对流边界条件，对流流体的温度为300K，换热系数为500。
为观察其温度变化，在铜块中心设置一域点探针，获取改点的温度。在COMSOL Multiphysics和自行编写的FETD热学求解器中计算该程序，对比结果。其结果参考图2和图3所示。
根据计算结果可以看出，自行编写的FETD热学求解器得到的结果与COMSOL Multiphysics中的结果基本吻合，误差在可以接受的范围内。
第二个算例为搭载散热片的体热源，参考图4所示。有四个2cm×2cm×2cm的体热源，与9片厚度为0.5cm的薄散热片，通过厚度为1cm的底座连接在一起。整个物体处于对流边界条件中，由于此模型设计的相关参数较多，因此不再赘述其具体参数，将于附录中给出具体的参数。
在COMSOL Multiphysics和自行编写的FETD程序中对此算例进行仿真。对比仿真结果得出，COMSOL Multiphysics和自行编写的FETD程序的计算结果大体一致。通过两个算例的对比验证，可以得出本发明的仿真结构是可靠的。
第三实施例；
本发明提供一种用于执行第一实施例或第二实施例任意一项所述的时域热传导仿真方法中步骤的计算机可读存储介质。
除非另有定义，否则这里所使用的全部术语(包括技术术语和科学术语)都具有与本发明所属领域的普通技术人员通常理解的意思相同的意思。还将理解的是，除非这里明确定义，否则诸如在通用字典中定义的术语这类术语应当被解释为具有与它们在相关领域语境中的意思相一致的意思，而不以理想的或过于正式的含义加以解释。
本文中应采用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应 用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (9)

  1. 一种时域热传导仿真方法，其特征在于，包括以下步骤：
    S1，利用有限元网格剖分工具几何建模和网格剖分，读取几何信息；
    S2，设置固体传热学分析中各个区域材料相关参数；
    S3，构建热学时域有限元的基本方程计算出单元矩阵中每一个元素的值，热学采用基于四面体单元的线性点基标量基函数，获得每个单元的单元矩阵；
    S3，通过每个单元、节点、边的编号，将单元矩阵中的元素与全局矩阵的元素对应，通过遍历叠加利用单元矩阵，得到全局矩阵；
    S4，通过热学时域有限元的基本方程对每一个时间步进行迭代运算。
  2. 如权利要求1所述的时域热传导仿真方法，其特征在于，还包括：
    S5，通过验证工具对热学时域有限元的基本方程结果进行验证。
  3. 如权利要求1所述的时域热传导仿真方法，其特征在于：
    步骤S1中，通过有限元网格剖分工具进行几何建模和网格剖分，通过第三方仿真软件读取几何信息。
  4. 如权利要求1所述的时域热传导仿真方法，其特征在于：
    步骤S2中，各个区域材料相关参数包括：热导率，质量密度，比热容、边界条件以及观测数据所需的域点探针。
  5. 如权利要求1所述的时域热传导仿真方法，其特征在于：热学时域有限元的基本方程表达如下；
    T为温度，C为材料比热容，t为时间，Δt为仿真时间间隔，[CT]为阻尼矩阵，[KT]为刚度矩阵，{pT}为激励源，{T}为离散形式的温度场标量，{T}t为t时刻离散形式的温度场标量，{T}t-Δt为t-Δt时刻离散形式的温度场标量。
  6. 如权利要求1所述的时域热传导仿真方法，其特征在于：步骤S3中，每个单元的单元矩阵包括质量矩阵、阻尼矩阵和刚度矩阵。
  7. 如权利要求1所述的时域热传导仿真方法，其特征在于：步骤S4 中，进行迭代运算采用LUPQ矩阵分解法。
  8. 如权利要求1所述的时域热传导仿真方法，其特征在于：步骤S5中，通过构建算例采用第三方仿真软件作为验证工具。
  9. 一种用于执行权利要求1-8任意一项所述的时域热传导仿真方法中步骤的计算机可读存储介质。
PCT/CN2023/114983 2022-08-26 2023-08-25 时域热传导仿真方法及存储介质 WO2024041642A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202211032489.1 2022-08-26
CN202211032489.1A CN115374673A (zh) 2022-08-26 2022-08-26 时域热传导仿真方法及存储介质

Publications (1)

Publication Number Publication Date
WO2024041642A1 true WO2024041642A1 (zh) 2024-02-29

Family

ID=84068291

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2023/114983 WO2024041642A1 (zh) 2022-08-26 2023-08-25 时域热传导仿真方法及存储介质

Country Status (2)

Country Link
CN (1) CN115374673A (zh)
WO (1) WO2024041642A1 (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115374673A (zh) * 2022-08-26 2022-11-22 宁波德图科技有限公司 时域热传导仿真方法及存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090228246A1 (en) * 2008-03-05 2009-09-10 Livermore Software Technology Corporation Methods and systems of engineering analysis using a hybrid approach with FEM and adaptive SPH
CN102033985A (zh) * 2010-11-24 2011-04-27 南京理工大学 基于*-矩阵算法的高效时域电磁仿真方法
CN109190169A (zh) * 2018-08-02 2019-01-11 电子科技大学 一种三维时域电磁学杂交时域间断伽辽金数值方法
CN109492341A (zh) * 2018-12-25 2019-03-19 南京邮电大学 表面等离激元波导的光热效应仿真方法
CN114117864A (zh) * 2021-12-03 2022-03-01 厦门大学 自适应时间步长有限元法在电子器件热仿真中的应用方法
CN115374673A (zh) * 2022-08-26 2022-11-22 宁波德图科技有限公司 时域热传导仿真方法及存储介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090228246A1 (en) * 2008-03-05 2009-09-10 Livermore Software Technology Corporation Methods and systems of engineering analysis using a hybrid approach with FEM and adaptive SPH
CN102033985A (zh) * 2010-11-24 2011-04-27 南京理工大学 基于*-矩阵算法的高效时域电磁仿真方法
CN109190169A (zh) * 2018-08-02 2019-01-11 电子科技大学 一种三维时域电磁学杂交时域间断伽辽金数值方法
CN109492341A (zh) * 2018-12-25 2019-03-19 南京邮电大学 表面等离激元波导的光热效应仿真方法
CN114117864A (zh) * 2021-12-03 2022-03-01 厦门大学 自适应时间步长有限元法在电子器件热仿真中的应用方法
CN115374673A (zh) * 2022-08-26 2022-11-22 宁波德图科技有限公司 时域热传导仿真方法及存储介质

Also Published As

Publication number Publication date
CN115374673A (zh) 2022-11-22

Similar Documents

Publication Publication Date Title
Gao et al. Element differential method for solving general heat conduction problems
Sovinec et al. Nonlinear magnetohydrodynamics simulation using high-order finite elements
Shadid et al. Scalable implicit incompressible resistive MHD with stabilized FE and fully-coupled Newton–Krylov-AMG
Bochev et al. Finite element methods of least-squares type
Alauzet et al. Feature-based and goal-oriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace
Bazyar et al. Scaled boundary finite-element method for solving non-homogeneous anisotropic heat conduction problems
Antoniadis et al. Assessment of high-order finite volume methods on unstructured meshes for RANS solutions of aeronautical configurations
WO2024041642A1 (zh) 时域热传导仿真方法及存储介质
Diosady et al. Design of a variational multiscale method for high reynolds number compressible flows
Chamoin et al. Certified real‐time shape optimization using isogeometric analysis, PGD model reduction, and a posteriori error estimation
Fang et al. Isogeometric boundary element analysis for two-dimensional thermoelasticity with variable temperature
Giacomini et al. HDGlab: An open-source implementation of the hybridisable discontinuous Galerkin method in MATLAB
Zhang DOF‐gathering stable generalized finite element methods for crack problems
Li et al. Proper orthogonal decomposition with SUPG-stabilized isogeometric analysis for reduced order modelling of unsteady convection-dominated convection-diffusion-reaction problems
Liang et al. A new alternating iteration strategy based on the proper orthogonal decomposition for solving large-scaled transient nonlinear heat conduction problems
Veilleux et al. A stable Spectral Difference approach for computations with triangular and hybrid grids up to the 6th order of accuracy
Camata et al. Reordering and incomplete preconditioning in serial and parallel adaptive mesh refinement and coarsening flow solutions
Williams et al. Nodal points and the nonlinear stability of high-order methods for unsteady flow problems on tetrahedral meshes
Zang et al. Isogeometric boundary element method for steady-state heat transfer with concentrated/surface heat sources
Zang et al. Isogeometric boundary element for analyzing steady-state heat conduction problems under spatially varying conductivity and internal heat source
Zhan et al. Three-dimensional high-order finite-volume method based on compact WENO reconstruction with hybrid unstructured grids
Chen et al. Shape optimization of fluid cooling channel based on Darcy reduced-order isogeometric analysis
Ramšak et al. 3D multidomain BEM for solving the Laplace equation
Li et al. A new approach to solve multi-medium nonlinear transient heat conduction problems using interface integration BEM
Guventurk et al. An arbitrary Lagrangian‐Eulerian framework with exact mass conservation for the numerical simulation of 2D rising bubble problem

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23856727

Country of ref document: EP

Kind code of ref document: A1