CN109753691B - 针对钠冷快堆单排组件耦合热变形的模拟方法 - Google Patents

Abstract

一种针对钠冷快堆单排组件耦合热变形的模拟方法，1、确定单排组件的组件个数、目标组件的物性参数以及几何参数，2、将单排组件按照由堆芯朝堆外的方向用数字编号进行标记，3、获取每一根组件在所有轴向位置的自由热变形位移，4、对每一根组件在上、下垫块处的热变形位移进行取值，5、推导出四种不同的接触力与垫块之间的位移模式，6、列出每一根组件在上、下垫块处的平衡方程，7、列出与所有组件相关的变形协调方程，8、将组件的位移参数进行合并处理，9、获得上、下垫块位置的总位移；10、求解每相邻两根组件之间上、下垫块的接触力，11、通过组件的位移叠加原则求得每一根组件的最终位移；本发明方法具备计算量小，收敛性好、计算结果准确的优点。

Description

针对钠冷快堆单排组件耦合热变形的模拟方法

技术领域

本发明涉及钠冷快堆单排组件耦合热变形技术领域，具体涉及一种针对钠冷快堆单排组件耦合热变形的模拟方法。

背景技术

钠冷快堆中的组件在整个钠冷快堆的运行过程中占有非常重要的地位。在钠冷快堆的运行过程中，燃料组件的活性区中会发生裂变反应释放大量的热量，而燃料组件周围充当冷却剂的钠会通过流动的方式将燃料组件内裂变反应释放的热量带出堆芯外。在我国自主设计研发的钠冷快堆中，由于中子通量分布不均、燃料富集度差异以及堆芯冷井等多重因素的影响，快堆中组件释热功率分布是不均匀的，并且在堆芯内部不同组件区域的温度差异非常大。对于堆芯中的某一根组件而言，由于其与毗邻组件的间隙非常小，因此该组件的温度场分布会受到毗邻组件温度场分布的影响，导致在该组件的不同高度以及组件的不同周向方向上温度分布也会有较大的区别。在不均匀温度场的作用下，组件局部的热膨胀特性会出现较大的差异，从而使得组件发生热弯曲变形。组件发生热弯曲变形后就会通过垫块与周围的毗邻组件发生接触。如果弯曲变形量较大的话随之产生的接触力也会变大，这样会给钠冷快堆的安全运行带来不利影响，如增加组件插拔力、组件倒换困难、加速组件破损、堆芯象限功率倾斜等。

在钠冷快堆中，中子源组件位于堆芯的中部，因此靠近堆芯的环形区域会有较大的中子通量，因而会释放更多的裂变能，从而该区域的组件会具有较高的温度。越靠近堆芯的外围，中子通量降低，裂变反应释放的裂变能降低，因而外围区域的组件温度也会相应降低。并且这种从堆芯中部到堆外的温度分布趋势沿着堆芯的周向方向上存在着一定的对称性。因此可以独立取出从堆芯中部到堆芯外围的一排组件进行研究，并且会发现这排组件的温度呈现出由堆芯向堆外逐渐降低的趋势。由于温度的不均匀分布，这排组件的会朝着温度降低的方向发生热弯曲变形，并且由于组件之间的距离非常近，组件与组件之间会发生接触，并产生很大的接触力。此外组件之间的相互接触也会阻碍组件原本的热弯曲变形，使得组件最终的变形状态更加难以预测。

由于单排组件的受限热变形过程是一个热力强耦合的过程，因此采用传统的有限元分析软件来分析该过程计算量巨大，且极易发散，而单根组件自由热变形分析程序又无法模拟及预测单排组件热变形的动态耦合过程。基于上述原因，有必要研发钠冷快堆单排组件的耦合热变形模拟技术来预测组件在快堆堆芯中的热变形行为，从而为快堆组件的安全分析提供充分的理论依据。

发明内容

为了克服上述现有技术存在的问题，本发明的目的在于提供针对钠冷快堆单排组件的耦合热变形的模拟方法，通过求解组件的变形协调方程组来预测单排组件在热梯度作用下的耦合变形行为。

为了实现上述目的，本发明采取了以下技术方案：

一种针对的钠冷快堆单排组件耦合热变形的模拟方法，该方法包括以下步骤：

步骤1：根据用户需求，输入单排组件的组件个数、目标组件的物性参数以及几何参数，物性参数指组件的杨氏模量和热膨胀系数，几何参数指组件的长度、套管厚度、外套管的对边距、组件的上下垫块到组件固支端的距离、相邻组件之间的轴心距以及相邻组件上下垫块之间的间隙；

步骤2：将单排组件按照由堆芯朝堆外的方向用数字编号进行标记，并以最靠近堆芯的组件位置为零点，根据相邻组件的轴心距确定单排组件中每一根组件在弯曲的一维方向的初始位置坐标；

步骤3：基于单排组件的三维温度场数据获取单排组件中的每一根组件在一维方向上的自由热变形挠度，针对每一根组件热变形位移的计算方法如下:

(1)输入目标组件的三维温度场数据：三维温度场数据以矩阵的形式给出，其中列标为组件轴向不同位置温度采样点的标高，行标为在组件某一固定标高处进行温度采点的具体位置标记；获得三维温度场数据以后采用插值方法获得轴向均匀分布的温度场数据从而满足计算的输入条件；

(2)对于目标组件的截面进行离散化处理，将其分为12个关于z轴对称的单元，并通过下式获得12个离散的截面单元的静矩即一阶矩；

其中：

A₁到A₁₂为每一个离散单元的横截面积；

z为离散单元的积分域到z轴的距离；

由于组件的截面为中心对称图形，因此在同一横截面上各离散单元之间存在如上所示的关系，因此只需求解相邻的三个离散单元的静矩值即可通过几何关系推得所有12个离散单元的静矩值；

其中相邻的三个离散单元静矩的计算方法如下所示：

其中：

y为离散单元的积分域到y轴的距离；

a₁为套管外对边距，a₂为套管内对边距；

y₁和y₂分别为离散单元7的左、右积分边界；

y₃和y₄分别为离散单元8的右侧子单元的左、右积分边界；

y₅为离散单元9的右积分边界；

(3)获得所有离散单元的静矩以后通过下式求解离散单元截面的总热弯矩；

其中：

M_T为离散单元所在截面的总热弯矩；

x_i为离散单元的轴向坐标；

E为杨氏模量；

z为离散单元的积分域到z轴的距离；

ε_T为轴向位置为x_i的横截面的热应变；

a_i为轴向位置为x_i的横截面的热膨胀系数；

i为轴向坐标的索引号；

j为同一截面上离散单元的索引号；

a_ij为每一个离散单元的热膨胀系数；

T_f为参考温度；

T_i为轴向位置为x_i的横截面的热态温度；

T_ij为每一个离散单元的热态温度；

A为套管的横截面积；

为离散单元的静矩即一阶矩；

然后遍历所有的轴向节点，求得所有轴向节点所在截面的总热弯矩；

(4)通过下式获得每一个轴向节点的横向热变形位移即挠度；

其中：

M_T为离散单元截面的总热弯矩；

w为横向热变形位移即挠度；

x_i为离散单元的轴向坐标；

EI为抗弯刚度；

h为步长即轴向节点间距；

步骤4：根据上、下垫块在组件轴向方向的高度位置，通过对每一根组件进行自由热变形挠度的计算来获取每一根组件在上、下垫块处横向方向上热变形位移的值；

步骤5：根据组件受集中载荷作用下组件横向位移与接触力的关系推导出四种接触力与垫块之间的位移模式，具体方法如下：

由于组件可以视为在集中载荷作用下的悬臂梁，因此组件横向位移与接触力的关系如式所示：

其中a为施加集中外载荷的位置，l为组件试验件的总长度，EI为组件的抗弯刚度，F为集中载荷，y(x)为集中载荷作用下，组件在特定轴向位置的横向位移。

将四种位移模式按照下表进行分类，并给出每种位移模式下的接触力与位移的关系：

表1位移模式分类表

其中h_up为组件固支端到上垫块中部的距离，h_down为组件固支端到下垫块中部的距离；

步骤6：根据每一根组件的初始位置、热变形位移以及接触力位移列出每一根组件在上、下垫块处关于垫块横向位移的平衡方程，具体方法如下：

由于上、下垫块处的最终平衡位置是三部分的共同作用结果：初始位置、热变形位移以及接触力位移；因此针对N根单排受限组件可以列出2N个平衡方程，其中关于上垫块有N个，下垫块同样存在N个；

上垫块的N个平衡方程如下所示：

下垫块的N平衡方程如下所示：

以上两方程中，N为组件编号，为第N个上垫块的最终平衡位置，y_Nu为第N个上垫块的初始位置，为组件热变形导致的第N个上垫块处的位移；同理，为第N个下垫块最终平衡位置，y_Nd为第N个下垫块的初始位置，为组件热变形导致的第N个下垫块处的位移；

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移；

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移；

步骤7：根据相邻组件之间的变形协调特性进一步列出两两相邻的组件关于上、下垫块处的原始变形协调方程，具体方法如下：

组件在发生弯曲变形的时候，垫块之间会相互接触，因此每两两相邻的两根组件的在垫块高度的位置均要满足变形协调条件；

由每两两相邻组件的平衡方程作差相减即可得到组件的变形协调方程，共计2N-2个，其中上垫块N-1个，下垫块N-1个；

上垫块的N-1变形协调方程如下所示：

下垫块的N-1变形协调方程如下所示：

以上两式中，g_u为相邻组件上垫块之间的间隙，g_d为相邻组件下垫块之间的间隙；

步骤8：将组件的变形协调方程中不同接触力对相同垫块的位移项进行合并，即为变形协调方程中括号中的项，然后将括号中的项构造为新的变量，从而获得修正后的变形协调方程便于进行求解，具体方法如下所示：

合并后上垫块处的变形协调方程可以修正为：

合并后下垫块处的变形协调方程修正为：

其中：为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移之和，为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移之和；

而以及由不同接触力产生的位移来表示，如下式：

步骤9：将修正后的变形协调方程写成矩阵的形式，并采用追赶法进行求解；通过求解求出步骤8中所有的yz_u以及yz_d的值；上垫块处的变形协调方程写成矩阵的形式如下：

下垫块处的变形协调方程写成矩阵的形式如下：

步骤10：根据步骤5中四种接触力与垫块之间的位移模式列出2X2的线性方程组来求解每相邻两根组件之间上、下垫块的接触力，进而获得单排所有组件的接触力，具体方法如下：

对于N根组件中的任意两根组件i与i+1，接触力表达式可以通过下式求得：

其中F_u为上垫块接触力，F_d为下垫块接触力，c矩阵为接触力与位移的关系矩阵，其中的每个元素可以由步骤5中的四种位移模式获得；

通过上述方法可以计算每两两相邻组件之间的接触力，进而算出所有的接触力；

步骤11：通过组件的位移叠加原则求得每一根组件所有轴向位置的最终耦合热变形位移，具体计算公式如下：

每根组件的最终耦合热变形位移的具体求解公式如下所示：

上式中：是第i根组件的最终耦合热变形位移，y_i是第i根组件的初始位置，是第i根组件的热变形位移，是第i-1根组件和第i根组件在上垫块处的接触力，是第i根组件和第i+1根组件在上垫块处的接触力，是第i-1根组件和第i根组件在下垫块处的接触力，是第i根组件和第i+1根组件在下垫块处的接触力。

与现有技术相比，本发明有如下突出特点：

1.目前传统的热力耦合分析过程都采用有限元商用软件，例如ABAQUS进行分析，针对钠冷快堆单排组件的分析对象不可避免地会出现计算量大、收敛性差等固有的劣势，而本发明方法采用求解大型稀疏线性方程组的方法来求解接触力，进而获取热变形位移，相比于传统方法具有计算量小，收敛性好等独到的优点；

2.采用传统的ABAQUS有限元软件进行分析时候需要进行包括几何建模、网格划分、参数设置等一系列繁琐的前处理预设过程，而本发明方法则免去了商业软件存在的前处理设置等一系列的麻烦；

3.本发明模拟方法基于Euler-Bernoulli梁理论进行构建，因此相比于传统的有限元方法具有计算效率更高，计算结果更加准确的特点；

4.本发明模拟方法不仅仅适用钠冷快堆组件，针对于不同的计算对象以及物理问题，本模拟方法可以在该基础之上进行扩展。

附图说明

图1为单排组件排列示意图。

图2为组件横截面离散单元示意图。

图3为相邻组件的变形协调示意图。具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明：

本发明的目的在于提供针对钠冷快堆单排组件耦合热变形的模拟方法，通过组件之间相对位置关系构建平衡方程以及相邻组件之间的变形协调方程，进而求解大型稀疏变形协调方程组实现接触力的求解，最终通过接触力数据获得每根组件的实际位移情况，来预测单排组件在热梯度作用下的耦合变形行为。单排组件耦合热变形过程中组件的排列如图1所示。该方法具体流程包括以下方面：

步骤1：根据用户需求，输入单排组件的组件个数、目标组件的物性参数以及几何参数，物性参数指组件的杨氏模量和热膨胀系数，几何参数指组件的长度、套管厚度、外套管的对边距、组件的上下垫块到组件固支端的距离、相邻组件之间的轴心距以及相邻组件上下垫块之间的间隙；

步骤2：将单排组件按照由堆芯朝堆外的方向用数字编号进行标记，并以最靠近堆芯的组件位置为零点，根据相邻组件的轴心距确定单排组件中每一根组件在弯曲的一维方向的初始位置坐标；

步骤3：基于单排组件的三维温度场数据获取单排组件中的每一根组件在一维方向上的自由热变形挠度，针对每一根组件热变形位移的计算方法如下:

(5)输入目标组件的三维温度场数据：三维温度场数据以矩阵的形式给出，其中列标为组件轴向不同位置温度采样点的标高，

行标为在组件某一固定标高处进行温度采点的具体位置标记；获得三维温度场数据以后采用插值方法获得轴向均匀分布的温度场数据从而满足计算的输入条件；

(6)对于目标组件的截面进行离散化处理，将其分为12个关于z轴对称的单元，并通过下式获得12个离散的截面单元的静矩即一阶矩；

其中：

A₁到A₁₂为每一个离散单元的横截面积；

z为离散单元的积分域到z轴的距离；

参考图2中的组件横截面离散单元示意图可知，由于组件的截面为中心对称图形，因此在同一横截面上各离散单元之间存在如上所示的关系，因此只需求解相邻的三个离散单元的静矩值即可通过几何关系推得所有12个离散单元的静矩值；

其中相邻的三个离散单元静矩的计算方法如下所示：

其中：

y为离散单元的积分域到y轴的距离；

a₁为套管外对边距，a₂为套管内对边距；

y₁和y₂分别为离散单元7的左、右积分边界；

y₃和y₄分别为离散单元8的右侧子单元的左、右积分边界；

y₅为离散单元9的右积分边界；

对于组件横截面的离散

(7)获得所有离散单元的静矩以后通过下式求解离散单元截面的总热弯矩；

其中：

M_T为离散单元所在截面的总热弯矩；

x_i为离散单元的轴向坐标；

E为杨氏模量；

z为离散单元的积分域到z轴的距离；

ε_T为轴向位置为x_i的横截面的热应变；

a_i为轴向位置为x_i的横截面的热膨胀系数；

i为轴向坐标的索引号；

j为同一截面上离散单元的索引号

a_ij为每一个离散单元的热膨胀系数；

T_f为参考温度；

T_i为轴向位置为x_i的横截面的热态温度；

T_ij为每一个离散单元的热态温度；

A为套管的横截面积；

为离散单元的静矩即一阶矩；

然后遍历所有的轴向节点，求得所有轴向节点所在截面的总热弯矩；

(8)通过下式获得每一个轴向节点的横向热变形位移即挠度；

其中：

M_T为离散单元截面的总热弯矩；

w为横向热变形位移即挠度；

x_i为离散单元的轴向坐标；

EI为抗弯刚度；

h为步长即轴向节点间距；

步骤4：根据上、下垫块在组件轴向方向的高度位置，通过对每一根组件进行自由热变形挠度的计算来获取每一根组件在上、下垫块处横向方向上热变形位移的值；

步骤5：根据组件受集中载荷作用下组件横向位移与接触力的关系推导出四种接触力与垫块之间的位移模式，具体方法如下：

由于组件可以视为在集中载荷作用下的悬臂梁，因此组件横向位移与接触力的关系如式所示：

其中a为施加集中外载荷的位置，l为组件试验件的总长度，EI为组件的抗弯刚度，F为集中载荷，y(x)为集中载荷作用下，组件在特定轴向位置的横向位移。

将四种位移模式可以按照下表进行分类，并可给出每种位移模式下的接触力与位移的关系：

表1位移模式分类表

其中h_up为组件固支端到上垫块中部的距离，h_down为组件固支端到下垫块中部的距离。

步骤6：根据每一根组件的初始位置、热变形位移以及接触力位移列出每一根组件在上、下垫块处关于垫块横向位移的平衡方程，具体方法如下：

由于上、下垫块处的最终平衡位置是三部分的共同作用结果：初始位置、热变形位移以及接触力位移。因此针对N根单排受限组件可以列出2N个平衡方程，其中关于上垫块有N个，下垫块同样存在N个。

上垫块的N个平衡方程如下所示：

下垫块的N平衡方程如下所示：

以上两方程中，N为组件编号，为第N个上垫块的最终平衡位置，y_Nu为第N个上垫块的初始位置，为组件热变形导致的第N个上垫块处的位移；同理，为第N个下垫块最终平衡位置，y_Nd为第N个下垫块的初始位置，为组件热变形导致的第N个下垫块处的位移。

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移；

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移；

步骤7：根据相邻组件之间的变形协调特性进一步列出两两相邻的组件关于上、下垫块处的原始变形协调方程，具体方法如下：

由图3中相邻组件的截面示意图可知，组件在发生弯曲变形的时候，垫块之间会相互接触，因此每两两相邻的两根组件的在垫块高度的位置均要满足变形协调条件。

由每两两相邻组件的平衡方程作差相减即可得到组件的变形协调方程，共计2N-2个，其中上垫块N-1个，下垫块N-1个。

上垫块的N-1变形协调方程如下所示：

下垫块的N-1变形协调方程如下所示：

以上两式中，g_u为相邻组件上垫块之间的间隙，g_d为相邻组件下垫块之间的间隙。

步骤8：将组件的变形协调方程中不同接触力对相同垫块的位移项进行合并，即为变形协调方程中括号中的项，然后将括号中的项构造为新的变量，从而获得修正后的变形协调方程便于进行求解，具体方法如下所示：

合并后上垫块处的变形协调方程可以修正为：

合并后下垫块处的变形协调方程修正为：

其中：为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移之和，为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移之和；

而以及由不同接触力产生的位移来表示，如下式：

步骤9：将修正后的变形协调方程写成矩阵的形式，并采用追赶法进行求解；通过求解求出步骤8中所有的yz_u以及yz_d的值。上垫块处的变形协调方程写成矩阵的形式如下：

下垫块处的变形协调方程写成矩阵的形式如下：

步骤10：根据步骤5中四种接触力与垫块之间的位移模式列出2X2的线性方程组来求解每相邻两根组件之间上、下垫块的接触力，进而获得单排所有组件的接触力，具体方法如下：

对于N根组件中的任意两根组件i与i+1，接触力表达式可以通过下式求得：

其中F_u为上垫块接触力，F_d为下垫块接触力，c矩阵为接触力与位移的关系矩阵，其中的每个元素可以由步骤5中的四种位移模式获得。

通过上述方法可以计算每两两相邻组件之间的接触力，进而算出所有的接触力。

步骤11：通过组件的位移叠加原则求得每一根组件所有轴向位置的最终耦合热变形位移，具体计算公式如下：

每根组件的最终耦合热变形位移的具体求解公式如下所示：

上式中：是第i根组件的最终耦合热变形位移，y_i是第i根组件的初始位置，是第i根组件的热变形位移，是第i-1根组件和第i根组件在上垫块处的接触力，是第i根组件和第i+1根组件在上垫块处的接触力，是第i-1根组件和第i根组件在下垫块处的接触力，是第i根组件和第i+1根组件在下垫块处的接触力。

Claims (1)

1.一种针对钠冷快堆单排组件耦合热变形的模拟方法，其特征在于：包括以下步骤：

步骤1：根据用户需求，输入单排组件的组件个数、目标组件的物性参数以及几何参数，物性参数指组件的杨氏模量和热膨胀系数，几何参数指组件的长度、套管厚度、外套管的对边距、组件的上下垫块到组件固支端的距离、相邻组件之间的轴心距以及相邻组件上下垫块之间的间隙；

步骤2：将单排组件按照由堆芯朝堆外的方向用数字编号进行标记，并以最靠近堆芯的组件位置为零点，根据相邻组件的轴心距确定单排组件中每一根组件在弯曲的一维方向的初始位置坐标；

步骤3：基于单排组件的三维温度场数据获取单排组件中的每一根组件在一维方向上的自由热变形挠度，针对每一根组件热变形位移的计算方法如下:

(1)输入目标组件的三维温度场数据：三维温度场数据以矩阵的形式给出，其中列标为组件轴向不同位置温度采样点的标高，行标为在组件某一固定标高处进行温度采点的具体位置标记；获得三维温度场数据以后采用插值方法获得轴向均匀分布的温度场数据从而满足计算的输入条件；

(2)对于目标组件的截面进行离散化处理，将其分为12个关于z轴对称的单元，并通过下式获得12个离散的截面单元的静矩即一阶矩；

其中：

A₁到A₁₂为每一个离散单元的横截面积；

z为离散单元的积分域到z轴的距离；

由于组件的截面为中心对称图形，因此在同一横截面上各离散单元之间存在如上所示的关系，因此只需求解相邻的三个离散单元的静矩值即可通过几何关系推得所有12个离散单元的静矩值；

其中相邻的三个离散单元静矩的计算方法如下所示：

其中：

y为离散单元的积分域到y轴的距离；

a₁为套管外对边距，a₂为套管内对边距；

y₁和y₂分别为离散单元7的左、右积分边界；

y₃和y₄分别为离散单元8的右侧子单元的左、右积分边界；

y₅为离散单元9的右积分边界；

(3)获得所有离散单元的静矩以后通过下式求解离散单元截面的总热弯矩；

其中：

M_T为离散单元所在截面的总热弯矩；

x_i为离散单元的轴向坐标；

E为杨氏模量；

z为离散单元的积分域到z轴的距离；

e_T为轴向位置为x_i的横截面的热应变；

a_i为轴向位置为x_i的横截面的热膨胀系数；

i为轴向坐标的索引号；

j为同一截面上离散单元的索引号；

a_ij为每一个离散单元的热膨胀系数；

T_f为参考温度；

T_i为轴向位置为x_i的横截面的热态温度；

T_ij为每一个离散单元的热态温度；

A为套管的横截面积；

为离散单元的静矩即一阶矩；

然后遍历所有的轴向节点，求得所有轴向节点所在截面的总热弯矩；

(4)通过下式获得每一个轴向节点的横向热变形位移即挠度；

其中：

M_T为离散单元截面的总热弯矩；

w为横向热变形位移即挠度；

x_i为离散单元的轴向坐标；

EI为抗弯刚度；

h为步长即轴向节点间距；

步骤4：根据上、下垫块在组件轴向方向的高度位置，通过对每一根组件进行自由热变形挠度的计算来获取每一根组件在上、下垫块处横向方向上热变形位移的值；

步骤5：根据组件受集中载荷作用下组件横向位移与接触力的关系推导出四种接触力与垫块之间的位移模式，具体方法如下：

由于组件被视为在集中载荷作用下的悬臂梁，因此组件横向位移与接触力的关系如式所示：

其中a为施加集中外载荷的位置，l为组件试验件的总长度，EI为组件的抗弯刚度，F为集中载荷，y(x)为集中载荷作用下，组件在特定轴向位置的横向位移；

四种位移模式分别表示为：(1)上垫块存在接触力，上垫块处位移，其分类简称为uu，表达式为(2)下垫块存在接触力，上垫块处位移，其分类简称为du，表达式为(3)上垫块存在接触力，下垫块处位移，其分类简称为ud，表达式为(4)下垫块存在接触力，下垫块处位移，其分类简称为dd，表达式为

其中h_up为组件固支端到上垫块中部的距离，h_down为组件固支端到下垫块中部的距离；

步骤6：根据每一根组件的初始位置、热变形位移以及接触力位移列出每一根组件在上、下垫块处关于垫块横向位移的平衡方程，具体方法如下：

由于上、下垫块处的最终平衡位置是三部分的共同作用结果：初始位置、热变形位移以及接触力位移；因此针对N根单排受限组件能够列出2N个平衡方程，其中关于上垫块有N个，下垫块同样存在N个；

上垫块的N个平衡方程如下所示：

下垫块的N平衡方程如下所示：

以上两方程中，N为组件编号，为第N个上垫块的最终平衡位置，y_Nu为第N个上垫块的初始位置，为组件热变形导致的第N个上垫块处的位移；同理，为第N个下垫块最终平衡位置，y_Nd为第N个下垫块的初始位置，为组件热变形导致的第N个下垫块处的位移；

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的上垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移；

为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i根组件和第i+1根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的上垫块与上垫块之间的接触力造成的第i根组件的下垫块的位移，为第i-1根组件和第i根组件的下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移；

步骤7：根据相邻组件之间的变形协调特性进一步列出两两相邻的组件关于上、下垫块处的原始变形协调方程，具体方法如下：

组件在发生弯曲变形的时候，垫块之间会相互接触，因此每两两相邻的两根组件的在垫块高度的位置均要满足变形协调条件；

由每两两相邻组件的平衡方程作差相减即得到组件的变形协调方程，共计2N-2个，其中上垫块N-1个，下垫块N-1个；

上垫块的N-1变形协调方程如下所示：

下垫块的N-1变形协调方程如下所示：

以上两式中，g_u为相邻组件上垫块之间的间隙，g_d为相邻组件下垫块之间的间隙；

步骤8：将组件的变形协调方程中不同接触力对相同垫块的位移项进行合并，即为变形协调方程中括号中的项，然后将括号中的项构造为新的变量，从而获得修正后的变形协调方程便于进行求解，具体方法如下所示：

合并后上垫块处的变形协调方程修正为：

合并后下垫块处的变形协调方程修正为：

其中：为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的上垫块的位移之和，为第i根组件和第i+1根组件的上垫块与上垫块之间的接触力以及下垫块与下垫块之间的接触力造成的第i根组件的下垫块的位移之和；

而以及由不同接触力产生的位移来表示，如下式：

步骤9：将修正后的变形协调方程写成矩阵的形式，并采用追赶法进行求解；通过求解求出步骤8中所有的yz_u以及yz_d的值；上垫块处的变形协调方程写成矩阵的形式如下：

下垫块处的变形协调方程写成矩阵的形式如下：

步骤10：根据步骤5中四种接触力与垫块之间的位移模式列出2×2的线性方程组来求解每相邻两根组件之间上、下垫块的接触力，进而获得单排所有组件的接触力，具体方法如下：

对于N根组件中的任意两根组件i与i+1，接触力表达式通过下式求得：

其中为组件i的上垫块接触力，为组件i的下垫块接触力，c矩阵为接触力与位移的关系矩阵，其中的每个元素可以由步骤5中的四种位移模式获得；

通过上述方法计算每两两相邻组件之间的接触力，进而算出所有的接触力；

步骤11：通过组件的位移叠加原则求得每一根组件所有轴向位置的最终耦合热变形位移，具体计算公式如下：

每根组件的最终耦合热变形位移的具体求解公式如下所示：

上式中：是第i根组件的最终耦合热变形位移，y_i是第i根组件的初始位置，是第i根组件的热变形位移，是第i-1根组件和第i根组件在上垫块处的接触力，是第i根组件和第i+1根组件在上垫块处的接触力，是第i-1根组件和第i根组件在下垫块处的接触力，是第i根组件和第i+1根组件在下垫块处的接触力。