CN108932392B - 基于改进三重互易边界元法的瞬态温度计算方法 - Google Patents

基于改进三重互易边界元法的瞬态温度计算方法 Download PDF

Info

Publication number
CN108932392B
CN108932392B CN201810768920.6A CN201810768920A CN108932392B CN 108932392 B CN108932392 B CN 108932392B CN 201810768920 A CN201810768920 A CN 201810768920A CN 108932392 B CN108932392 B CN 108932392B
Authority
CN
China
Prior art keywords
equation
boundary
domain
time
frequency domain
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.)
Active
Application number
CN201810768920.6A
Other languages
English (en)
Other versions
CN108932392A (zh
Inventor
郭帅平
黄伟
宾光富
李学军
杨大炼
冯和英
廖子豪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan University of Science and Technology
Original Assignee
Hunan University of Science and Technology
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 Hunan University of Science and Technology filed Critical Hunan University of Science and Technology
Priority to CN201810768920.6A priority Critical patent/CN108932392B/zh
Publication of CN108932392A publication Critical patent/CN108932392A/zh
Application granted granted Critical
Publication of CN108932392B publication Critical patent/CN108932392B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E60/00—Enabling technologies; Technologies with a potential or indirect contribution to GHG emissions mitigation

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)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于改进三重互易边界元法的瞬态温度计算方法，包括以下步骤：采用倍数时刻方案确定计算时刻序列；基于Laplace变换得到瞬态热传导问题的频域方程，并基于分析时刻t和Laplace数值逆变换方法确定频域方程的频域参数；基于修正亥姆霍兹基本解得到频域方程的边界‑域积分方程；基于改进三重互易公式将边界‑域积分方程包含的域积分转化为等效边界积分，得到纯边界积分方程；基于边界元法理论表面网格与域内插值点，求解频域边界积分方程；对频域解进行Laplace逆变换，求解时域温度分布。本发明提出了低自由度表达的改进三重互易法计算公式，大大地减少了三重互易法的计算时间，同时节约了存储空间。

Description

基于改进三重互易边界元法的瞬态温度计算方法
技术领域
本发明涉及一种温度计算方法，特别涉及一种基于改进三重互易边界元法的瞬态温度计算方法。
背景技术
三维瞬态热传导问题广泛存在于国民经济各个领域，如核电、化工、建筑、热成型、热加工等，瞬态温度的仿真分析对工业生产中各环节温度的了解与控制至关重要。边界元法是进行热传导分析的常用数值方法，针对含非零初始条件或热源的热传导问题，其边界积分方程中存在域积分，这削减了边界元法降维的优势；针对边界积分方程中存在的域积分，三重互易方法被提出，并被成功地应用于热传导、弹性静力等问题中，但原三重互易法大大增加了求解需要的内存空间和计算时间，计算效率仍有待提高；同时针对瞬态热传导问题，分析时刻方案仍需进一步优化。
发明内容
为了解决上述技术问题，本发明提供一种算法简单的基于改进三重互易边界元法的瞬态温度计算方法。
本发明解决上述问题的技术方案是：一种基于改进三重互易边界元法的瞬态温度计算方法，包括以下步骤：
S1：基于计算时长TT和最小时间步长Δt，采用倍数时刻方案确定计算时刻序列；
S2：基于Laplace变换得到瞬态热传导问题的频域方程，并基于分析时刻t和Laplace数值逆变换方法确定频域方程的频域参数；
S3：基于修正亥姆霍兹基本解得到频域方程的边界-域积分方程；
S4：基于改进三重互易公式将边界-域积分方程包含的域积分转化为等效边界积分，得到纯边界积分方程；
S5：基于边界元法理论表面网格与域内插值点，求解频域边界积分方程；
S6：对频域边界积分方程的解进行Laplace逆变换，求解时域温度分布，得到计算时刻序列中所有时刻的温度值；
S7：判断所得所有温度值中相邻时刻的最大温度变化值是否小于Δumax，Δumax为相邻时刻间允许最大温度变化值，若否，则返回步骤S1，基于倍数时刻方案进行新的计算时刻***，重复S2-S6步骤计算新***时刻的温度值，直至满足所有相邻时刻的最大温度变化值小于Δumax。
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S1具体步骤为：
S11：假设初始时刻为0，针对时刻ti，采用基于Laplace变换的频域法进行求解，计算ti时刻的频域参数sij为
Figure BDA0001729754870000021
/>
β＝log2，整数NL表示逆变换精度参数，一个时刻ti对应2NL个频域参数sij，每一个频域参数sij对应一个需要求解的频域方程，针对时刻ti与tl＝2ti，根据上式可知，时刻ti与tl的频域参数sij与slj有一半重合，考虑到成2倍关系的时刻中tl时刻可节约一半的计算量，提出倍数时刻方案，即取计算时刻
ti＝2i·Δt i＝0,1,2,......,ti≤TT
式中Δt表示最小时间步长。
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S2具体步骤为：
瞬态热传导问题的控制方程为：
Figure BDA0001729754870000031
式中a为热扩散系数，a＝k/(cρ)，k表示材料的热传导系数，c和ρ分别表示材料的比热和密度；Ω表示模型所占的几何区域，x为Ω中的任意点，称为场点，u(x,t)表示在t时刻点x处的温度，w(x,t)表示热源密度函数，对方程两边均进行Laplace变换，得到
Figure BDA0001729754870000032
式中
Figure BDA0001729754870000033
分别表示Laplace变换后的温度函数、热源密度函数，u0(x)为初始温度分布函数，s为频域参数。
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S3中，运用边界积分方程方法将频域方程转化为弱形式的积分方程，即：
Figure BDA0001729754870000034
其中u*(x,y)为修正亥姆霍兹基本解，且
Figure BDA0001729754870000035
r＝r(x,y)为源点y到场点x的距离，λ＝s/a；
运用高斯散度定理，得到边界-域积分方程：
Figure BDA0001729754870000036
其中
Figure BDA0001729754870000037
为源点y的温度，Γ为Ω的边界；/>
Figure BDA0001729754870000038
为热流密度函数，n为边界Γ上的外法线方向；/>
Figure BDA0001729754870000039
为热流密度基本解；c(y)为关于源点y位置的函数，当源点y∈Γ时，c(y)＝0.5,当源点y∈Ω时，c(y)＝1。
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S4具体步骤如下：
S41：针对瞬态热传导问题中热源及初始温度分布引起边界积分方程中的域积分，基于三重互易法，记域积分函数为：
Figure BDA0001729754870000041
采用满足如下方程的b1(x)对域函数b(x)进行近似，
b1(xnm)＝b(xnm)
Figure BDA00017297548700000410
/>
Figure BDA0001729754870000042
其中δ(x,xm)为狄拉克δ函数，xnm为所有边界节点和域内插值点，nm＝1,2,3,…,(N+M)，其中N为边界节点的个数，M为域内插值点的个数，xm表示域内插值点，m＝1,2,3…M，b1(x)、b2(x)、b3(xm)为待求的未知函数；
S42：对S41中的近似函数b1(x)满足的方程，基于Laplace基本解及边界元法进行求解b1(x)、b2(x)、b3(xm)，即得到如下边界积分方程：
Figure BDA0001729754870000043
Figure BDA0001729754870000044
Figure BDA0001729754870000045
其中
Figure BDA0001729754870000046
为Laplace方程基本解，/>
Figure BDA0001729754870000047
为高阶基本解，满足
Figure BDA0001729754870000048
对于θ＝1和2有/>
Figure BDA0001729754870000049
ym为域内源点，从域内插值点xm中依次选取；
S43：对S42中的边界积分方程，对第一、第二个方程中的源点y依次选N个边界节点，并对第三方程中的源点y依次选M个域内插值点，将共得到(2N+M)个方程，并基于边界元法理论，得到如下(2N+M)维线性代数方程组
G1D1+H2B2-G2D2-G2DB3D＝H1B1
H1B2-G1D2-G1DB3D＝0
Figure BDA0001729754870000051
G1、G2、H1、H2为方程系数矩阵，为N维方阵，其第
Figure BDA0001729754870000052
行ψ1列元素分别为：
Figure BDA0001729754870000053
Figure BDA0001729754870000054
其中
Figure BDA0001729754870000055
为第ψ1个边界节点/>
Figure BDA0001729754870000056
的形函数；系数矩阵/>
Figure BDA0001729754870000057
为(M×N)维矩阵，其第/>
Figure BDA0001729754870000058
行ψ2列元素为：
Figure BDA0001729754870000059
Figure BDA00017297548700000510
G1D、G2D为(N×M)维系数矩阵，其第
Figure BDA00017297548700000511
行ψ3列元素分别为：
Figure BDA00017297548700000512
/>
Figure BDA00017297548700000513
为M维系数方矩，其第/>
Figure BDA00017297548700000514
行ψ4列元素分别为：
Figure BDA00017297548700000515
其中B1、B2、D1、D2为N维列向量，分别为函数b1(x)、b2(x)、d1(x)、d2(x)在N个边界节点上的取值，其中B1、B1D为已知；M维列向量B1D、B3D为函数b1(x)和b3(x)在M个域内插值点上的取值；
令B2＝0，得到求解总方程组如下
Figure BDA00017297548700000516
S44：对S43中(2N+M)维的线性代数方程组，基于自由度凝集方法进行降维求解，先记右端向量为：
H1B1＝R1,
Figure BDA00017297548700000517
将方程组拆散重组，进行自由度凝集，得低维矩阵表达形式：
Figure BDA0001729754870000061
进一步推导另外两个未知向量的表达式，并简记为：
B3D＝TB1D+Τ1B1
D2＝Τ2B3D
D1＝Τ3B1+Τ4D2+Τ5B3D
其中
Figure BDA0001729754870000062
Figure BDA0001729754870000063
Figure BDA0001729754870000064
Figure BDA0001729754870000065
Figure BDA0001729754870000066
Figure BDA0001729754870000067
Figure BDA0001729754870000068
通过上述降维的表达式，求解得到b1(x)高阶函数d1(x)、d2(x)、b3(x)在边界节点及域内插值点上的取值向量D1、D2、B3D。
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S5具体步骤如下：
S51：将S44中得到的近似函数b1(x)代替频域边界积分方程中的域函数b(x)，并基于互易理论得到不含域内积分的纯边界积分方程；
Figure BDA0001729754870000069
/>
S52：选择所有边界节点作为源点，基于边界元法理论进行边界单元内的积分运算，得到如下矩阵方程
HU-GQ＝R
其中矩阵H、G为系数矩阵；U为温度
Figure BDA00017297548700000610
在所有边界节点上的取值组成的向量矩阵，记为/>
Figure BDA0001729754870000071
Q表示热流密度/>
Figure BDA0001729754870000072
在所有边界节点上的取值组成的向量矩阵，记为
Figure BDA0001729754870000073
右端R为域积分对应的向量：
Figure BDA0001729754870000074
其中矩阵H、G第
Figure BDA0001729754870000075
行ψ1列元素的为：
Figure BDA0001729754870000076
Figure BDA0001729754870000077
Figure BDA0001729754870000078
为狄拉克δ函数；矩阵GD第/>
Figure BDA0001729754870000079
行ψ4列元素为：
Figure BDA00017297548700000710
求解上述线性代数方程组，即得到频域方程的解
Figure BDA00017297548700000711
上述基于改进三重互易边界元法的瞬态温度计算方法，所述步骤S6具体步骤如下：
S61：采用只包含实数域运算的Gaver functionals方法进行Laplace逆变换，基于该方法及计算时间ti确定频域参数sij及待求频域方程，基于步骤S3-S5得到2NL个频域解
Figure BDA00017297548700000712
进行Laplace逆变换得到时域解u(ti)，逆变换公式为：
Figure BDA00017297548700000713
Figure BDA00017297548700000714
Figure BDA00017297548700000715
其中
Figure BDA00017297548700000716
为逆变换递推数组第Λ行第ν列元素，uε(ti)为时域解；
S62：对S61中得到的时域解采用Wynn’s rho algorithm进行外推，得到更高精度的时域解，即
Figure BDA00017297548700000717
Figure BDA00017297548700000718
Figure BDA00017297548700000719
Figure BDA00017297548700000720
为递推数组第/>
Figure BDA00017297548700000721
行第τ列元素，最终得到ti时刻的更高精度的时域解u(ti)。
本发明的有益效果在于：
1、本发明的步骤S1中采用倍数时刻方案确定计算时刻序列；依次***计算时刻的目的在于保证得到温度随时间变化规律的同时，减少计算时刻数量，节约计算时间。倍数时刻方案能做到如此灵活地***计算时刻而不受前后时刻影响的原因，是因为频域法计算各时刻温度值时具有相互独立的特点，在已经计算完温度的两个时刻间增加新的计算时刻十分便利，这使得计算效率得到大幅提高。
2、本发明提出了低自由度表达的改进三重互易法计算公式，从而得到近似函数b1(x)及其相关函数在边界节点及域内插值点上的取值D1、D2、B2、B3D；该降维的新求解形式，其计算量为原(2N+M)维矩阵方程计算时间的1/6～1/8，大大地减少了三重互易法的计算时间，同时节约了维度为(2N+M)矩阵的存储空间，极大地提升了边界元法能求解的问题的复杂度。
附图说明
图1为本发明的整体流程图。
图2为本发明的步骤S1中倍数时刻方案的示意图。
图3是本发明一种实施算例的几何模型与三重互易边界元法离散模型。
图4是本发明实施算例中改进TRM计算精度与边界节点个数N及域内插值点个数M的关系图。
图5是本发明实施算例中改进三重互易公式的计算效率与原三重互易公式的比较图。
图6是本发明实施算例中采用倍数时刻方案温度随时间变化的曲线图。
具体实施方式
下面结合附图和实施例对本发明作进一步的说明。
如图1所示，一种基于改进三重互易边界元法的瞬态温度计算方法，包括以下步骤：
S1：基于计算时长TT和最小时间步长Δt，采用倍数时刻方案确定计算时刻序列。具体步骤为：
S11：假设初始时刻为0，针对时刻ti，采用基于Laplace变换的频域法进行求解，计算ti时刻的频域参数sij为
Figure BDA0001729754870000091
β＝log2，整数NL表示逆变换精度参数，一个时刻ti对应2NL个频域参数sij，每一个频域参数sij对应一个需要求解的频域方程，针对时刻ti与tl＝2ti，根据上式可知，时刻ti与tl的频域参数sij与slj有一半重合，考虑到成2倍关系的时刻中tl时刻可节约一半的计算量，提出倍数时刻方案，即取计算时刻
ti＝2i·Δt i＝0,1,2,......,ti≤TT
式中Δt表示最小时间步长。
为减少计算时刻间的间隔导致相邻时刻温度变化过大，依次补充:
tii＝3·2iΔt i＝0,1,2,......,tii≤TT
tiii＝5·2iΔt i＝0,1,2,......,tiii≤TT
得到的初始时刻序列{ti}U{tii}U{tiii}，如图2所示。
对S11中的时间点选择方案，能节约近50％的计算量，且时间轴上的计算时间能够较好地展示温度随时间变化的规律，以Δt为基数，相邻时刻的时间间隔依次为1,1,1,1,1,1,2,2,2,4,4,4,…,2kk,2kk,2kk,2kk+1,2kk+1,2kk+1,……倍Δt，kk＝0,1,2，…。若温度仿真分析要求相邻时刻间允许最大温度变化值为Δumax，经过S2-S6的各时刻温度计算后，若计算时刻序列中存在多个相邻时刻间隔温度变化大于Δumax，可根据实际情况逐步选择继续***时刻序列{7·2iΔt}，{9·2iΔt}，…，{(2nn+1)·2iΔt}，nn＝3,4,5，…，直到仅存在若干相邻时刻间温度变化较大，此时可直接在该间隔中点***计算时刻，如图2所示，以满足相邻时刻间温度变化小于Δumax的要求。
依次***计算时刻的目的在于保证得到温度随时间变化规律的同时，减少计算时刻数量，节约计算时间。倍数时刻方案能做到如此灵活地***计算时刻而不受前后时刻的影响的原因，是因为频域法计算各时刻温度值时具有相互独立的特点，在已经计算完温度的两个时刻间增加新的计算时刻十分便利，这使得计算效率得到大幅提高。
S2：基于Laplace变换得到瞬态热传导问题的频域方程，并基于分析时刻t和Laplace数值逆变换方法确定频域方程的频域参数。具体步骤为：
瞬态热传导问题的控制方程为：
Figure BDA0001729754870000101
式中a为热扩散系数，a＝k/(cρ)，k表示材料的热传导系数，c和ρ分别表示材料的比热和密度；Ω表示模型所占的几何区域，x为Ω中的任意点，称为场点，u(x,t)表示在t时刻点x处的温度，w(x,t)表示热源密度函数，对方程两边均进行Laplace变换，得到
Figure BDA0001729754870000102
式中
Figure BDA0001729754870000103
分别表示Laplace变换后的温度函数、热源密度函数，u0(x)为初始温度分布函数，s为频域参数。
S3：基于修正亥姆霍兹基本解得到频域方程的边界-域积分方程。
运用边界积分方程方法将频域方程转化为弱形式的积分方程，即：
Figure BDA0001729754870000104
其中u*(x,y)为修正亥姆霍兹基本解，且
Figure BDA0001729754870000111
r＝r(x,y)为源点y到场点x的距离，λ＝s/a；
运用高斯散度定理，得到边界-域积分方程：
Figure BDA0001729754870000112
其中
Figure BDA0001729754870000113
为源点y的温度，Γ为Ω的边界；/>
Figure BDA0001729754870000114
为热流密度函数，n为边界Γ上的外法线方向；/>
Figure BDA0001729754870000115
为热流密度基本解；c(y)为关于源点y位置的函数，当源点y∈Γ时，c(y)＝0.5,当源点y∈Ω时，c(y)＝1。
S4：基于改进三重互易公式将边界-域积分方程包含的域积分转化为等效边界积分，得到纯边界积分方程。具体步骤如下：
S41：针对瞬态热传导问题中热源及初始温度分布引起边界积分方程中的域积分，基于三重互易法，记域积分函数为：
Figure BDA0001729754870000116
采用满足如下方程的b1(x)对域函数b(x)进行近似，
b1(xnm)＝b(xnm)
Figure BDA00017297548700001111
Figure BDA0001729754870000117
其中δ(x,xm)为狄拉克δ函数，xnm为所有边界节点和域内插值点，nm＝1,2,3,…,(N+M)，其中N为边界节点的个数，M为域内插值点的个数，xm表示域内插值点，m＝1,2,3…M，b1(x)、b2(x)、b3(xm)为待求的未知函数。
S42：对S41中的近似函数b1(x)满足的方程，基于Laplace基本解及边界元法进行求解b1(x)、b2(x)、b3(xm)，即得到如下边界积分方程：
Figure BDA0001729754870000118
Figure BDA0001729754870000119
Figure BDA00017297548700001110
其中
Figure BDA0001729754870000121
为Laplace方程基本解，/>
Figure BDA0001729754870000122
为高阶基本解，满足
Figure BDA0001729754870000123
对于θ＝1和2有/>
Figure BDA0001729754870000124
ym为域内源点，从域内插值点xm中依次选取。
S43：对S42中的边界积分方程，对第一、第二个方程中的源点y依次选N个边界节点，并对第三方程中的源点y依次选M个域内插值点，将共得到(2N+M)个方程，并基于边界元法理论，得到如下(2N+M)维线性代数方程组
G1D1+H2B2-G2D2-G2DB3D＝H1B1
H1B2-G1D2-G1DB3D＝0
Figure BDA0001729754870000125
G1、G2、H1、H2为方程系数矩阵，为N维方阵，其第
Figure BDA0001729754870000126
行ψ1列元素分别为：
Figure BDA0001729754870000127
Figure BDA0001729754870000128
其中
Figure BDA0001729754870000129
为第ψ1个边界节点/>
Figure BDA00017297548700001210
的形函数；系数矩阵/>
Figure BDA00017297548700001211
为(M×N)维矩阵，其第/>
Figure BDA00017297548700001212
行ψ2列元素为：
Figure BDA00017297548700001213
/>
Figure BDA00017297548700001214
G1D、G2D为(N×M)维系数矩阵，其第
Figure BDA00017297548700001215
行ψ3列元素分别为：
Figure BDA00017297548700001216
Figure BDA00017297548700001217
为M维系数方矩，其第/>
Figure BDA00017297548700001218
行ψ4列元素分别为：
Figure BDA00017297548700001219
其中B1、B2、D1、D2为N维列向量，分别为函数b1(x)、b2(x)、d1(x)、d2(x)在N个边界节点上的取值，其中B1、B1D为已知；M维列向量B1D、B3D为函数b1(x)和b3(x)在M个域内插值点上的取值；
令B2＝0，得到求解总方程组如下
Figure BDA0001729754870000131
上式为传统的三重互易法求解公式，针对复杂工程问题，变量的自由度N及M较大，上式方程的维度将进一步增大，除了频域方程的求解以外，基于三重互易法的近似函数求解需要额外的存储空间和巨大的计算时间。
S44：针对上述问题，本发明根据上述方程组的特点，基于自由度凝集方法，将方程进行转化，先记右端向量为：
H1B1＝R1,
Figure BDA0001729754870000132
将方程组拆散重组，进行自由度凝集，得低维矩阵表达形式：
Figure BDA0001729754870000133
进一步推导另外两个未知向量的表达式，并简记为：
B3D＝TB1D+Τ1B1
D2＝Τ2B3D
D1＝Τ3B1+Τ4D2+Τ5B3D
其中
Figure BDA0001729754870000134
Figure BDA0001729754870000135
Figure BDA0001729754870000136
Figure BDA0001729754870000137
Figure BDA0001729754870000138
Figure BDA0001729754870000139
Figure BDA00017297548700001310
通过上述降维的表达式，求解得到b1(x)高阶函数d1(x)、d2(x)、b3(x)在边界节点及域内插值点上的取值向量D1、D2、B3D。该降维的新求解形式，其计算量约为S43中矩阵求解形式的1/6～1/8，极大减少了计算时间，同时节约了维度为(2N+M)的矩阵的存储空间，极大地提升了边界元法能求解的问题的复杂度。
S5：基于边界元法理论表面网格与域内插值点，求解频域边界积分方程。具体步骤如下：
S51：将S44中得到的近似函数b1(x)代替频域边界积分方程中的域函数b(x)，并基于互易理论得到不含域内积分的纯边界积分方程；
Figure BDA0001729754870000141
S52：基于边界元理论对模型边界进行离散，得到边界单元与节点，同时在Ω内仅生成插值点；分别取所有边界节点作为源点建立方程组，基于边界元法理论进行边界单元内的积分运算，得到如下矩阵方程
HU-GQ＝R
其中矩阵H、G为系数矩阵；U为温度
Figure BDA0001729754870000142
在所有边界节点上的取值组成的向量矩阵，记为/>
Figure BDA0001729754870000143
Q表示热流密度/>
Figure BDA0001729754870000144
在所有边界节点上的取值组成的向量矩阵，记为
Figure BDA0001729754870000145
右端R为域积分对应的向量：
Figure BDA0001729754870000146
其中矩阵H、G第
Figure BDA0001729754870000147
行ψ1列元素的为：
Figure BDA0001729754870000148
Figure BDA0001729754870000149
Figure BDA00017297548700001410
为狄拉克δ函数；矩阵GD第/>
Figure BDA00017297548700001411
行ψ4列元素为：
Figure BDA00017297548700001412
求解上述线性代数方程组，即得到频域方程的解
Figure BDA00017297548700001413
S6：对频域解进行Laplace逆变换，求解时域温度分布，得到计算时刻序列中所有时刻的温度值。具体步骤如下：
S61：上面描述了如何在不进行域内网格划分的前提下进行频域方程快速求解，为了得到时域解，应用数值Laplace逆变换对频域解进行转化。本发明采用只包含实数域运算的Gaver functionals方法进行Laplace逆变换，基于该方法及计算时间ti确定频域参数sij及待求频域方程，基于步骤S3-S5得到2NL个频域解
Figure BDA0001729754870000151
进行Laplace逆变换得到时域解u(ti)，逆变换公式为：
Figure BDA0001729754870000152
Figure BDA0001729754870000153
/>
Figure BDA0001729754870000154
其中
Figure BDA0001729754870000155
为逆变换递推数组第Λ行第ν列元素，uε(ti)为时域解，其下标ε为自然数，表示精度；ε越大精度越高，ε取值从小到大，得到一个精度从低到高的时域解序列uε(ti)；
S62：对S61中得到的时域解采用Wynn’s rho algorithm进行外推，得到更高精度的时域解，即
Figure BDA0001729754870000156
Figure BDA0001729754870000157
Figure BDA0001729754870000158
Figure BDA0001729754870000159
为递推数组第/>
Figure BDA00017297548700001510
行第τ列元素，最终得到ti时刻的更高精度的时域解u(ti)。
S7：判断所得所有温度值中相邻时刻的最大温度变化值是否小于Δumax，Δumax为相邻时刻间允许最大温度变化值，若否，则返回步骤S1，基于倍数时刻方案进行新的计算时刻***，重复S2-S6步骤计算新***时刻的温度值，直至满足所有相邻时刻的最大温度变化值小于Δumax。
实施例
考虑如图3所示几何模型，其瞬态温度解析场为
u(x,y,z,t)＝x6+30x4t+180x2t2+y2+z+120t3+2t
给定的初始温度场为
u0＝x6+y2+z
热扩散系数a＝1m2/s，分析中给定纽曼边界条件，求解瞬态温度分布。分析中计算得到的温度值与解析解进行对比，误差计算公式如下：
Figure BDA0001729754870000161
其中
Figure BDA0001729754870000162
表示采用本文方法得到的数值解，/>
Figure BDA0001729754870000163
为问题的解析解。
给定分析时间范围为10s，最小时间步长为0.2s，基于倍数时刻方案，得到待求的初步计算时刻序列为{0.2,0.4,0.6,0.8,1.0,1.2,1.6,2.0,2.4,3.2,4.0,4.8,6.4,9.6}，单位为秒。
针对上述计算时刻进行边界元分析，基于Laplace变换的三重互易边界元法计算精度与三个参数有关，第一个为逆变换参数NL，第二个参数为边界节点个数N，第三个参数为域内插值点的个数M；其中参数NL直接决定逆变换的计算精度，在本次仿真中取NL＝8保证逆变换的精度，便于考虑另外两个参数对仿真结果的影响；参数N决定了未知函数的插值精度和域积分函数插值精度，参数M决定了域积分函数插值精度，通常情况下N和M越大计算精度越高，域内插值点的个数可由域内函数的变化梯度来决定。
本算例中边界单元采用二次四边形单元，图4中左图显示了计算误差随边界节点个数N的变化趋势，如图4中右图显示了计算误差随域内插值点个数M的变化趋势，从图中可以看出改进TRM的计算误差随着边界节点和域内插值点个数的增加而减小，显示了改进TRM计算的稳定性，同时在本次仿真分析中，当边界节点个数为2407、域内插值点个数为1475时，计算结果已趋于稳定。
为了进一步显示改进TRM算法在计算速度和存储空间上的优势，研究了TRM插值的计算时间随边界节点N增加的变化趋势，如图5所示，圆形图标为采用原TRM公式下的计算时间曲线，正方形图标为采用改进TRM公式下的计算时间曲线，从图中可以看出改进TRM公式的计算时间显著下降；三角形图标为两者的时间比，其坐标轴为右侧的坐标轴，从图中可以看出改进TRM的计算时间大约为原TRM方法的1/8；其次，原TRM由于需要组装大矩阵进行求解，将受到存储空间的限制，在N＝11000时已经出现内存不足，无法计算，而改进TRM继续增大N，在N＝18000时仍然可以计算，说明改进TRM公式大大减低了TRM插值所需要的矩阵存储空间，可计算的自由度增加了近一倍。总之，改进TRM大大减少了TRM插值的时间，同时节约了大量的存储空间，为计算更大自由度的问题提供了解决方案，且节约的时间和空间随N的增加更加明显。
图6显示了点(-1.8,0.2,0)的温度随时间变化历程，图中黑色正方形点为新的时间步长方案下的结果，计算的时间序列为{0.2,0.4,0.6,0.8,1.0,1.2,1.6,2.0,2.4,3.2,4.0,4.8,6.4,9.6}，共有14个时刻，14个计算时刻原本需要计算224个频域方程，由于倍数时刻方案考虑了逆变换的特点，此处的14个计算时刻只计算了112个频域方程，节约了50％的计算量；由于该算例中温度变化较大，依次***时刻序列{7·2iΔt}，{9·2iΔt}，…，{(2nn+1)·2iΔt}，nn＝3,4,5，…，最终达到等距时间间隔，如图中三角形所示，共有50个时刻，原本需要计算800个频域方程，考虑了逆变换的特点后只计算了496个频域方程，节约了38％的计算量；从图中可以看出，初始倍数时刻方案得到的温度时间曲线展示出较好的光滑性，且仅用了50％的计算量，即使逐步***新的计算时刻直至等时间间隔，仍然节约了大量的频域方程计算时间，计算量为原时间方案的62％。

