CN118410663A - 一种高效求解非线性热边界问题的改良弧长法 - Google Patents
一种高效求解非线性热边界问题的改良弧长法 Download PDFInfo
- Publication number
- CN118410663A CN118410663A CN202410343630.2A CN202410343630A CN118410663A CN 118410663 A CN118410663 A CN 118410663A CN 202410343630 A CN202410343630 A CN 202410343630A CN 118410663 A CN118410663 A CN 118410663A
- Authority
- CN
- China
- Prior art keywords
- load
- increment
- vector
- follows
- unit
- Prior art date
- Legal status (The legal status 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 status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 99
- 238000004364 calculation method Methods 0.000 claims abstract description 90
- 238000012546 transfer Methods 0.000 claims abstract description 20
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 19
- 239000011159 matrix material Substances 0.000 claims description 99
- 238000012937 correction Methods 0.000 claims description 58
- 239000000463 material Substances 0.000 claims description 44
- 230000010354 integration Effects 0.000 claims description 15
- 238000002076 thermal analysis method Methods 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 9
- 238000004458 analytical method Methods 0.000 claims description 7
- 239000006185 dispersion Substances 0.000 claims description 7
- 230000004907 flux Effects 0.000 claims description 6
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000009191 jumping Effects 0.000 claims description 3
- 238000012805 post-processing Methods 0.000 claims description 3
- 238000003672 processing method Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 238000011217 control strategy Methods 0.000 abstract description 4
- 238000004088 simulation Methods 0.000 abstract description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000013461 design Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005272 metallurgy Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种高效求解非线性热边界问题的改良弧长法,属于高性能计算仿真技术领域。所述方法首先需要建立待分析的几何模型,其次对待分析的几何模型进行任意单元类型的网格划分,然后设置待分析模型的边界条件及荷载情况,最后将上述条件输入至算法求解器中。在求解中,通过判定相邻增量步预测阶段所求得的切线增量的二范数大小。依据判定结果选择算法预测阶段弧长计算方法。从而实现传热非线性问题可靠求解的前提下,对算法求解过程进行加速。本发明将传热学理论、有限元方法、非线性算法控制策略结合起来,在现有弧长法的理论基础上进行改良,以计算机卓越的数值计算能力为载体来求解热学非线性问题,不仅提升了计算结果的数值精度,而且提升了热学非线性问题求解速度。
Description
技术领域
本发明属于高性能计算仿真技术领域,具体是一种高效求解非线性热边界问题的改良弧长法。
背景技术
热传递无处不在。传热不仅是常见的自然现象,而且广泛存在于工程技术领域。在能源动力、化工制药、材料冶金、机械制造、部门中均存在热传递影响突出的问题
目前,虽然传热学的理论已逐渐完善,但是由于传热问题处理的问题多为非线性的实际问题,一般通过将非线性问题细化为多个荷载增量步,通过将传统的以直代曲的求解思想应用到每个单个荷载增量步中,以多次的线性求解的迭代方式来对非线性问题进行逼近。而随着电子计算机的软硬件技术的飞速发展以及数值仿真理论和算法的不断优化。在非线性热分析方面求解的数值精度、算法收敛的可靠性、求解速度的要求不断提高,以牛顿迭代法为原型的求解算法不断被改良升级,出现了一批适用性很强的非线性热问题求解的算法。并且通过改良荷载以及温度的控制策略、细化荷载增量步、在求解过程中调节求解速度的方式,在保证计算结果数值精度的前提下,在算法层面进行了速度的提升。
因此,本发明建立一种高效求解非线性热边界问题的改良弧长法,利用发展成熟的有限元数值计算理论、传热学理论、非线性求解策略,着力于求解热学非线性问题,在准确求解非线性热问题同时,提升求解速度,不仅可以应用于航空航天等高端装备大区域传热分析,还可以在为芯片设计、微纳制造等高精端领域的非线性问题提供更加精确的求解方法,能够极好的促进材料和结构设计的热学以及热力学性能优化。
发明内容
本发明公开了一种高效求解非线性热边界问题的改良弧长法,该方法将传热学理论、有限元方法、非线性算法控制策略结合起来,在现有弧长法的理论基础上进行改良,以计算机卓越的数值计算能力为载体来求解热学非线性问题,不仅提升了计算结果的数值精度,而且提升了热学非线性问题求解速度。
本发明解决上述技术问题所提供的技术方案是:一种高效求解非线性热边界问题的改良弧长法,包括以下步骤:
S1、建立待定热分析对象的几何模型或者几何装配体模型;
S2、对待分析几何模型进行空间离散化,即网格划分,生成计算所需的网格数据;
S3、对划分的网格模型进行边界条件、载荷、材料导热系数等参数设定;
S4、设定改良弧长法首个增量步初始荷载因子数值的大小λinitial、迭代容差值大小t、荷载因子符号S、总荷载因子数值大小λ、Hammer积分点个数N;
S5、使用改良弧长法针对既定荷载边界条件下的问题求解,需计算当前参考荷载热流率向量{Qref}。进入荷载增量步预测阶段求解。计算单元的切线传导矩阵、当前构型下各个单元的热导率、切线预测增量求解。判断相邻荷载增量步切线增量的二范数大小,依据判定结果确定加速方法;
S6、进入校正求解阶段,首先更新当前外荷载向量{Qe}、内部热流率向量{Qa}、残差向量其次判断当前荷载增量步是否达到收敛标准。最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小;
S7、完成热分析求解,收集最终温度场的温度向量结果;
S8、计算单元传导矩阵、并计算温度梯度、热流,实现结果后处理;
S9、计算节点的平均热流完成既定荷载边界条件下热分析非线性问题求解;记录求解器总耗时。
进一步的技术方案是,所述步骤S1具体实现方法为:
建立宏观尺度任意规模的几何模型或者几何装配体模型,建立一个连续的求解域;
进一步的技术方案是,所述步骤S2具体实现方法为:
S21、离散参数的设定;
结合大规模的几何连续模型的几何特征,设定几何模型的全局网格种子、部分特征几何边的局部种子,在部分几何特征较复杂的特征边处设置局部细化参数;
S22、对传热现象问题进行数学建模,会产生偏微分方程(PDE),数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解,其离散方程
K(T)=Q,Q∈Rn
式中:n是离散化系统的总自由度,K是在当前边界荷载问题下系统的线性方程组系数矩阵,Q为荷载向量,可使用线性代数方法求解未知量u;
S23、建立离散化几何模型;
设定网格划分参数,使用网格划分工具开展几何模型网格划分,可实现不同形状的大规模几何模型的空间离散,将模型离散为以四面体为单元的空间连续离散体,同时将模型离散后的数据信息写出至本地文件;
进一步的技术方案是,所述步骤S3具体实现方法为:
S31、设置下列材料参数,材料导热系数K、材料导热系数与温度相关的方程K(T)、环境温度Ta等,并将求解所需材料参数输出到本地文件;
S32、对离散后的模型施加边界条件及载荷,指定边界约束的几何位置、载荷大小以及载荷施加的几何位置。并将受约束载荷的几何位置集合信息输出到本地文件;
进一步的技术方案是,所述步骤S4具体实现方法为:
S41、改良弧长法求解时需设置初始荷载因子值的大小,荷载因子初始值会影响求解收敛性。迭代容差值是校正阶段求解时迭代步收敛标准。荷载因子符号S默认为正。荷载因子最大值λmax为1,即满载情况一参考荷载向量{Qref}等于{Qe}外部热流率向量。预设初始荷载因子值为0.1以及迭代容差值为10-3。按照四面体数值积分精度要求,预设Hammer积分的积分点个数N为4个。
进一步的技术方案是,所述步骤S5具体实现方法为:
S51、进入改良弧长法求解器之前,需要读取模型信息,其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储。首先根据既定载荷以及边界情况,计算当前构型下的参考荷载热流率向量{Qref}。其次从本地文件中读取既定的初始荷载因子参数λinitial、迭代容差t、荷载因子最大值λmax、数值积分点个数N。
具体计算步骤如下:
首先在第i个荷载增量步下第j个迭代步,读取既定荷载以及边界情况,将受约束自由度的边界情况及载荷大小进行存储;
其次依据有限元理论对荷载数值大小在所约束区域上进行计算。
式中:为第i个荷载增量步下第j个迭代步的整体传导矩阵;Ω为单元域;其中h为对流换热系数;N为单元形函数;Γ为模型离散后荷载加载区域。G为体热流率密度大小,q为面热流密度大小,Ta为环境温度。
S52、进入改良弧长法求解器荷载增量步主循环,首先遍历模型单元,对离散后模型每个单元进行节点温度的提取,依据材料热导率与温度关系式K(T)计算当前单元高斯点上导热系数K。其次采用数值积分的方法对第i个荷载增量步的第j个迭代步下的单元切线传导矩阵进行求解。
具体计算步骤如下:
首先计算当前单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s,t)下四面体单元的形函数N;
其次形函数在母单元局部坐标系下对各个坐标分量求偏导,获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量
然后将局部坐标系(r,s,t)下的形函数与整体坐标系(x,y,z)下的形函数通过形函数偏导向量通过Jacobian[J]建立联系,其计算公式为:
Ni~l为单元形函数,下标i,j,k,l是四面体单元中节点编号;xi~l为单元节点在整体坐标下x轴的坐标分量;
最后获得梯度矩阵[B]为:
各项同性材料导热矩阵[D]为:
式中:Kx、Ky、Kz是材料的固有属性,数值大小由材料的热导率指定。
单元有限元线性热传导矩阵的计算公式为:
式中:为第i个荷载增量步的第j个迭代步下的单元导热矩阵,Ω为单元域。
以单个四面体单元为例有限元方程具体形式如下:
针对非线性问题时,第i个荷载增量步下第j个迭代步下的单元切线传导矩阵计算具体步骤如下:
当模型单元类型为一阶四面体单元,单元高斯点温度计算公式为:
Tgauss=NiTi+NjTj+NkTk+NlTl
依据材料参数设置步骤中的热导率与温度的关系式K(T),将高斯点的温度代入到K(T)中可得当前构型下单元导热系数。
各项同性材料导热矩阵线性项[D]为:
式中:K(T)x、K(T)y、K(T)z是根据当前材料所处单元的高斯积分点的温度值通过热导率与温度的关系式K(T)计算所得。
各项同性材料导热矩阵非线性项[Dn]为:
梯度矩阵[B]为;
至此,非线性求解时在第i个荷载增量步下第j个迭代步的单元切线传导矩阵线性项及非线性项计算方法如下:
线性项的计算公式为[Kl]e=∫Ω[B]T[D][B]dΩ,具体形式如下:
非线性项的计算公式为[Knl]e=∫Ω[B]T[Dn][B][T][N]dΩ,具体形式如下:
因此单元的切线传导矩阵为并对各个矩阵进行存储;
S53、将上述整体坐标系下的单元切线传导矩阵,保存其非0元素的位置索引及其对应位置矩阵元素的值,使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储;完成求解域在第i个荷载增量步下第j个迭代步的整体切线传导矩阵的组装。
S54、通过置大数法对第i个荷载增量步下第j个迭代步的整体切线传导矩阵外部荷载热流率向量{Qe}进行处理。依据边界条件,在受约束的点的自由度上放置一个大数值如108甚至更大,其意义为让受约束的自由度在求解时等于设定的边界值dofvalue。dofvalue从模型离散后的文本文件将受约束的自由度和边界数值读取可得。
置大数法具体形式如下:
S55、配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小。将第i个荷载增量步下第j个迭代步处理过后的整体切线传导矩阵外部热流率向量{Qe}传入,调用语句进行线性方程组的求解,求得预测阶段切线增量方程形式为:
S56、判断当前荷载增量步是否为首个荷载增量步,若当前增量步为首个荷载增量步时,则不需要计算弧长大小。若当前增量步不为首个荷载增量步时,需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量的二范数的大小,依据判定结果对算法进行加速。其中上标i为当前荷载增量步,下标1代表为预测阶段迭代步。
弧长计算公式为:ΔSi={ΔTi-1}T{ΔTi-1}式中上标i代表当前荷载增量步;
对于环境温度Ta非负的热传导分析时,当该增量步预测阶段切线增量二范数满足上述要求时,此种情况可将弧长计算公式进行如下修改:
(n为系统自由度)
弧长ΔS具体的计算表达式如下:
(n为系统自由度)
式中:{ΔT}T上标T为转置含义,ΔT为荷载增量步温度增量;上标i代表荷载增量步步数,下标1代表预测阶段迭代步;
S57、若当前荷载增量步为首个荷载增量步,荷载因子大小为初始荷载因子大小λinitial,记录荷载增量步预测阶段切线增量向量。若当前荷载增量步不为首个荷载增量步荷载因子具体计算步骤如下:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号。
S58、预测阶段荷载因子计算完成,保存当前荷载增量步预测阶段切线增量向量,更新温度向量,以及荷载因子λ大小。结束荷载增量步预测阶段求解。
计算公式为:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号
进一步的技术方案是,所述步骤S6具体实现方法为:
S61、首先依据更新的第i个荷载增量步下第j个迭代步的荷载因子更新外部荷载热流率向量以及总的荷载因子λ。具体计算步骤如下;
式中:i为荷载增量步;j为迭代步;为荷载因子增量步增量;为荷载因子迭代步增量
S62、其次遍历单元,计算第i个荷载增量步下第j个迭代步的单元的线性传导矩阵并组装成整体线性传导矩阵将整体线性传导矩阵乘以节点温度向量{T}得到单元内部热流率向量
具体计算步骤如下:
式中:整体线性传导矩阵的计算同S52所述,此处不再赘述。下标j-1是由于预测阶段并未进行计算。
S63、更新当前构型下的残差向量判断残差向量是否小于既定残差。此处的残差向量要依据边界条件对受约束的自由度位置进行处理,自由度受约束位置的内力为0;
第i个荷载增量步下第j个迭代步的残差向量计算步骤如下:
若则程序跳出校正阶段循环,若荷载因子λ小于λmax时,程序执行下一荷载增量步计算。反之,完成非线性热学问题的求解。若则进入迭代校正阶段。
S64、迭代校正阶段主要是为了消除系统内外热流率的不平衡。首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是温度迭代校正向量以及残差迭代校正向量求解,然后依据弧长约束方程更新荷载因子大小,最后更新总的位移向量以及总荷载因子大小。
其主要计算步骤如下:
算法中采用的约束方程为柱面约束方程,其方程形式为:
{ΔTi-1}T{ΔTi-1}=ΔS2
式中:i代表荷载增量步;ΔTi-1表示上一荷载增量步温度增量;
对于第i个荷载增量步下第j个迭代步有限元方程如下:
而对于第i个荷载增量步下第j个迭代步的外热流率载荷如下:
而在第i个荷载增量步下第j个迭代步的温度增量如下:
因此在第i个荷载增量步下第j个迭代步方程的最终形式如下:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
因此校正阶段求解公式为:
单元及整体切线传导矩阵的计算步骤同S52,此处不再赘述。
将外部热流率参考向量{Qref},更新后的残差向量依据边界条件进行处理,对自由度被约束的位置采用置大数法进行处理。处理方法同S54,此处不再赘述。
配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小。传入即可求得温度迭代校正向量同理可求得残差迭代校正向量
S65、迭代校正阶段荷载因子求解
在荷载增量步i下第j迭代步的温度增量可表示为如下形式:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
温度增量代入到柱面约束方程中可得以当前迭代步荷载因子为变量的一元二次方程。因此在校正阶段时关于荷载因子的一元二次方程为:
则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为:
式中:i代表荷载增量步步数;j代表迭代步步数;代表预测阶段切线增量;代表温度迭代校正向量;代表残差迭代校正向量;
S66、依据S56中的判定结果选择校正阶段计算方式,计算荷载因子更新温度向量以及总荷载因子λ大小。直至残差满足小于预设迭代容差停止迭代。
进一步的技术方案是,所述步骤S7具体实现方法为:
S71、当总荷载因子数值大小大于预设最大荷载因子大小时,跳出荷载增量步主循环,记录当前构型下总的温度向量。求解器求解完成并存储当前模型各个节点的温度向量。
进一步的技术方案是,所述步骤S8具体实现方法为:
S81、整个求解域的温度场计算完成后获得温度场向量,遍历单元计算求解温度梯度矩阵{g},计算公式如下,
式中:Ti、Tj、Tk、Tl为单元节点温度,中间矩阵为梯度矩阵[B]矩阵,Ni、Nj、Nk、Nl为单元形函数;
S82、计算单元热流梯度{q}e向量,整个求解域的温度梯度计算完成后,计算单元热流,热传导满足傅里叶定律,公式如下,
式中:qx,qy,qz为x,y,z方向的热通量;k是热导率,dT/dx,dT/dy,dT/dz是温度梯度;
进一步的技术方案是,所述步骤S9具体实现方法为:
S91、根据单元热通量{q}e计算节点平均热流,平均热流矢量计算公式为:
式中:为节点平均热流矢量;{q}e为包含该节点单元的单元热流大小,n为包含该节点的单元个数。
本发明的有益效果是:本发明基于LiToMesh完成几何模型的离散建立有限元模型,并以四面体单元类型对几何模型进行离散。然后将模型信息导入到自主研发的求解器中。通过判定相邻增量步预测阶段所求得的切线增量的二范数大小,依据判定结果对算法进行加速,从而实现传热非线性问题可靠求解的前提下,不使用并行技术对算法求解过程进行加速。不仅提升了计算结果的数值精度,而且提升了热学非线性问题求解速度。对求解速度有高要求的传热计算场景提供了解决方案。
附图说明
图1是本发明的一种高效求解非线性热边界问题的改良弧长法的流程图;
图2是本发明的几何计算模型示意图;
图3是离散单元类型的示意图;
图4是几何模型离散后的有限元网格图;
图5是非线性热流载荷与温度的关系图;
图6是求解所得的温度场结果;
图7是三维空间下各方向节点热流。
具体实施方式
下面将结合附图对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例是发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图1所示,本发明提供了一种高效求解非线性热边界问题的改良弧长法,包括以下步骤:
S1、建立待定热分析对象的几何模型或者几何装配体模型;
具体的实现步骤为:
建立宏观尺度任意规模的几何模型或者几何装配体模型,建立一个连续的求解域;几何模型示意图如图2所示;
S2、对待分析几何模型进行空间离散化,即网格划分,生成计算所需的网格数据;
具体的实现步骤为:
S21、结合大规模的几何连续模型的几何特征,设定几何模型的全局网格种子、部分特征几何边的局部种子,在部分几何特征较复杂的特征边处设置局部细化参数;
S22、对传热现象问题进行数学建模,会产生偏微分方程(PDE),数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解,其离散方程
K(T)=Q,Q∈Rn
式中:n是离散化系统的总自由度,K是在当前边界荷载问题下系统的线性方程组系数矩阵,Q为荷载向量,可使用线性代数方法求解未知量u;
S23、建立离散化几何模型;
设定网格划分参数,使用网格划分工具开展几何模型网格划分,可实现不同形状的大规模几何模型的空间离散,将模型离散为以四面体为单元的空间连续离散体,同时将模型离散后的数据信息写出至本地文件;单元模型如图3所示,离散后的模型如图4所示
S3、对划分的网格模型进行边界条件、载荷、材料导热系数等参数设定;
具体包括如下步骤:
S31、设置下列材料参数,材料导热系数K、材料导热系数与温度相关的方程K(T)、环境温度Ta等,并将求解所需材料参数输出到本地文件;
S32、对离散后的模型施加边界条件及载荷,指定边界约束的几何位置、载荷大小以及载荷施加的几何位置。并将受约束载荷的几何位置集合信息输出到本地文件;载荷大小与高斯点温度的非线性函数如图5所示;
S4、设定改良弧长法首个增量步初始荷载因子数值的大小λinitial、迭代容差值大小t、荷载因子符号S、总荷载因子数值大小λ、Hammer积分点个数N;
具体包括如下步骤:
S41、改良弧长法求解时需设置初始荷载因子值的大小,荷载因子初始值会影响求解收敛性。迭代容差值是校正阶段求解时迭代步收敛标准。荷载因子符号S默认为正。荷载因子最大值λmax为1,即满载情况一参考荷载向量{Qref}等于{Qe}外部热流率向量。预设初始荷载因子值为0.1以及迭代容差值为10-3。按照四面体数值积分精度要求,预设Hammer积分的积分点个数N为4个。
S5、使用改良弧长法针对既定荷载边界条件下的问题求解,需计算当前参考荷载热流率向量{Qref}。进入荷载增量步预测阶段求解。计算单元的切线传导矩阵、当前构型下各个单元的热导率、切线预测增量求解。判断相邻荷载增量步切线增量的二范数大小,依据判定结果确定加速方法;
具体包括如下步骤:
S51、进入改良弧长法求解器之前,需要读取模型信息,其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储。首先根据既定载荷以及边界情况,计算当前构型下的参考荷载热流率向量{Qref}。其次从本地文件中读取既定的初始荷载因子参数λinitial、迭代容差t、荷载因子最大值λmax、数值积分点个数N。
具体计算步骤如下:
首先在第i个荷载增量步下第j个迭代步,读取既定荷载以及边界情况,将受约束自由度的边界情况及载荷大小进行存储;
其次依据有限元理论对荷载数值大小在所约束区域上进行计算。
式中:为第i个荷载增量步下第j个迭代步的整体传导矩阵;Ω为单元域;其中h为对流换热系数;N为单元形函数;Γ为模型离散后荷载加载区域。G为体热流率密度大小,q为面热流密度大小,Ta为环境温度。
S52、进入改良弧长法求解器荷载增量步主循环,首先遍历模型单元,对离散后模型每个单元进行节点温度的提取,依据材料热导率与温度关系式K(T)计算当前单元高斯点上导热系数K。其次采用数值积分的方法对第i个荷载增量步的第j个迭代步下的单元切线传导矩阵进行求解。
具体计算步骤如下:
首先计算当前单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s,t)下四面体单元的形函数N;
其次形函数在母单元局部坐标系下对各个坐标分量求偏导,获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量
然后将局部坐标系(r,s,t)下的形函数与整体坐标系(x,y,z)下的形函数通过形函数偏导向量通过Jacobian[J]建立联系,其计算公式为:
Ni~l为单元形函数,下标i,j,k,l是四面体单元中节点编号;xi~l为单元节点在整体坐标下x轴的坐标分量;
最后获得梯度矩阵[B]为:
各项同性材料导热矩阵[D]为:
式中:Kx、Ky、Kz是材料的固有属性,数值大小由材料的热导率指定。
单元有限元线性热传导矩阵的计算公式为:
式中:为第i个荷载增量步的第j个迭代步下的单元导热矩阵,Ω为单元域。
以单个四面体单元为例有限元方程具体形式如下:
针对非线性问题时,第i个荷载增量步下第j个迭代步下的单元切线传导矩阵计算具体步骤如下:
当模型单元类型为一阶四面体单元,单元高斯点温度计算公式为:
Tgauss=NiTi+NjTj+NkTk+NlTl
依据材料参数设置步骤中的热导率与温度的关系式K(T),将高斯点的温度代入到K(T)中可得当前构型下单元导热系数。
各项同性材料导热矩阵线性项[D]为:
式中:K(T)x、K(T)y、K(T)z是根据当前材料所处单元的高斯积分点的温度值通过热导率与温度的关系式K(T)计算所得。
各项同性材料导热矩阵非线性项[Dn]为:
梯度矩阵[B]为;
至此,非线性求解时在第i个荷载增量步下第j个迭代步的单元切线传导矩阵线性项及非线性项计算方法如下:
线性项的计算公式为[Kl]e=∫Ω[B]T[D][B]dΩ,具体形式如下:
非线性项的计算公式为[Knl]e=∫Ω[B]T[Dn][B][T][N]dΩ,具体形式如下:
因此单元的切线传导矩阵为并对各个矩阵进行存储;
S53、将上述整体坐标系下的单元切线传导矩阵,保存其非0元素的位置索引及其对应位置矩阵元素的值,使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储;完成求解域在第i个荷载增量步下第j个迭代步的整体切线传导矩阵的组装。
S54、通过置大数法对第i个荷载增量步下第j个迭代步的整体切线传导矩阵外部荷载热流率向量{Qe}进行处理。依据边界条件,在受约束的点的自由度上放置一个大数值如108甚至更大,其意义为让受约束的自由度在求解时等于设定的边界值dofvalue。dofvalue从模型离散后的文本文件将受约束的自由度和边界数值读取可得。
置大数法具体形式如下:
S55、配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小。将第i个荷载增量步下第j个迭代步处理过后的整体切线传导矩阵外部热流率向量{Qe}传入,调用语句进行线性方程组的求解,求得预测阶段切线增量方程形式为:
S56、判断当前荷载增量步是否为首个荷载增量步,若当前增量步为首个荷载增量步时,则不需要计算弧长大小。若当前增量步不为首个荷载增量步时,需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量的二范数的大小,依据判定结果对算法进行加速。其中上标i为当前荷载增量步,下标1代表为预测阶段迭代步。
弧长计算公式为:ΔSi={ΔTi-1}T{ΔTi-1}式中上标i代表当前荷载增量步;
对于环境温度Ta非负的热传导分析时,当该增量步预测阶段切线增量二范数满足上述要求时,此种情况可将弧长计算公式进行如下放大:
(n为系统自由度)
弧长ΔS具体的计算表达式如下:
(n为系统自由度)
式中:{ΔT}T上标T为转置含义,ΔT为荷载增量步温度增量;上标i代表荷载增量步步数,下标1代表预测阶段迭代步;
S57、若当前荷载增量步为首个荷载增量步,荷载因子大小为初始荷载因子大小λinitial,记录荷载增量步预测阶段切线增量向量。若当前荷载增量步不为首个荷载增量步荷载因子具体计算步骤如下:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号。
S58、预测阶段荷载因子计算完成,保存当前荷载增量步预测阶段切线增量向量,更新温度向量,以及荷载因子λ大小。结束荷载增量步预测阶段求解。
计算公式为:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号
S6、进入校正求解阶段,首先更新当前外荷载向量{Qe}、内部热流率向量{Qa}、残差向量其次判断当前荷载增量步是否达到收敛标准。最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小;
具体包括如下步骤:
S61、首先依据更新的第i个荷载增量步下第j个迭代步的荷载因子更新外部荷载热流率向量以及总的荷载因子λ。具体计算步骤如下;
式中:i为荷载增量步;j为迭代步;为荷载因子增量步增量;为荷载因子迭代步增量
S62、其次遍历单元,计算第i个荷载增量步下第j个迭代步的单元的线性传导矩阵并组装成整体线性传导矩阵将整体线性传导矩阵乘以节点温度向量{T}得到单元内部热流率向量
具体计算步骤如下:
式中:整体线性传导矩阵的计算同S52所述,此处不再赘述。下标j-1是由于预测阶段并未进行计算。
S63、更新当前构型下的残差向量判断残差向量是否小于既定残差。此处的残差向量要依据边界条件对受约束的自由度位置进行处理,自由度受约束位置的内力为0;
第i个荷载增量步下第j个迭代步的残差向量计算步骤如下:
若则程序跳出校正阶段循环,若荷载因子λ小于λmax时,程序执行下一荷载增量步计算。反之,完成非线性热学问题的求解。若则进入迭代校正阶段。
S64、迭代校正阶段主要是为了消除系统内外热流率的不平衡。首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是温度迭代校正向量以及残差迭代校正向量求解,然后依据弧长约束方程更新荷载因子大小,最后更新总的位移向量以及总荷载因子大小。
其主要计算步骤如下:
算法中采用的约束方程为柱面约束方程,其方程形式为:
{ΔTi-1}T{ΔTi-1}=ΔS2
式中:i代表荷载增量步;ΔTi-1表示上一荷载增量步温度增量;
对于第i个荷载增量步下第j个迭代步有限元方程如下:
而对于第i个荷载增量步下第j个迭代步的外热流率载荷如下:
而在第i个荷载增量步下第j个迭代步的温度增量如下:
因此在第i个荷载增量步下第j个迭代步方程的最终形式如下:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
因此校正阶段求解公式为:
单元及整体切线传导矩阵的计算步骤同S52,此处不再赘述。
将外部热流率参考向量{Qref},更新后的残差向量依据边界条件进行处理,对自由度被约束的位置采用置大数法进行处理。处理方法同S54,此处不再赘述。
配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小。传入即可求得温度迭代校正向量同理可求得残差迭代校正向量
S65、迭代校正阶段荷载因子求解
在荷载增量步i下第j迭代步的温度增量可表示为如下形式:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
温度增量代入到柱面约束方程中可得以当前迭代步荷载因子为变量的一元二次方程。因此在校正阶段时关于荷载因子的一元二次方程为:
则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为:
式中:i代表荷载增量步步数;j代表迭代步步数;代表预测阶段切线增量;代表温度迭代校正向量;代表残差迭代校正向量;
S66、依据S56中的判定结果选择校正阶段计算方式,计算荷载因子更新温度向量以及总荷载因子λ大小。直至残差满足小于预设迭代容差停止迭代。
S7、完成热分析求解,收集最终温度场的温度向量结果;
具体包括如下步骤:
S71、当总荷载因子数值大小大于预设最大荷载因子大小时,跳出荷载增量步主循环,记录当前构型下总的温度向量。求解器求解完成并存储当前模型各个节点的温度向量。温度场结果如图6所示;
S8、计算单元传导矩阵、并计算温度梯度、热流,实现结果后处理;
具体包括如下步骤:
S81、整个求解域的温度场计算完成后获得温度场向量,遍历单元计算求解温度梯度矩阵{g},计算公式如下,
式中:Ti、Tj、Tk、Tl为单元节点温度,中间矩阵为梯度矩阵[B]矩阵,Ni、Nj、Nk、Nl为单元形函数;
S82、计算单元热流梯度{q}e向量,整个求解域的温度梯度计算完成后,计算单元热流,热传导满足傅里叶定律,公式如下,
式中:qx,qy,qz为x,y,z方向的热通量;k是热导率,dT/dx,dT/dy,dT/dz是温度梯度;
S9、计算节点的平均热流完成既定荷载边界条件下热分析非线性问题求解;记录求解器总耗时。
具体包括如下步骤:
S91、根据单元热通量{q}e计算节点平均热流,平均热流矢量计算公式为:
式中:为节点平均热流矢量;{q}e为包含该节点单元的单元热流大小,n为包含该节点的单元个数。节点热流计算结果如图7所示;各方法求解耗时统计如表1所示。
表1改良弧长法求解非线性传热问题耗时表
本发明公开了一种高效求解非线性热边界问题的改良弧长法,该方法将传热学理论、有限元方法、非线性算法控制策略结合起来,在现有弧长法的理论基础上进行改良,以计算机卓越的数值计算能力为载体来求解热学非线性问题,不仅提升了计算结果的数值精度,而且提升了热学非线性问题求解速度。
以上所述,并非对本发明作任何形式上的限制,虽然本发明已通过上述实施例揭示,然而并非用以限定本发明,任何熟悉本专业的技术人员,在不脱离本发明技术方案范围内,可利用上述揭示的技术内容做出变动或修饰为等同变化的等效实施例。但凡是未脱离本发明技术方案的内容,依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰,均仍属于本发明技术方案的范围内。
Claims (10)
1.一种高效求解非线性热边界问题的改良弧长法,其特征在于,包括以下步骤:
S1、建立待定热分析对象的几何模型或者几何装配体模型;
S2、对待分析几何模型进行空间离散化,即网格划分,生成计算所需的网格数据;
S3、对划分的网格模型进行边界条件、载荷、材料导热系数等参数设定;
S4、设定改良弧长法首个增量步初始荷载因子数值的大小λinitial、迭代容差值大小t、荷载因子符号S、总荷载因子数值大小λ、Hammer积分点个数N;
S5、使用改良弧长法针对既定荷载边界条件下的问题求解,需计算当前参考荷载热流率向量{Qref};进入荷载增量步预测阶段求解;计算单元的切线传导矩阵、当前构型下各个单元的热导率、切线预测增量求解;判断相邻荷载增量步切线增量的二范数大小,依据判定结果确定加速方法;
S6、进入校正求解阶段,首先更新当前外荷载向量{Qe}、内部热流率向量{Qa}、残差向量其次判断当前荷载增量步是否达到收敛标准;然后依据判断结果进入迭代校正求解流程;最终迭代阶段完成后判断总荷载因子大小是否超过预设最大荷载因子数值大小;
S7、完成热分析求解,收集最终温度场的温度向量结果;
S8、计算单元传导矩阵、并计算温度梯度、热流,实现结果后处理;
S9、计算节点的平均热流完成既定荷载边界条件下热分析非线性问题求解;记录求解器总耗时。
2.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S1具体实现方法为:建立宏观尺度任意规模的几何模型或者几何装配体模型,建立一个连续的求解域。
3.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S2具体实现方法为:
S21、离散参数的设定;
结合大规模的几何连续模型的几何特征,设定几何模型的全局网格种子、部分特征几何边的局部种子,在部分几何特征较复杂的特征边处设置局部细化参数;
S22、对传热现象问题进行数学建模,会产生偏微分方程(PDE),数值分析即有限元方法计算使得偏微分方程能够被近似求解原来模型域的数值解,其离散方程
K(T)=Q,Q∈Rn
式中:n是离散化系统的总自由度,K是在当前边界荷载问题下系统的线性方程组系数矩阵,Q为荷载向量,可使用线性代数方法求解未知量T;
S23、建立离散化几何模型;
设定网格划分参数,使用网格划分工具开展几何模型网格划分,可实现不同形状的大规模几何模型的空间离散,将模型离散为以四面体为单元的空间连续离散体,同时将模型离散后的数据信息写出至本地文件。
4.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S3具体实现方法为:
S31、设置下列材料参数,材料导热系数K、材料导热系数与温度相关的方程K(T)、环境温度Ta等,并将求解所需材料参数输出到本地文件;
S32、对离散后的模型施加边界条件及载荷,指定边界约束的几何位置、载荷大小以及载荷施加的几何位置;并将受约束载荷的几何位置集合信息输出到本地文件。
5.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S4具体实现方法为:
S41、改良弧长法求解时需设置初始荷载因子值的大小,荷载因子初始值会影响求解收敛性;迭代容差值是校正阶段求解时迭代步收敛标准;荷载因子符号S默认为正;荷载因子最大值λmax为1,即满载情况—参考荷载向量{Qref}等于{Qe}外部热流率向量;预设初始荷载因子值为0.1以及迭代容差值为10-3;按照四面体数值积分精度要求,预设Hammer积分的积分点个数N为4个。
6.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S5具体实现方法为:
S51、进入改良弧长法求解器之前,需要读取模型信息,其中包括离散后模型的单元节点信息、材料信息、边界及初始荷载因子并存储;首先根据既定载荷以及边界情况,计算当前构型下的参考荷载热流率向量{Qref};其次从本地文件中读取既定的初始荷载因子参数λinitial、迭代容差t、荷载因子最大值λmax、数值积分点个数N;
具体计算步骤如下:
首先在第i个荷载增量步下第j个迭代步,读取既定荷载以及边界情况,将受约束自由度的边界情况及载荷大小进行存储;
其次依据有限元理论对荷载数值大小在所约束区域上进行计算;
式中:为第i个荷载增量步下第j个迭代步的整体传导矩阵;Ω为单元域;其中h为对流换热系数;N为单元形函数;Γ为模型离散后荷载加载区域;G为体热流率密度大小,q为面热流密度大小,Ta为环境温度;
S52、进入改良弧长法求解器荷载增量步主循环,首先遍历模型单元,对离散后模型每个单元进行节点温度的提取,依据材料热导率与温度关系式K(T)计算当前单元高斯点上导热系数K;其次采用数值积分的方法对第i个荷载增量步的第j个迭代步下的单元切线传导矩阵进行求解;
具体计算步骤如下:
首先计算当前单元在第i个荷载增量步下第j迭代步的局部坐标系(r,s,t)下四面体单元的形函数N;
其次形函数在母单元局部坐标系下对各个坐标分量求偏导,获得母单元在局部坐标系下形函数对各个坐标轴的偏导向量
然后将局部坐标系(r,s,t)下的形函数与整体坐标系(x,y,z)下的形函数通过形函数偏导向量通过Jacobian[J]建立联系,其计算公式为:
Ni~l为单元形函数,下标i,j,k,l是四面体单元中节点编号;xi~l为单元节点在整体坐标下x轴的坐标分量;
最后获得梯度矩阵[B]为:
各项同性材料导热矩阵[D]为:
式中:Kx、Ky、Kz是材料的固有属性,数值大小由材料的热导率指定;
单元有限元线性热传导矩阵的计算公式为:
式中:为第i个荷载增量步的第j个迭代步下的单元导热矩阵,Ω为单元域;
以单个四面体单元为例有限元方程具体形式如下:
针对非线性问题时,第i个荷载增量步下第j个迭代步下的单元切线传导矩阵计算具体步骤如下:
当模型单元类型为一阶四面体单元,单元高斯点温度计算公式为:
Tgauss=NiTi+NjTj+NkTk+NlTl
依据材料参数设置步骤中的热导率与温度的关系式K(T),将高斯点的温度代入到K(T)中可得当前构型下单元导热系数;
各项同性材料导热矩阵线性项[D]为:
式中:K(T)x、K(T)y、K(T)z是根据当前材料所处单元的高斯积分点的温度值通过热导率与温度的关系式K(T)计算所得;
各项同性材料导热矩阵非线性项[Dn]为:
梯度矩阵[B]为;
至此,非线性求解时在第i个荷载增量步下第j个迭代步的单元切线传导矩阵线性项及非线性项计算方法如下:
线性项的计算公式为[Kl]e=∫Ω[B]T[D][B]dΩ,具体形式如下:
非线性项的计算公式为[Knl]e=∫Ω[B]T[Dn][B][T][N]dΩ,具体形式如下:
因此单元的切线传导矩阵为并对各个矩阵进行存储;
S53、将上述整体坐标系下的单元切线传导矩阵,保存其非0元素的位置索引及其对应位置矩阵元素的值,使用稀疏矩阵CSR方式中的COO存储方式进行稀疏矩阵的存储;完成求解域在第i个荷载增量步下第j个迭代步的整体切线传导矩阵的组装;
S54、通过置大数法对第i个荷载增量步下第j个迭代步的整体切线传导矩阵外部荷载热流率向量{Qe}进行处理;依据边界条件,在受约束的点的自由度上放置一个大数值如108甚至更大,其意义为让受约束的自由度在求解时等于设定的边界值dofvalue;dofvalue从模型离散后的文本文件将受约束的自由度和边界数值读取可得;
置大数法具体形式如下:
S55、配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小;将第i个荷载增量步下第j个迭代步处理过后的整体切线传导矩阵外部热流率向量{Qe}传入,调用语句进行线性方程组的求解,求得预测阶段切线增量方程形式为:
S56、判断当前荷载增量步是否为首个荷载增量步,若当前增量步为首个荷载增量步时,则不需要计算弧长大小;若当前增量步不为首个荷载增量步时,需判断上一个荷载增量步预测阶段的切线增量与当前荷载增量步的切线增量的二范数的大小,依据判定结果对算法进行加速;其中上标i为当前荷载增量步,下标1代表为预测阶段迭代步;
弧长计算公式为:ΔSi={ΔTi-1}T{ΔTi-1}式中上标i代表当前荷载增量步;
对于环境温度Ta非负的热传导分析时;当该增量步预测阶段切线增量二范数满足上述要求时,此种情况可将弧长计算公式进行如下放大:(n为系统自由度)
弧长ΔS具体的计算表达式如下:
(n为系统自由度)
式中:{ΔT}T上标T为转置含义,ΔT为荷载增量步温度增量;上标i代表荷载增量步步数,下标1代表预测阶段迭代步;
S57、若当前荷载增量步为首个荷载增量步,荷载因子大小为初始荷载因子大小λinitial,记录荷载增量步预测阶段切线增量向量;若当前荷载增量步不为首个荷载增量步荷载因子具体计算步骤如下:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号;
S58、预测阶段荷载因子计算完成,保存当前荷载增量步预测阶段切线增量向量,更新温度向量,以及荷载因子λ大小;结束荷载增量步预测阶段求解;
计算公式为:
式中:为预测阶段所求得的切线增量;其中i为荷载增量步;T为转置符号。
7.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S6具体实现方法为:
S61、首先依据更新的第i个荷载增量步下第j个迭代步的荷载因子更新外部荷载热流率向量以及总的荷载因子λ;具体计算步骤如下;
式中:i为荷载增量步;j为迭代步;为荷载因子增量步增量;为荷载因子迭代步增量
S62、其次遍历单元,计算第i个荷载增量步下第j个迭代步的单元的线性传导矩阵并组装成整体线性传导矩阵将整体线性传导矩阵乘以节点温度向量{T}得到单元内部热流率向量
具体计算步骤如下:
式中:整体线性传导矩阵的计算同S52所述,此处不再赘述;下标j-1是由于预测阶段并未进行计算;
S63、更新当前构型下的残差向量判断残差向量是否小于既定残差;此处的残差向量要依据边界条件对受约束的自由度位置进行处理,自由度受约束位置的内力为0;
第i个荷载增量步下第j个迭代步的残差向量计算步骤如下:
若则程序跳出校正阶段循环,若荷载因子λ小于λmax时,程序执行下一荷载增量步计算;反之,完成非线性热学问题的求解;若则进入迭代校正阶段;
S64、迭代校正阶段主要是为了消除系统内外热流率的不平衡;首先是计算计算第i个荷载增量步下第j个迭代步的整体切线传导矩阵其次是温度迭代校正向量以及残差迭代校正向量求解,然后依据弧长约束方程更新荷载因子大小,最后更新总的位移向量以及总荷载因子大小;
其主要计算步骤如下:
算法中采用的约束方程为柱面约束方程,其方程形式为:
式中:i代表荷载增量步;ΔTi-1表示上一荷载增量步温度增量;
对于第i个荷载增量步下第j个迭代步有限元方程如下:
而对于第i个荷载增量步下第j个迭代步的外热流率载荷如下:
而在第i个荷载增量步下第j个迭代步的温度增量如下:
因此在第i个荷载增量步下第j个迭代步方程的最终形式如下:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
因此校正阶段求解公式为:
单元及整体切线传导矩阵的计算步骤同S52,此处不再赘述;
将外部热流率参考向量{Qref},更新后的残差向量依据边界条件进行处理,对自由度被约束的位置采用置大数法进行处理;处理方法同S54,此处不再赘述;
配置第三方库Eigen,设定线性方程组求解器为共轭梯度求解器<BICGSTAB>,初始化共轭梯度求解器,设置求解器数值精度大小;传入即可求得温度迭代校正向量同理可求得残差迭代校正向量
S65、迭代校正阶段荷载因子求解
在荷载增量步i下第j迭代步的温度增量可表示为如下形式:
式中:i代表荷载增量步步数;j代表迭代步步数;代表温度迭代校正向量;代表残差迭代校正向量;
温度增量代入到柱面约束方程中可得以当前迭代步荷载因子为变量的一元二次方程;因此在校正阶段时关于荷载因子的一元二次方程为:
则第i个荷载增量步下第j个迭代步校正阶段的荷载因子计算公式为:
式中:i代表荷载增量步步数;j代表迭代步步数;代表预测阶段切线增量;代表温度迭代校正向量;代表残差迭代校正向量;
S66、依据S56中的判定结果选择校正阶段计算方式,计算荷载因子更新温度向量以及总荷载因子λ大小;直至残差满足小于预设迭代容差停止迭代。
8.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S7具体实现方法为:
S71、当总荷载因子数值大小大于预设最大荷载因子大小时,跳出荷载增量步主循环,记录当前构型下总的温度向量;求解器求解完成并存储当前模型各个节点的温度向量。
9.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S8具体实现方法为:
S81、整个求解域的温度场计算完成后获得温度场向量,遍历单元计算求解温度梯度矩阵{g},计算公式如下,
式中:Ti、Tj、Tk、Tl为单元节点温度,中间矩阵为梯度矩阵[B]矩阵,Ni、Nj、Nk、Nl为单元形函数;
S82、计算单元热流梯度{q}e向量,整个求解域的温度梯度计算完成后,计算单元热流,热传导满足傅里叶定律,公式如下,
式中:qx,qy,qz为x,y,z方向的热通量;k是热导率,dT/dx,dT/dy,dT/dz是温度梯度。
10.根据权利要求1所述的一种高效求解非线性热边界问题的改良弧长法,其特征在于,所述步骤S9具体实现方法为:
S91、根据单元热通量{q}e计算节点平均热流,平均热流矢量计算公式为:
式中:为节点平均热流矢量;{q}e为包含该节点单元的单元热流大小,n为包含该节点的单元个数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202410343630.2A CN118410663A (zh) | 2024-03-25 | 2024-03-25 | 一种高效求解非线性热边界问题的改良弧长法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202410343630.2A CN118410663A (zh) | 2024-03-25 | 2024-03-25 | 一种高效求解非线性热边界问题的改良弧长法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN118410663A true CN118410663A (zh) | 2024-07-30 |
Family
ID=91990095
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202410343630.2A Pending CN118410663A (zh) | 2024-03-25 | 2024-03-25 | 一种高效求解非线性热边界问题的改良弧长法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN118410663A (zh) |
-
2024
- 2024-03-25 CN CN202410343630.2A patent/CN118410663A/zh active Pending
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Fan et al. | Locally optimal reach set over-approximation for nonlinear systems | |
CN112347687A (zh) | 一种自适应自由度电磁-温度多物理场耦合分析方法 | |
CN112906272B (zh) | 一种反应堆稳态物理热工全耦合精细数值模拟方法及系统 | |
CN113177290B (zh) | 基于深度代理模型归一化的卫星组件温度场预测方法 | |
CN112580855A (zh) | 基于自适应变异pso-bp神经网络的电缆群稳态温升预测方法 | |
Roman et al. | Fast eigenvalue calculations in a massively parallel plasma turbulence code | |
Zhuang et al. | Temperature-constrained topology optimization of nonlinear heat conduction problems | |
Zheng et al. | Concurrent topology optimization for thermoelastic structures with random and interval hybrid uncertainties | |
CN111241728B (zh) | 一种欧拉方程的间断伽辽金有限元数值求解方法 | |
CN110110406A (zh) | 一种基于Excel计算平台实现LS-SVM模型的边坡稳定性预测方法 | |
CN118410663A (zh) | 一种高效求解非线性热边界问题的改良弧长法 | |
CN116090260A (zh) | 一种反应堆全耦合的系统仿真方法 | |
CN115422808A (zh) | 一种基于Krylov子空间的变压器温度场模型降阶方法 | |
CN115270487A (zh) | 一种燃料组件级辐照热流固紧密耦合计算装置及计算方法 | |
CN113868916B (zh) | 基于lfvpso-bpnn的多回路沟槽敷设电缆温升预测方法 | |
Mastrippolito et al. | RBF-based mesh morphing improvement using Schur complement applied to rib shape optimization | |
Moukalled et al. | The discretization process | |
CN114880792A (zh) | 一种基于形变预测的全方位多角度优化方法 | |
Amrit et al. | Efficient multi-objective aerodynamic optimization by design space dimension reduction and co-kriging | |
CN118171519A (zh) | 一种高效求解几何非线性问题的组合法 | |
Wang et al. | Unsteady Aerodynamic Modeling Based on POD‐ARX | |
CN115700577A (zh) | 载流子输运仿真方法、装置、介质及电子设备 | |
Boschitsch et al. | High accuracy computation of fluid-structure interaction in transonic cascades | |
Ding et al. | A Novel Method for Output Characteristics Calculation of Electromagnetic Devices using Multi-kernel RBF Neural Network | |
CN117874941B (zh) | 三维超导管内电缆导体冷却和失超模拟方法、设备、介质 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |