Page 128 - 2025年第56卷第2期

- 2
- 1
式中：(Δ u n - Δ u n ) 为外荷载增量引起的增量位移差，与当前增量步的接触状态无关，可直接回代得
1
2
1
到；(Δu n - Δu n ) 为当前荷载步增量的总体位移差，与当前增量步的接触状态有关。式（6）中 (Δu n -
Δu n )和 Δf 均未知，需要结合接触条件进行迭代求解。
2
n
2.2 瞬态温度场和应力场的求解方程 瞬态温度场的有限元支配方程可表述为：
R T + HT = P （7）
̇
式中：R 为整体热容矩阵；H 为整体热传导矩阵；P 为整体温度荷载矩阵；T 为各个时刻整体节点温度
值。考虑通水冷却时，温度场的计算采用朱伯芳院士提出的等效算法 ［20］ 。
依据文献［20］知，在只考虑温度荷载的情况下，温度产生的位移有限元方程为：
Ku = R e （8）
e
式中：K 为单元刚度矩阵；u 为由变温引起的节点位移；R 为由变温引起的等效节点荷载。考虑到弹
e
性体内总应变 ε 包括弹性应变 ε 和变温初始应变 ε ，计算温度应力场的方程为：
0
（9）
σ = DBu - Dε 0
2.3 考虑有限元混合法的温度应力耦合场计算方程 联合式（2）（9）得到考虑温度荷载的第 n+1 个荷载
增量步的有限元平衡方程为：
)
Δu n = K -1 ( F n + 1 + f n - ∫ B σ n dΩ + K Δf n + K R n e （10）
T
-1
-1
对于耦合场应力的求解，是根据叠加原理将混合算法求解的应力与考虑温度荷载求解的应力进行
叠加。
（11）
σ T = DBu - Dε 0
2.4 温控接触计算问题的高效算法 接触问题需要结合接触定解条件和约束条件来迭代求解，详细过
程见文献［13］。其中每次迭代需要求解接触面上可能接触对的柔度矩阵元素 c ，其含义是由自由度 j 上
ij
的单位力 F 在自由度 i 上产生的位移 u ，可由单位力法求得：
j
i
Ku = F （12）
系统刚度 K 为大型矩阵，无法直接求逆。利用克洛特法对 K 进行分解：
Ku = LUu = F （13）
其中 L、U 分别为下三角阵和单位上三角矩阵。令
Uu = v （14）
则
Lv = F （15）
给定荷载向量 F，先利用式（15）求解 v，然后求解 u，从而得到每个接触对的柔度矩阵。
对于考虑温度及接触耦合效应的结构静力场问题，假定接触对个数为 n ，则每次迭代需要进行 3×
p
n 次的 LU 回代，再考虑施工期温度场仿真时极为可观的时间增量步数和迭代次数，则总体计算量过
p
于庞大。由于接触点对仅占全部节点的一小部分，因此采用如下的矩阵方程简化计算方法。
（1）对每个施工面貌下的全部激活节点进行重新排序，将接触点置于最后，再进行整体刚度集成
和 LU 分解；
（2）由于 F 的上半部分为零元素，故可记为：
] （16）
F' = [ 0 … F' c
式（15）可写成
O
é ù
ê ê é□ O ù é ù ú = ê ú ú ú （17）
ê ê ú
ú ú ê ê ê ú ú ú
ê
ê
ê ê
ú ú
v 1
ë□ L c û ë û ë û
F c
v c
式 中 ： O 指 代 全 零 子 矩 阵 ； □ 指 代 其 他 子 矩 阵 ； F' c 指 代 F 后 半 非 零 部 分 。 由 矩 阵 运 算 易 知 ， v =
1
O，故
O
O
ê ê é□ O ù é ù ú ú é ù ú ú （18）
ê
ú ú ê ú
ê ê ú
ú ú ê ú = ê ú
ê
ê
ê ê
ë□ L c û ë û ë û
F c
v c
— 268 —