Claims (4)

1.一种基于改进三重互易边界元法的瞬态温度计算方法，包括以下步骤：
S1：基于计算时长TT和最小时间步长Δt，采用倍数时刻方案确定计算时刻序列；
所述步骤S1具体步骤为：
S11：假设初始时刻为0，针对时刻ti，采用基于Laplace变换的频域法进行求解，计算ti时刻的频域参数sij为
Figure FDA0004190974350000011
β＝log2，整数NL表示逆变换精度参数，一个时刻ti对应2NL个频域参数sij，每一个频域参数sij对应一个需要求解的频域方程，针对时刻ti与tl＝2ti，根据上式可知，时刻ti与tl的频域参数sij与slj有一半重合，考虑到成2倍关系的时刻中tl时刻可节约一半的计算量，提出倍数时刻方案，即取计算时刻
ti＝2i·Δt i＝0,1,2,......；ti≤TT
式中Δt表示最小时间步长；
S2：基于Laplace变换得到瞬态热传导问题的频域方程，并基于分析时刻t和Laplace数值逆变换方法确定频域方程的频域参数；
所述步骤S2具体步骤为：
瞬态热传导问题的控制方程为：
Figure FDA0004190974350000012
式中a为热扩散系数，a＝k/(cρ)，k表示材料的热传导系数，c和ρ分别表示材料的比热和密度；Ω表示模型所占的几何区域，x为Ω中的任意点，称为场点，u(x,t)表示在t时刻点x处的温度，w(x,t)表示热源密度函数，对方程两边均进行Laplace变换，得到
Figure FDA0004190974350000021
式中
Figure FDA0004190974350000022
分别表示Laplace变换后的温度函数、热源密度函数，u0(x)为初始温度分布函数，s为频域参数；
S3：基于修正亥姆霍兹基本解得到频域方程的边界-域积分方程；
S4：基于改进三重互易公式将边界-域积分方程包含的域积分转化为等效边界积分，得到纯边界积分方程；
所述步骤S4具体步骤如下：
S41：针对瞬态热传导问题中热源及初始温度分布引起边界积分方程中的域积分，基于三重互易法，记域积分函数为：
Figure FDA0004190974350000023
采用满足如下方程的b1(x)对域函数b(x)进行近似，
b1(xnm)＝b(xnm)
Figure FDA0004190974350000024
Figure FDA0004190974350000025
其中δ(x,xm)为狄拉克δ函数，xnm为所有边界节点和域内插值点，nm＝1,2,3,…,N+M，其中N为边界节点的个数，M为域内插值点的个数，xm表示域内插值点，m＝1,2,3…M，b1(x)、b2(x)、b3(xm)为待求的未知函数；
S42：对S41中的近似函数b1(x)满足的方程，基于Laplace基本解及边界元法进行求解b1(x)、b2(x)、b3(xm)，即得到如下边界积分方程：
Figure FDA0004190974350000026
Figure FDA0004190974350000027
Figure FDA0004190974350000028
其中
Figure FDA0004190974350000029
为Laplace方程基本解，/>
Figure FDA00041909743500000210
为高阶基本解，满足/>
Figure FDA0004190974350000031
对于θ＝1和2有/>
Figure FDA0004190974350000032
ym为域内源点，从域内插值点xm中依次选取；
S43：对S42中的边界积分方程，对第一、第二个方程中的源点y依次选N个边界节点，并对第三方程中的源点y依次选M个域内插值点，将共得到2N+M个方程，并基于边界元法理论，得到如下2N+M维线性代数方程组
G1D1+H2B2-G2D2-G2DB3D＝H1B1
H1B2-G1D2-G1DB3D＝0
Figure FDA0004190974350000033
G1、G2、H1、H2为方程系数矩阵，为N维方阵，其第
Figure FDA0004190974350000034
行ψ1列元素分别为：
Figure FDA0004190974350000035
Figure FDA0004190974350000036
其中
Figure FDA0004190974350000037
为第ψ1个边界节点/>
Figure FDA0004190974350000038
的形函数；系数矩阵/>
Figure FDA0004190974350000039
为M×N维矩阵，其第/>
Figure FDA00041909743500000310
行ψ2列元素为：
Figure FDA00041909743500000311
Figure FDA00041909743500000312
G1D、G2D为N×M维系数矩阵，其第
Figure FDA00041909743500000313
行ψ3列元素分别为：
Figure FDA00041909743500000314
Figure FDA00041909743500000315
为M维系数方矩，其第/>
Figure FDA00041909743500000316
行ψ4列元素分别为：
Figure FDA00041909743500000317
其中B1、B2、D1、D2为N维列向量，分别为函数b1(x)、b2(x)、d1(x)、d2(x)在N个边界节点上的取值，其中B1、B1D为已知；M维列向量B1D、B3D为函数b1(x)和b3(x)在M个域内插值点上的取值；
令B2＝0，得到求解总方程组如下
Figure FDA0004190974350000041
S44：对S43中2N+M维的线性代数方程组，基于自由度凝集方法进行降维求解，先记右端向量为：
Figure FDA0004190974350000042
将方程组拆散重组，进行自由度凝集，得低维矩阵表达形式：
Figure FDA0004190974350000043
进一步推导另外两个未知向量的表达式，并简记为：
B3D＝TB1D+Τ1B1
D2＝Τ2B3D
D1＝Τ3B1+Τ4D2+Τ5B3D
其中
Figure FDA0004190974350000044
Figure FDA0004190974350000045
Figure FDA0004190974350000046
Figure FDA0004190974350000047
Figure FDA0004190974350000048
Figure FDA0004190974350000049
Figure FDA00041909743500000410
通过上述降维的表达式，求解得到b1(x)高阶函数d1(x)、d2(x)、b3(x)在边界节点及域内插值点上的取值向量D1、D2、B3D；
S5：基于边界元法理论表面网格与域内插值点，求解频域边界积分方程；
S6：对频域边界积分方程的解进行Laplace逆变换，求解时域温度分布，得到计算时刻序列中所有时刻的温度值；
S7：判断所得所有温度值中相邻时刻的最大温度变化值是否小于Δumax，Δumax为相邻时刻间允许最大温度变化值，若否，则返回步骤S1，基于倍数时刻方案进行新的计算时刻***，重复S2-S6步骤计算新***时刻的温度值，直至满足所有相邻时刻的最大温度变化值小于Δumax。
2.根据权利要求1所述的基于改进三重互易边界元法的瞬态温度计算方法，其特征在于，所述步骤S3中，运用边界积分方程方法将频域方程转化为弱形式的积分方程，即：
Figure FDA0004190974350000051
/>
其中u*(x,y)为修正亥姆霍兹基本解，且
Figure FDA0004190974350000052
r＝r(x,y)为源点y到场点x的距离，λ＝s/a；
运用高斯散度定理，得到边界-域积分方程：
Figure FDA0004190974350000053
其中
Figure FDA0004190974350000054
为源点y的温度，Γ为Ω的边界；/>
Figure FDA0004190974350000055
为热流密度函数，n为边界Γ上的外法线方向；/>
Figure FDA0004190974350000056
为热流密度基本解；c(y)为关于源点y位置的函数，当源点y∈Γ时，c(y)＝0.5,当源点y∈Ω时，c(y)＝1。
3.根据权利要求2所述的基于改进三重互易边界元法的瞬态温度计算方法，其特征在于，所述步骤S5具体步骤如下：
S51：将S44中得到的近似函数b1(x)代替频域边界积分方程中的域函数b(x)，并基于互易理论得到不含域内积分的纯边界积分方程；
Figure FDA0004190974350000057
S52：选择所有边界节点作为源点，基于边界元法理论进行边界单元内的积分运算，得到如下矩阵方程
HU-GQ＝R
其中矩阵H、G为系数矩阵；U为温度
Figure FDA0004190974350000058
在所有边界节点上的取值组成的向量矩阵，记为
Figure FDA0004190974350000061
Q表示热流密度/>
Figure FDA0004190974350000062
在所有边界节点上的取值组成的向量矩阵，记为
Figure FDA0004190974350000063
右端R为域积分对应的向量：
Figure FDA0004190974350000064
其中矩阵H、G第
Figure FDA0004190974350000065
行ψ1列元素的为：
Figure FDA0004190974350000066
Figure FDA0004190974350000067
Figure FDA0004190974350000068
为狄拉克δ函数；矩阵GD第/>
Figure FDA0004190974350000069
行ψ4列元素为：
Figure FDA00041909743500000610
求解上述线性代数方程组，即得到频域方程的解
Figure FDA00041909743500000611
4.根据权利要求3所述的基于改进三重互易边界元法的瞬态温度计算方法，其特征在于，所述步骤S6具体步骤如下：
S61：采用只包含实数域运算的Gaver functionals方法进行Laplace逆变换，基于该方法及计算时间ti确定频域参数sij及待求频域方程，基于步骤S3-S5得到2NL个频域解
Figure FDA00041909743500000612
进行Laplace逆变换得到时域解u(ti)，逆变换公式为：/>
Figure FDA00041909743500000613
Figure FDA00041909743500000614
Figure FDA00041909743500000615
其中
Figure FDA00041909743500000616
为逆变换递推数组第Λ行第ν列元素，uε(ti)为时域解；
S62：对S61中得到的时域解采用Wynn’s rho algorithm进行外推，得到更高精度的时域解，即
Figure FDA00041909743500000617
Figure FDA00041909743500000618
Figure FDA00041909743500000619
Figure FDA00041909743500000620
为递推数组第/>
Figure FDA00041909743500000621
行第τ列元素，最终得到ti时刻的更高精度的时域解u(ti)。/>
CN201810768920.6A 2018-07-13 2018-07-13 基于改进三重互易边界元法的瞬态温度计算方法 Active CN108932392B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810768920.6A CN108932392B (zh) 2018-07-13 2018-07-13 基于改进三重互易边界元法的瞬态温度计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810768920.6A CN108932392B (zh) 2018-07-13 2018-07-13 基于改进三重互易边界元法的瞬态温度计算方法

Publications (2)

Publication Number Publication Date
CN108932392A CN108932392A (zh) 2018-12-04
CN108932392B true CN108932392B (zh) 2023-05-26

Family

ID=64447651

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810768920.6A Active CN108932392B (zh) 2018-07-13 2018-07-13 基于改进三重互易边界元法的瞬态温度计算方法

Country Status (1)

Country Link
CN (1) CN108932392B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109829130B (zh) * 2019-01-29 2023-06-16 武汉轻工大学 二重积分计算方法、装置、终端设备及可读存储介质
CN110069837B (zh) * 2019-04-04 2023-06-09 天津职业技术师范大学(中国职业培训指导教师进修中心) 横观各向同性多层涂层体系三维温度场的求解方法
CN110032787B (zh) * 2019-04-04 2023-06-09 天津职业技术师范大学(中国职业培训指导教师进修中心) 各向同性多层涂层体系二维温度场的求解方法
CN110070591B (zh) * 2019-04-25 2023-01-20 湖南科技大学 一种计算机绘图用的多边形填充方法
CN110162849B (zh) * 2019-05-07 2022-08-16 南京理工大学 一种混杂纤维混凝土的建模方法
CN111597745B (zh) * 2020-05-14 2024-05-24 北京工业大学 一种用于求解瞬态热传导问题的时空数值方法
CN112013978B (zh) * 2020-09-03 2022-04-08 安徽大学 一种温度传感器动态温度测量的自动补偿方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101833597A (zh) * 2010-04-08 2010-09-15 西北工业大学 用快速边界元法得到大型复杂飞行器电场分布的方法
CN105760345A (zh) * 2016-02-02 2016-07-13 河海大学 扩散型动态数据重构的奇异边界法
CN107391895A (zh) * 2017-09-15 2017-11-24 湖南科技大学 基于转移矩阵的轴对称热声谐振管频率计算方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101833597A (zh) * 2010-04-08 2010-09-15 西北工业大学 用快速边界元法得到大型复杂飞行器电场分布的方法
CN105760345A (zh) * 2016-02-02 2016-07-13 河海大学 扩散型动态数据重构的奇异边界法
CN107391895A (zh) * 2017-09-15 2017-11-24 湖南科技大学 基于转移矩阵的轴对称热声谐振管频率计算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
金属热处理瞬态温度场的边界元法研究;杨书申;《郑州纺织工学院学报》;20000905(第03期);全文 *

Also Published As

Publication number Publication date
CN108932392A (zh) 2018-12-04

Similar Documents

Publication Publication Date Title
CN108932392B (zh) 基于改进三重互易边界元法的瞬态温度计算方法
Holten et al. Comparison of different methods for state estimation
Tan Recursion-transform method for computing resistance of the complex resistor network with three arbitrary boundaries
Gofuku et al. Point-tangent/point-normal B-spline curve interpolation by geometric algorithms
Leonard et al. The NIRVANA scheme applied to one‐dimensional advection
Ha et al. Topological shape optimization of heat conduction problems using level set approach
Swarztrauber et al. Generalized discrete spherical harmonic transforms
Cai et al. Numerical regularized moment method of arbitrary order for Boltzmann-BGK equation
Fox et al. Hyperbolic quadrature method of moments for the one-dimensional kinetic equation
Gherardini et al. Reconstructing quantum entropy production to probe irreversibility and correlations
Abinandanan et al. An extended Cahn-Hilliard model for interfaces with cubic anisotropy
Frolkovič et al. High-resolution flux-based level set method
Choi et al. Sampling-based RBDO of ship hull structures considering thermo-elasto-plastic residual deformation
Chen The differential quadrature element method irregular element torsion analysis model
Kindelan et al. Radial basis function interpolation in the limit of increasingly flat basis functions
Isik et al. Bernstein series solution of a class of Lane‐Emden type equations
Liu et al. Explicit matrix representation for NURBS curves and surfaces
Azencot et al. Advection‐based function matching on surfaces
Wang et al. The Error Estimates of the Interpolating Element‐Free Galerkin Method for Two‐Point Boundary Value Problems
Zhu et al. A proper orthogonal decomposition analysis method for transient nonlinear heat conduction problems. Part 2: Advanced algorithm
Chen et al. A novel approach to uncertainty analysis using methods of hybrid dimension reduction and improved maximum entropy
Luo et al. An improved finite volume scheme for compressible flows on unstructured grids
Løvgren et al. The reduced basis element method for fluid flows
Hanukah et al. Improving mass matrix and inverse mass matrix computations of hexahedral elements
Firsov et al. FVTD—integral equation hybrid for Maxwell's equations

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
GR01 Patent grant
GR01 Patent grant