科学技术与工程

类型	主要形式	处理方式
Dirichlet 边界^[17]	固定约束	直接消去法
指定位移	对角元素改1,其他元素置0
地震动作用	地震作用的施加 (一致激励、非一致激励、人工边界地震动输入)	本质上是一种运动着的边界约束条件,直接代入计算
主从约束	单级约束, (刚性约束、柔性约束)	将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]
多级约束^[18],
时变系数约束
非线性约束	大位移大变形问题,	转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]
曲线或曲面滑动支承,
接触碰撞问题

类型

主要形式

处理方式

Dirichlet
边界^[17]

固定约束

直接消去法

指定位移

对角元素改1,其他元素置0

地震动
作用

地震作用的施加
(一致激励、非一致激励、
人工边界地震动输入)

本质上是一种运动着的边界约束条件,直接代入计算

主从
约束

单级约束,
(刚性约束、柔性约束)

将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]

多级约束^[18],

时变系数约束

非线性
约束

大位移大变形问题,

转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]

曲线或曲面滑动支承,

接触碰撞问题

类型	主要形式	处理方式
Dirichlet 边界^[17]	固定约束	直接消去法
指定位移	对角元素改1,其他元素置0
地震动作用	地震作用的施加 (一致激励、非一致激励、人工边界地震动输入)	本质上是一种运动着的边界约束条件,直接代入计算
主从约束	单级约束, (刚性约束、柔性约束)	将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]
多级约束^[18],
时变系数约束
非线性约束	大位移大变形问题,	转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]
曲线或曲面滑动支承,
接触碰撞问题

类型

主要形式

处理方式

Dirichlet
边界^[17]

固定约束

直接消去法

指定位移

对角元素改1,其他元素置0

地震动
作用

地震作用的施加
(一致激励、非一致激励、
人工边界地震动输入)

本质上是一种运动着的边界约束条件,直接代入计算

主从
约束

单级约束,
(刚性约束、柔性约束)

将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]

多级约束^[18],

时变系数约束

非线性
约束

大位移大变形问题,

转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]

曲线或曲面滑动支承,

接触碰撞问题

模型	N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
1	1	1	24	23	1	—	—	—
2	2	1	24	23	1	12	12	1
		2	24	22	2

模型

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

—

模型	N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
1	1	1	24	23	1	—	—	—
2	2	1	24	23	1	12	12	1
		2	24	22	2

模型

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

—

N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
3	1	3	3	0	4 4	3 3	1 1
2	6	5	1
3	8	6	2

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

4
4

3
3

1
1

N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
3	1	3	3	0	4 4	3 3	1 1
2	6	5	1
3	8	6	2

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

4
4

3
3

1
1

第一阶段:
1.初始化 1-1.初始化K_II块、v_alue_K_II[ $K I I n u m$ ]、r_ow_K_II[ $K I I n u m$ ]、 c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber]; 1-2.临时变量N_um初始化为0。(记录元素位置)。
2.遍历E₀(1,1,N_E) 2-1.获取当前单元E_LEM(E₀)_i; 2-2.获取当前单元独立自由度个数D_i及相应K_II子块; 2-3.判断单元是否循环完毕,若循环完毕,集成结束;若循环没有完毕,进入操作3。
3.遍历K_II子块 3-1.获取子块矩阵元素v_alue、r_ow、c_ol; 3-2.判断r_ow是否小于c_ol,若小于进行3-3,若大于返回3-1; 3-3.计算o_ff=U_BW(i-1)+j-i+1,并判断o_ff是否存在于B_and_K_II[o_ff,N_umber]中,若存在,根据o_ff索引出N_umber,将v_alue_K_II[N_umber]进行数值叠加;若不存在,存储相应数据:v_alue_K_II[N_um]、r_ow_K_II[N_um]=r_ow、c_ol.height_K_II[c_ol]++、 B_and_K_II[o_ff,N_umber](默认第一次遍历不存在,直接存储即可)。
4.单元循环完毕,返回操作操作2

第一阶段:

1.初始化
1-1.初始化K_II块、v_alue_K_II[

K I I n u m

]、r_ow_K_II[

K I I n u m

]、
c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber];
1-2.临时变量N_um初始化为0。(记录元素位置)。

2.遍历E₀(1,1,N_E)
2-1.获取当前单元E_LEM(E₀)_i;
2-2.获取当前单元独立自由度个数D_i及相应K_II子块;
2-3.判断单元是否循环完毕,若循环完毕,集成结束;若循环没有完毕,进入操作3。

3.遍历K_II子块
3-1.获取子块矩阵元素v_alue、r_ow、c_ol;
3-2.判断r_ow是否小于c_ol,若小于进行3-3,若大于返回3-1;
3-3.计算o_ff=U_BW(i-1)+j-i+1,并判断o_ff是否存在于B_and_K_II[o_ff,N_umber]中,若存在,根据o_ff索引出N_umber,将v_alue_K_II[N_umber]进行数值叠加;若不存在,存储相应数据:v_alue_K_II[N_um]、r_ow_K_II[N_um]=r_ow、c_ol.height_K_II[c_ol]++、
B_and_K_II[o_ff,N_umber](默认第一次遍历不存在,直接存储即可)。

4.单元循环完毕,返回操作操作2

第一阶段:
1.初始化 1-1.初始化K_II块、v_alue_K_II[ $K I I n u m$ ]、r_ow_K_II[ $K I I n u m$ ]、 c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber]; 1-2.临时变量N_um初始化为0。(记录元素位置)。
2.遍历E₀(1,1,N_E) 2-1.获取当前单元E_LEM(E₀)_i; 2-2.获取当前单元独立自由度个数D_i及相应K_II子块; 2-3.判断单元是否循环完毕,若循环完毕,集成结束;若循环没有完毕,进入操作3。
3.遍历K_II子块 3-1.获取子块矩阵元素v_alue、r_ow、c_ol; 3-2.判断r_ow是否小于c_ol,若小于进行3-3,若大于返回3-1; 3-3.计算o_ff=U_BW(i-1)+j-i+1,并判断o_ff是否存在于B_and_K_II[o_ff,N_umber]中,若存在,根据o_ff索引出N_umber,将v_alue_K_II[N_umber]进行数值叠加;若不存在,存储相应数据:v_alue_K_II[N_um]、r_ow_K_II[N_um]=r_ow、c_ol.height_K_II[c_ol]++、 B_and_K_II[o_ff,N_umber](默认第一次遍历不存在,直接存储即可)。
4.单元循环完毕,返回操作操作2

第一阶段:

1.初始化
1-1.初始化K_II块、v_alue_K_II[

K I I n u m

]、r_ow_K_II[

K I I n u m

]、
c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber];
1-2.临时变量N_um初始化为0。(记录元素位置)。

4.单元循环完毕,返回操作操作2

第二阶段:
1.初始化 1-1.临时变量N初始化为0。(记录循环次数);初始化 $R o w_K I I$ [ $K I I n u m$ ]、 $V a l u e_K I I$ [ $K I I n u m$ ]、 $C o l . o f f_K I I$ [D_i+1]。
2.遍历 $B a n d_K I I$ [o_ff,N_umber] 2-1.获取当前o_ff对应的编号N_umber; 2-2.获取 $r o w_K I I$ [N_umber]、 $v a l u e_K I I$ [N_umber]; 2-3. $R o w_K I I$ [N]= $r o w_K I I$ [N_umber]、 $V a l u e_K I I$ [N]= $v a l u e_K I I$ [N_umber];判断N是否小于D_i+1,若小于,将 $c o l . h e i g h t_K I I$ [D_i]的前N个数值进行叠加并赋值给数组 $C o l . o f f_K I I$ [N+1]。 2-4.变量N++,返回操作2-1。

第二阶段:

1.初始化
1-1.临时变量N初始化为0。(记录循环次数);初始化

R o w_K I I

[

K I I n u m

]、

V a l u e_K I I

[

K I I n u m

]、

C o l . o f f_K I I

[D_i+1]。

2.遍历

B a n d_K I I

[o_ff,N_umber]
2-1.获取当前o_ff对应的编号N_umber;
2-2.获取

r o w_K I I

[N_umber]、

v a l u e_K I I

[N_umber];
2-3.

R o w_K I I

[N]=

r o w_K I I

[N_umber]、

V a l u e_K I I

[N]=

v a l u e_K I I

[N_umber];判断N是否小于D_i+1,若小于,将

c o l . h e i g h t_K I I

[D_i]的前N个数值进行叠加并赋值给数组

C o l . o f f_K I I

[N+1]。
2-4.变量N++,返回操作2-1。

第二阶段:
1.初始化 1-1.临时变量N初始化为0。(记录循环次数);初始化 $R o w_K I I$ [ $K I I n u m$ ]、 $V a l u e_K I I$ [ $K I I n u m$ ]、 $C o l . o f f_K I I$ [D_i+1]。
2.遍历 $B a n d_K I I$ [o_ff,N_umber] 2-1.获取当前o_ff对应的编号N_umber; 2-2.获取 $r o w_K I I$ [N_umber]、 $v a l u e_K I I$ [N_umber]; 2-3. $R o w_K I I$ [N]= $r o w_K I I$ [N_umber]、 $V a l u e_K I I$ [N]= $v a l u e_K I I$ [N_umber];判断N是否小于D_i+1,若小于,将 $c o l . h e i g h t_K I I$ [D_i]的前N个数值进行叠加并赋值给数组 $C o l . o f f_K I I$ [N+1]。 2-4.变量N++,返回操作2-1。

第二阶段:

1.初始化
1-1.临时变量N初始化为0。(记录循环次数);初始化

R o w_K I I

[

K I I n u m

]、

V a l u e_K I I

[

K I I n u m

]、

C o l . o f f_K I I

[D_i+1]。

2.遍历

B a n d_K I I

[o_ff,N_umber]
2-1.获取当前o_ff对应的编号N_umber;
2-2.获取

r o w_K I I

[N_umber]、

v a l u e_K I I

[N_umber];
2-3.

R o w_K I I

[N]=

r o w_K I I

[N_umber]、

V a l u e_K I I

[N]=

v a l u e_K I I

[N_umber];判断N是否小于D_i+1,若小于,将

c o l . h e i g h t_K I I

[D_i]的前N个数值进行叠加并赋值给数组

C o l . o f f_K I I

[N+1]。
2-4.变量N++,返回操作2-1。

位置	OpenSees/ cm	CSC存储/ cm	变带宽存储/cm	误差1/ %	误差2/ %
0#墩顶	1.079 30	1.079 00	1.079 00	0.027 0	0.027 0
1#墩顶	1.079 10	1.079 07	1.079 07	0.002 7	0.002 7
2#墩顶	1.081 44	1.081 41	1.081 41	0.002 7	0.002 7
3#墩顶	1.084 60	1.084 57	1.084 57	0.002 7	0.002 7
4#墩顶	1.102 53	1.102 48	1.102 48	0.004 5	0.004 5

位置

OpenSees/
cm

CSC存储/
cm

变带宽存
储/cm

误差1/
%

误差2/
%

0#墩顶

1.079 30

1.079 00

0.027 0

1#墩顶

1.079 10

1.079 07

0.002 7

2#墩顶

1.081 44

1.081 41

0.002 7

3#墩顶

1.084 60

1.084 57

0.002 7

4#墩顶

1.102 53

1.102 48

0.004 5

位置	OpenSees/ cm	CSC存储/ cm	变带宽存储/cm	误差1/ %	误差2/ %
0#墩顶	1.079 30	1.079 00	1.079 00	0.027 0	0.027 0
1#墩顶	1.079 10	1.079 07	1.079 07	0.002 7	0.002 7
2#墩顶	1.081 44	1.081 41	1.081 41	0.002 7	0.002 7
3#墩顶	1.084 60	1.084 57	1.084 57	0.002 7	0.002 7
4#墩顶	1.102 53	1.102 48	1.102 48	0.004 5	0.004 5

位置

OpenSees/
cm

CSC存储/
cm

变带宽存
储/cm

误差1/
%

误差2/
%

0#墩顶

1.079 30

1.079 00

0.027 0

1#墩顶

1.079 10

1.079 07

0.002 7

2#墩顶

1.081 44

1.081 41

0.002 7

3#墩顶

1.084 60

1.084 57

0.002 7

4#墩顶

1.102 53

1.102 48

0.004 5

考虑约束的CSC格式总体矩阵存储算法

PDF下载

董正方 ¹ , 王佳豪 ¹ , 李运华 ¹^,^* , 程纬 ² , 王君杰 ²

科学技术与工程 | 论文·建筑科学 2025,25(10): 4265-4273

收起

科学技术与工程 | 论文·建筑科学 2025, 25(10): 4265-4273

考虑约束的CSC格式总体矩阵存储算法

全屏

董正方¹, 王佳豪¹, 李运华¹^,^*, 程纬², 王君杰²

作者信息

¹ 河南大学建筑工程学院, 开封 475004

² 同济大学桥梁工程系, 上海 200092

董正方(1980—),男,汉族,河南滑县人,博士,教授。研究方向:有限元软件;桥梁抗震。E-mail:dzf@henu.edu.cn。

通讯作者:

* 李运华(1980—),男,汉族,河南项城人,硕士,高级实验师。研究方向:材料学。E-mail:24561591@qq.com。

CSC Format Global Matrix Storage Algorithm Considering Constraints

Zheng-fang DONG¹, Jia-hao WANG¹, Yun-hua LI¹^,^*, Wei CHENG², Jun-jie WANG²

Affiliations

¹ School of Civil Engineering and Architecture, Henan University, Kaifeng 475004, China

² Department of bridge engineering, Tongji University, Shanghai 200092, China

出版时间: 2025-04-08 doi: 10.12404/j.issn.1671-1815.2403703

文章导航

摘要

收起

为了提高有限元分析效率,减少内存消耗,对总体矩阵的存储算法进行研究。将约束条件统一处理为单元层面的矩阵分块,后续以刚度矩阵为例,确定了约束条件下刚度矩阵非零元素分布规律,并提出了适用于二维、三维有限元刚度矩阵非零元素数量计算方法及相应计算公式,同时检验了正确性;推导了一维等带宽存储与原始方阵地址对应关系,并将其作为辅助数组应用于CSC格式总体矩阵集成过程中,提高了总体矩阵组集效率,避免了节点编号优化难题。最后,以北京市房山线长稻区间四跨单线梁桥验证所提出算法的正确性与实用性,结果表明:该算法分析时总体矩阵组集时间减少30%、存储空间节省68%以上、保证了计算精度的同时线性方程组计算效率提升超过71%。

关键词

有限单元法 / 刚度矩阵 / 非零元素 / 压缩存储 / Cholesky分解

Abstract

收起

In order to improve the efficiency of finite element analysis and reduce memory consumption, the storage algorithm of the overall matrix was studied. The constraints were uniformly processed as matrix partitions at the element level. Subsequently, taking the stiffness matrix as an example, the distribution pattern of non-zero elements under the constrained condition was determined. A calculation method and corresponding formulas suitable for the number of non-zero elements in the stiffness matrix of two-dimensional and three-dimensional finite elements were proposed, and the correctness was verified. The correspondence between one-dimensional equal bandwidth storage and the original square matrix address was deduced, and it was applied as an auxiliary array in the process of integrating the overall matrix in CSC format, which has improved the efficiency of the overall matrix assembly and avoided the difficulty of node numbering optimization. Finally, the correctness and practicability of the proposed algorithm were verified by using thefour-span single-line girder bridge in the Changdao section of the Fangshan Line in Beijing. The results show that the overall matrix assembly time is reduced by 30%, the storage space is saved by more than 68%, and the calculation efficiency of linear equations is improved by more than 71% while ensuring the calculation accuracy.

Key words

finite element method / stiffness matrix / non-zero elements / compressed storage / cholesky decomposition

引用本文

董正方, 王佳豪, 李运华, 程纬, 王君杰. 考虑约束的CSC格式总体矩阵存储算法. 科学技术与工程, 2025 , 25 (10) : 4265 -4273 . DOI: 10.12404/j.issn.1671-1815.2403703

Zheng-fang DONG, Jia-hao WANG, Yun-hua LI, Wei CHENG, Jun-jie WANG. CSC Format Global Matrix Storage Algorithm Considering Constraints[J]. Science Technology and Engineering, 2025 , 25 (10) : 4265 -4273 . DOI: 10.12404/j.issn.1671-1815.2403703

正文

收起

随着工程技术的快速发展和数值计算方法的不断研究,有限元法逐渐成为数值仿真的主流方法。近年来,工程规模不断扩大,所需要分析的工程问题日益复杂,人们必须考虑内存资源占用与计算时间成本的影响,如何进一步提高有限元分析效率已成为人们关注的热点。有限元分析最终归结为求解一个大型线性方程组,其求解效率取决于线性方程组的解法和储存结构^[1-3]。现有的一些存储方法,如等带宽存储、变带宽存储等都受到节点编号的影响,若节点编号设置不合理则会增加有限元分析时间和内存开销;若对节点编号进行优化,则会耗费一定的时间^[4-5]。

目前,节省内存开销和提高分析效率主要有两个方面。一方面是减少存储数量,只存储非零元素。通常根据ID表法、细胞元和广义相邻节点对的概念确定总体刚度矩阵中的非零元素数量和位置^[6-8];这些方法避免了节点编号优化的复杂问题,但需根据单元模型的拓扑结构形成索引数组,增加了计算时间;并且这些方法适用性较差,仅适用于实体单元、特定单元分析以及没有复杂约束的模型。另一方面是改进压缩存储方法和技术,如一维变带宽、三元组、十字链表法^[9]等。然而,随着模型复杂程度的增加,这些方法无法满足工程计算人员对内存资源的要求。因此,列压缩(compressed sparse column,CSC)和行压缩(compressed sparse row,CSR)成为主流^[10]。为了进一步优化内存开销,研究人员又提出了多种新的压缩存储技术,包括改进一维变带宽存储、MSR(modified sparse row)、改进S(B)CSR^[11][streamed (block) compressed sparse row]和负号CSR存储^[12]等;后续随着计算机技术的发展,并行计算、机器学习与存储技术及矩阵运算的结合也成为实现高效率存储与高效率计算的主要手段^[13-16];尽管如此,但上述方法仍有不足,如改进S(B)CSR在稀疏矩阵向量乘法(sparse matrix-vector multiplication,SpMV)和稀疏矩阵处理上表现优异,但局限性强,需要特殊的工具与方法;负号CSR存储目前只适用于三维有限元;而与并行计算、机器学习结合的方法在稀疏矩阵处理上仍有欠缺,主要在SpMV计算上表现优异等。

为提升有限元分析的效率并降低其在分析过程中的内存需求,现从自由度的角度出发,对受约束条件影响的总体矩阵进行分块,消除所有约束自由度,缩减矩阵规模;后续以刚度矩阵为例,根据其组成特点,提出一种适用于二维、三维刚度矩阵非零元素确定方法,并给出通用计算公式;将等带宽移植原则与CSC存储相结合,仅存储总体刚度矩阵中的非零元素,可以高效地直接集成CSC压缩格式的总体刚度矩阵。

1 总体矩阵非零元素确定

收起

1.1 考虑约束的总体矩阵分块

在有限元分析中,结构的约束形式及其处理方式对于在模拟实际物理行为时至关重要。常见的约束形式及其处理方式如表 1所示。

对于复杂的结构形式及不同的边界条件或约束形式,节点位移形态迥异,为了有效地处理这类问题,钟万勰^[25]提出了位移规格数的概念,根据节点位移自由度将其分为5种规格。为了明确辨析地震作用的影响,对几何可动自由度进一步分类为独立自由度与地震动自由度、将受主从约束及非线性约束影响的自由度归类为从自由度,同时将受Dirichlet边界限制的自由度细分为固定自由度与指定位移自由度,共5种规格。

由表1可知,处理固定、指定位移、从自由度等约束时,通常采用直接消去法、罚函数法、拉格朗日乘子法等在总体层面对矩阵进行处理。然而,对于大型矩阵,在总体层面进行处理时会出现工作繁琐、计算量大的情况。针对于此,采用在单元层面统一将约束条件处理转换为单元层面的矩阵分块与总体矩阵分块,在装配总体矩阵时将所有约束自由度消除。假设总体矩阵为对称矩阵,单元的全局位移自由度向量为

(1)δ_e={

u 1 u 2 ··· u n

}

式(1)中:u₁、u₂、···、u_n分别为单元自由度1、2、···、n的位移。

考虑一般情况,式(1)中位移元素包含了独立自由度、约束自由度(包括固定、指定位移、从自由度)和地震动自由度,经过多次行变换,即左乘转换矩阵H^T,且H^T=H^-1,将自由度向量分类排序为

(2)$\boldsymbol{\delta}_{\mathrm{e}}=\boldsymbol{H}\left\{\begin{array}{lll} \boldsymbol{\delta}_{\mathrm{e}}^{\mathrm{I}} & \boldsymbol{\delta}_{\mathrm{e}}^{\mathrm{C}} & \boldsymbol{\delta}_{\mathrm{e}}^{\mathrm{G}} \end{array}\right\}^{\mathrm{T}}$

式(2)中:

δ e I

、

δ e C

、

δ e G

分别单元的独立、约束和地震动自由度位移向量。

对于复杂的约束条件方程,一般的约束位移自由度可以表示为

(3)$\begin{aligned} u_{j}^{\mathrm{C}}= & \left(\alpha_{j 1} u_{1}^{\mathrm{I}}+\alpha_{j 2} u_{2}^{\mathrm{I}}+\cdots+\alpha_{j i} u_{i}^{\mathrm{I}}\right)+ \\ & \left(\beta_{j 1} u_{1}^{\mathrm{G}}+\beta_{j 2} u_{2}^{\mathrm{G}}+\cdots+\beta_{j i} u_{i}^{\mathrm{G}}\right)+ \\ & \left(\gamma_{j 1} u_{1}^{\mathrm{P}}+\gamma_{j 2} u_{2}^{\mathrm{P}}+\cdots+\gamma_{j i} u_{i}^{\mathrm{P}}\right)+ \\ & \left(\eta_{j 1} u_{1}^{\mathrm{F}}+\eta_{j 2} u_{2}^{\mathrm{F}}+\cdots+\eta_{j i} u_{i}^{\mathrm{F}}\right)+v_{j}^{\mathrm{C}} \end{aligned}$

式(3)中:

u I i

、

u G i

、

u P i

、

u F i

分别为与约束位移自由度j相关的独立、地震动、固定约束及指定位移约束自由度位移;α_ji、β_ji、γ_ji、η_ji为相应自由度关系系数,其中i=1,2,···,i;

v C j

为非其次项,这里是常数。固定约束部分为零,所以,约束自由度有以下转换方式。

(4)$\begin{aligned} \boldsymbol{\delta}_{j}^{\mathrm{C}}= & {\left[\begin{array}{c} u_{1}^{\mathrm{C}} \\ u_{2}^{\mathrm{C}} \\ \vdots \\ u_{k}^{\mathrm{C}} \end{array}\right]=\left[\begin{array}{cccc} \alpha_{11} & \alpha_{12} & \cdots & \alpha_{1 i} \\ \alpha_{21} & \alpha_{22} & \cdots & \alpha_{2 i} \\ \vdots & \vdots & & \vdots \\ \alpha_{k 1} & \alpha_{k 2} & \cdots & \alpha_{k i} \end{array}\right]\left[\begin{array}{c} u_{1}^{\mathrm{I}} \\ u_{2}^{\mathrm{I}} \\ \vdots \\ u_{k}^{\mathrm{I}} \end{array}\right]+} \\ & {\left[\begin{array}{cccc} \beta_{11} & \beta_{12} & \cdots & \beta_{1 i} \\ \beta_{21} & \beta_{22} & \cdots & \beta_{2 i} \\ \vdots & \vdots & & \vdots \\ \beta_{k 1} & \beta_{k 2} & \cdots & \beta_{k i} \end{array}\right]\left[\begin{array}{c} u_{1}^{\mathrm{G}} \\ u_{2}^{\mathrm{G}} \\ \vdots \\ u_{k}^{\mathrm{G}} \end{array}\right]+} \\ & {\left[\begin{array}{cccc} \gamma_{11} & \gamma_{12} & \cdots & \gamma_{1 i} \\ \gamma_{21} & \gamma_{22} & \cdots & \gamma_{2 i} \\ \vdots & \vdots & & \vdots \\ \gamma_{k 1} & \gamma_{k 2} & \cdots & \gamma_{k i} \end{array}\right]\left[\begin{array}{c} u_{1}^{\mathrm{P}} \\ u_{2}^{\mathrm{P}} \\ \vdots \\ u_{k}^{\mathrm{P}} \end{array}\right]+\left[\begin{array}{c} v_{1}^{\mathrm{C}} \\ v_{2}^{\mathrm{C}} \\ \vdots \\ v_{k}^{\mathrm{C}} \end{array}\right] } \\ = & \boldsymbol{A} \overline{\boldsymbol{\delta}}_{\mathrm{e}}^{\mathrm{I}}+\boldsymbol{B} \overline{\boldsymbol{\delta}}_{\mathrm{e}}^{\mathrm{G}}+\boldsymbol{C} \overline{\boldsymbol{\delta}}_{\mathrm{e}}^{\mathrm{P}}+\boldsymbol{V}^{\mathrm{C}} \end{aligned}$

式(4)中:A、B、C分别为约束独立自由度关系矩阵、约束地震动自由度关系矩阵和约束指定位自由度关系矩阵。

将式(4)代入式(2),单元位移自由度向量可表示为

(5)δ_e=H

δ e I δ e C δ e G

I I O A B O I G δ - e I δ - e G

O C δ - e P + V C O

式(5):O为零矩阵;I_I、I_G均为eye矩阵(不一定是方阵)。

将式(5)对时间求一阶导数与二阶导数,即

(6)

δ · e = H δ · e I δ · e C δ · e G = H I I O A B O I G δ - · e I δ - · e G δ · · e = H δ · · e I δ · · e C δ · · e G = H I I O A B O I G δ - · · e I δ - · · e G

对于动力问题,考虑瞬时势能泛函,即

(7)Π_e=

12 δ e T

K_eδ_e-

δ e T

(P_e-M_e

δ · · e

-C_e

δ · e

)

式(7)中:K_e、M_e、C_e分别为单元刚度、质量、阻尼矩阵;P_e为单元荷载向量。

将式(5)、式(6)代入式(7),根据动力最小势能原理^[26]可得引入了约束条件后的单元刚度、质量、阻尼矩阵为

(8)

K - e = I I T A T O T O T B T I G T H T K e H I I O A B O I G C - e = I I T A T O T O T B T I G T H T C e H I I O A B O I G M - e = I I T A T O T O T B T I G T H T M e H I I O A B O I G

将矩阵按自由度类型分块,即

(9)

K - e = K I I e K I G e K G I e K G G e C - e = C I I e C I G e C G I e C G G e M - e = M I I e M I G e M G I e M G G e

式(9)中:

K I I e

为单元刚度矩阵中独立自由度模块;

K G G e

为单元刚度矩阵中地震动自由度模块;

K I G e

、

K G I e

为单元刚度矩阵中地震动自由度对独立自由度影响模块。单元阻尼矩阵、单元质量矩阵与单元刚度矩阵分块相同。

K - e

、

C - e

、

K - e

都是方阵,阶数为单元独立自由度数与地震动自由度数总和。将各个单元矩阵块分类组装可以得到总体矩阵为

(10)

K = K I I K I G K G I K G G C = C I I C I G C G I C G G M = M I I M I G M G I M G G

式(10)中:K、M、C分别为总体刚度、质量、阻尼矩阵。

对于静力问题,不考虑惯性力和阻尼力的影响,刚度矩阵表达式与动力问题相同。

综上所述,可以看出无论是静力问题还是动力问题,经过约束处理后,只有独立自由度、地震动自由度对单元矩阵和总体矩阵有贡献,一定程度上减少了总体矩阵的规模及相应存储需求。并且从刚度矩阵、质量矩阵、阻尼矩阵计算方式与分块形式上可以看出,3种矩阵原理相同,后续以刚度矩阵为例。

1.2 非零元素确定

文献[7]采用“广义相邻节点对”确定了三维实体单元的非零子矩阵数量与位置,本文研究对其进行推广。在同一单元内,任意两个节点是相关的,进而任意两个自由度也相关。单元刚度矩阵的每个矩阵数据可由这些“相关自由度”所确定。“相关自由度”数量与总体刚度矩阵上三角部分存储的非零元素数量一致,且一一对应。所以,通过分析自由度的关联性,可以准确地确定刚度矩阵非零元素的数量。一旦确定了非零元素的数据量,即可精确地分配总体矩阵的存储空间,这对于总体矩阵存储算法的实现具有重要意义。

由刚度矩阵特性可知,单元之间共同节点只影响矩阵数值,不会改变矩阵非零元素数量。因此,当模型由单个单元组成时,非零元素数量完全取决于该单元的自由度数量;如果考虑模型由多个单元组成时,单元与单元之间有共节点,非零元素数量将由单元与单元之间增加自由度数量确定。可得具体计算公式为

(11)

K I I n u m

∑ i = 1 N D i (D i + 1) 2 - ∑ i = 2 N [D i - n i (i - 1) D] [D i - n i (i - 1) D + 1] 2, N ≠ 1 D i (D i + 1) 2, N = 1

(12)

K G G n u m

∑ i = 1 N G i (D i + 1) 2 - ∑ i = 2 N [G i - n i (i - 1) G] [G i - n i (i - 1) G + 1] 2, N ≠ 1 G i (G i + 1) 2, N = 1

(13)

K I G n u m

∑ i = 1 N G i D i 2 - ∑ i = 2 N [D i - n i (i - 1) D] [G i - n i (i - 1) G] 2, N ≠ 1 G i D i, N = 1

(14)K_Num=

∑ i = 1 N Y i (Y i + 1) 2 - ∑ i = 2 N [Y i - n i (i - 1)] [Y i - n i (i - 1) + 1] 2, N ≠ 1 Y i (Y i + 1) 2, N = 1

(15)K_Num=

K I I n u m

K G G n u m

K I G n u m

式中:

K x x n u m

为各子块非零元素数量;K_Num为总体刚度矩阵的非零刚度矩阵元素数量;D_i、G_i与Y_i(Y_i=D_i+G_i)为第i个单元的独立自由度、地震动自由度与或总自由度数量;

n i (i - 1) D

、

n i (i - 1) G

与n_i₍_i_-1)[n_i₍_i_-1)=

n i (i - 1) D

n i (i - 1) G

]为第i个单元相较于前个i-1单元增加的独立自由度、地震动自由度或总自由度数量;N为单元数量。

为验证式(11)~式(15)的正确性,建立如图1所示的八节点六面体单元模型。各参数取值如表2所示。

由式(11)计算可得,

K I I n u m 1

=276、

K G G n u m 1

=1、

K I G n u m 1

=3、K_num1=300;

K I I n u m 2

=463、

K G G n u m 1

=3、

K I G n u m 2

=56、K_num2=522。文献[7]中,对于图 1所示模型计算的非零元素数量分别为300和522,上述计算结果与文献[7]完全一致。

进一步,为验证式(11)~式(15)的通用性,建立梁单元与平面单元连接的二维模型,如图2所示。各参数取值如表3所示。

由式(11)计算可得,

K I I n u m

=33、

K G G n u m

=3、

K I G n u m

=14、 K_num=50。可以看出公式的计算结果与图 2非零元素总数相同,验证了公式的正确性与通用性。

2 CSC格式总体刚度矩阵直接集成

收起

2.1 等带宽一维数组定位原理

等带宽存储是有限元中常用的矩阵存储格式。如图3所示,矩阵按照等带宽存储时,矩阵数组中的每一行仍是原来U_BW个元素,只是矩阵的主对角元素都移到了矩阵数组中的第1列。

用i、 j表示方阵形式的行号和列号,k_ij表示矩阵元素。按等带宽存储后,行号记为r_ow,列号记为c_ol。此时,方阵形式与等带宽形式的二维矩阵元素对应关系为

(12)

r o w = i c o l = j - i + 1

将等带宽矩阵数组以一维数组B表示,即

(13)B={

k 11 k 12 ··· 0 k 99 000

}

按等带宽一维存储形式储存后,在一维数组中的偏移记为o_ff。此时,方阵形式与等带宽形式的一维矩阵元素对应关系为

(14)o_ff=U_BW(r_ow-1)+c_ol=U_BW(i-1)+j-i+1

只保存非零元素,此时矩阵元素k_ij、元素位置N_umber与偏移位置o_ff的对应关系如图4所示。

可以看出,方阵中的每一个矩阵元素在等带宽一维存储形式下都有一个独立o_ff与N_umber,可以通过o_ff来找寻方阵中的某一矩阵元素k_ij与元素位置N_umber;相反,也可以通过方阵中的某一矩阵元素k_ij精准定位元素位置N_umber与偏移信息o_ff。

因此,可以通过遍历单元刚度矩阵信息计算偏移信息o_ff,并根据o_ff大小进行排序获取位置信息N_umber,通过o_ff与N_umber的唯一性实现矩阵元素定位、排序以及总体矩阵的集成。为减少计算的工作量且保证偏移信息的唯一性,该算法不再计算带宽U_BW,将U_BW取模型自由度个数。

2.2 CSC存储

该方法采用三个一维数组存储矩阵数据,即(R_ow,C_ol.off,V_alue)。V_alue数组存储矩阵中的非零元素,矩阵按列划分为n个段,每段包含该列的所有非零元素,长度为矩阵非零元素个数;R_ow记录矩阵非零元素相应的行号,按列划分为n段,与V_alue的长度相同,数组的每段都是矩阵行中非零元素的行索引;C_ol.off为矩阵表示某列的起始元素在V_alue中的偏移位置,长度为矩阵的阶数加一。

以一个包含11个非零元素的6×6二维稀疏矩阵A为例,根据上述存储方法,该矩阵的CSC存储如图5所示。

2.3 刚度矩阵的集成

一般情况下,在刚度矩阵生成时,首先需要根据单元的拓扑结构来确定非零子矩阵的数量和位置,生成索引数组,根据刚度矩阵的生成规则和索引数组,将刚度矩阵数据对号入座进行组合。

本文研究将等带宽元素定位原理与CSC格式存储相结合,不需要根据单元的拓扑结构来形成索引数组;在单元刚度矩阵层面,将单元刚度矩阵的每个非零元素数据直接集成到总体刚度矩阵中。

CSC格式的总体刚度矩阵的集成主要分为两阶段,第一阶段与第二阶段相继作用完成CSC格式的矩阵组装:①将单元刚度矩阵每块矩阵数据无序放入总刚度矩阵一维数组中,形成等带宽辅助数组,并获取总体刚度矩阵数值、列高、每个数据的行号。等带宽辅助数组起定位作用;②根据等带宽一维定位辅助数组将无序的数组进行排列,形成最终的CSC存储格式。

K_II、K_GG、K_IG或K_GI块集成原理相同,以K_II块集成为例,第一阶段组集流程如表4所示。表4中

K I I n u m

数组长度,即非零元素数量;D_i为模型自由度数量;v_alue_K_II[

K I I n u m

]、c_ol.height_K_II[D_i]、r_ow_K_II[

K I I n u m

]为刚度矩阵的数值、行号及列高数组;B_and_K_II[o_ff,N_umber]为等带宽辅助数组(第二阶段组集完成后消除),该数组是一个特殊N_umber的二叉树结构,每一个o_ff数据对应一个编号,可以根据o_ff数据索引出N_umber编号,并且该数组不允许出现重复o_ff数据。

对所有单元遍历完毕即完成K_II块的组装,得到临时数组r_ow_K_II[

K I I n u m

]、v_alue_K_II[

K I I n u m

]、c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber]。第二阶段组集流程如表5所示。R_ow_K_II[

K I I n u m

]、V_alue_K_II[

K I I n u m

]、C_ol.off_K_II[D_i+1]为CSC存储格式的行号、数值、列偏移数组。

对

B a n d_K I I

[o_ff,N_umber]遍历结束即形成按顺序排列的

R o w_K I I

[

K I I n u m

]、

V a l u e_K I I

[

K I I n u m

]、

C o l . o f f_K I I

[D_i+1]数组。对于K_GG、K_IG或K_GI的集成,将上述集成过程的块替换为其他相应子块数据即可,并且在实际组装过程中,所有子块可以同时进行,所以该方法对单元遍历完成之后即可完成总体刚度矩阵的组装,减少了矩阵数据约束处理工作,避免了节点编号优化的难题。

3 算例分析

收起

以北京市轨道交通房山线长稻区间四跨单线梁桥为研究对象,大桥设计桥长120 m,跨径均为30 m,桥面净宽12.5 m,该区间桥梁采用简支箱梁,T型独柱桥墩和钻孔灌注桩(4桩或5桩)基础,桥面铺设单向轨道,桥墩编号从左到右为0#~4#。全桥有限元模型如图6所示,主箱梁、桥墩和桩采用梁单元;桩土作用采用土弹簧;支座、轨道与桥面连接采用两点弹簧单元;支座与箱梁、桩基与承台、轨道与桥面连接弹簧的下端均采用主从约束,目的使其位移保持一致。

为了验证本文算法正确性、实用性,将非零元素计算方法、组集方法与CSC存储方法应用于自行开发的有限元分析程序ADUTS中。开发平台为Visual Studio 2017、计算机配置为CPU为Intel(R) Core(TM) i7-10700、运行内存16 GB。从桩基础X方向输入图7所示的地震动时程,对全桥模型进行动力线性分析,通过比较组集效率、从内存占用和计算效率3个方面,将本文提出的组集方法、存储算法与ID表组集法、变带存储方法进行比较;为确保计算精度,将两种存储方法的计算结果与OpenSees软件进行对比分析。

对于组集效率方面,在不同的节点和自由度数量下,对采用不同组集方法的组集效率进行比较。结果如图8所示。

辅助数组法相较于ID表组集法效率提升30%左右,原因是辅助数组法减少了构建拓扑结构索引数组步骤;然而,随着模型的增大,提升幅度会逐渐降低,这是由于ID表形成时间相较于整个非零元素组集过程所占的时间比重相对较小。

在内存占用方面,采用双精度浮点型保存刚度矩阵数据,每个数据占用8个字节,整型保存行、列数组,每个数据占用4个字节;对于计算效率方面,以采用Cholesky分解法求解整个动力时程分析为例,时间间隔为0.02 s,共持续61.86 s,计算3 093步。分别对比变带宽存储方法与CSC存储方法的内存占用与求解时间,结果如图9和图10所示。

由图9和图10可以看出,CSC存储在内存占用和求解时间显著低于变带存储。这是由于CSC存储只存储矩阵中的非零元素,从而有效减少了所需存储的数据量;对于变带存储,存储的数据量与节点编号有关,若节点编号设置不合理,则会存储大量刚度矩阵中不必要的零元素。而且,在进行Cholesky分解以求解线性方程组的过程中,求解效率与存储的数据量大小密切相关,数据量越大,所需的分解时间越长。

在计算精度方面,提取整个动力分析过程2#墩顶节点整体坐标系X方向位移和各墩顶地震作用下的最大反应与OpenSees进行对比,最大反应由表6给出,位移曲线如图11所示。

通过对比分析,可以发现两种存储方法计算的结果与OpenSees吻合度较高;以OpenSees为基准,墩顶最大反应误差最大值为0.027%,均在合理范围内。

4 结论

收起

(1)经过约束处理后,只有独立自由度、地震动自由度总体矩阵有贡献,可以减少总体矩阵的规模及相应存储需求。

(2)基于自由度与刚度矩阵非零元素之间的关系,确定刚度矩阵非零元素的分布规律,提出二维、三维有限元刚度矩阵非零元素数量计算方法,给出各分块矩阵与总体刚度矩阵非零元素数量计算公式。

(3)通过等带宽定位原理与CSC格式总体矩阵集成相结合,可以在单元层面直接将每个矩阵元素集成到总体矩阵中,组集效率提高30%,避免了节点编号的优化的难题。

(4)与变带宽存储相比,CSC存储空间节省在68%以上;在采用直接法Cholesky分解来求解有限元线性方程组时,采用CSC存储进行求解可以显著提升线性方程组的求解效率,同时保证了计算精度;与变带宽存储相比,效率提升在71%以上,并且节点数量越多,效率提升越明显。

基金

收起

中国博士后科学基金(2024M750780)
中建六局科技研发计划(CSCEC6B-2024-Z-7)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Bakari

A I

, Dahiru

I A

. Comparison of Jacobi and Gauss-Seidel iterative methods for the solution of systems of linear equations[J]. Asian Research Journal of Mathematics, 2018, 8(3): 1-7.

[2]

纪国良, 丁勇, 周曼, 等. 工程计算中大型稀疏矩阵存储方法研究[J]. 数值计算与计算机应用, 2018, 39(3): 217-230.

Guoliang

, Ding

Yong

, Zhou

Man

, et al. Research on the storage method of large-scale sparse matrix in engineering calculation[J]. Journal on Numerical Methods and Computer Applications, 2018, 39(3): 217-230.

[3]

Czech Technical University in Prague, Czech Technical University in Prague. Evaluation of different approaches to solution of the direct solution of large, sparse systems of linear equations[J]. Advanced Materials Research, 2017, 4442: 97-101.

[4]

姜涛, 王安麟, 朱灯林. 有限元结点编号的综合带宽优化算法[J]. 机械设计, 2005, 22(11): 3-6.

Jiang

Tao

, Wang

Anlin

, Zhu

Denglin

. Synthetic bandwidth optimization algorithm of finite element node numbering[J]. Journal of Machine Design, 2005, 22(11): 3-6.

[5]

吴鸿庆. 结构有限元分析[M]. 北京: 中国铁道出版社, 2009.

Hongqing

. Finite element analysis of structure[M]. Beijing: China Railway Publishing House, 2009.

[6]

姚松, 田红旗. 有限元刚度矩阵的压缩存贮及组集[J]. 中南大学学报(自然科学版), 2006, 37(4): 826-830.

Yao

Song

, Tian

Hongqi

. Compression storage scheme and assembly procedure of global stiffness matrix finite element analysis[J]. Journal of Central South University(Science and Technology Edition), 2006, 37(4): 826-830.

[7]

王忠雷, 赵国群, 马新武. 三维有限元刚度矩阵的压缩存储算法[J]. 材料科学与工艺, 2012, 20(2): 96-100, 107.

Wang

Leizhong

, Zhao

Guoqun

, Ma

Xinwu

. Compressed storage algorithm of 3D-FEM stiffness matrix[J]. Materials Science and Technology, 2012, 20(2): 96-100, 107.

[8]

刘尧喜, 唐进元, 周炜, 等. 扩展有限元刚度矩阵的CSR存储实现[J]. 机械强度, 2019, 41(6): 1384-1390.

Liu

Yaoxi

, Tang

Jinyuan

, Zhou

Wei

, et al. CSR storage method of extended finite element stiffness matrix[J]. Journal of Mechanical Strength, 2019, 41(6): 1384-1390.

[9]

周张兰. 基于十字链表与三元组表的稀疏矩阵压缩存储实例研究[J]. 软件导刊, 2017, 16(11): 22-25.

Zhou

Zhanglan

. Study on the compressed storage to sparse matrix based on cross list and tiriple table[J]. Software Guide, 2017, 16(11): 22-25.

[10]

Liu

, Ye

, Ma

. Storage and solving of large sparse matrix linear equations[C]// Fourth International Conference on Computational and Information Sciences. New York: IEEE, 2012: 673-677.

[11]

Guo

D H

, Gropp

. Applications of the streamed storage format forsparse matrix operations[J]. The International Journal of High Performance Computing Applications, 2014, 28(1): 3-12.

[12]

Pei

D M

, Meng

F J

, Wang

H L

. An improved compression and storage algorithm of the stiffness matrix of 3D-FEM based on rider binary classification and negative sign CSR[J]. International Journal of u- and e-Service(Science and Technology Edition), 2015, 8(8): 215-224.

[13]

Xing

L Y

, Wang

Z S

, Ding

Z Z

, et al. An efficient sparse stiffness matrix vector multiplication using compressed sparse row storage format on AMD GPU[J]. Concurrency and Computation: Practice and Experience, 2022, 34(23). DOI: 10.1002/CPE.7186.

[14]

Saak

, Schulze

. Diagonally-addressed matrix nicknack: how to improve SpMV performance[J]. PAMM, 2023, 23(4): e202300228.

[15]

Bian

H D

, Huang

J Q

, Dong

R T

, et al. A simple and efficient storage format for SIMD-accelerated SpMV[J]. Cluster Computing, 2021, 24(4): 1-18.

[16]

刘丽, 陈长波. 带状稀疏矩阵乘法及高效GPU实现[J]. 计算机应用, 2023, 43(12): 3856-3867.

Liu

, Chen

Changbo

. Band sparse matrix multiplication and efficient GPU implementation[J]. Journal of Computer Applications, 2023, 43(12): 3856-3867.

[17]

Jendele

, Červenka

. On the solution of multi-point constraints-application to FE analysis of reinforced concrete structures[J]. Computers & Structures, 2009, 87(15/16): 970-980.

[18]

Boungard

, Wackerfuß

. Consideration of nonlinear multipoint constraints in finite element analyses based on a master-slave eli-mination scheme operating at the global level[J]. PAMM, 2023, 22(1): e202200311.

[19]

张军锋, 温珺博, 张金鹏, 等主从约束在有限元计算中的实现方式和过程[J]. 水利与建筑工程学报, 2023, 21(2): 195-200.

Zhang

Junfeng

, Wen

Junbo

, Zhang

Jinpeng

, et al. Realization method and process of master-slave constraint in finite element method[J]. Journal of Water Resources and Architectural Engineering, 2023, 21(2): 195-200.

[20]

Quyen

V T B

, Tien

D N

, Dung

N N

, et al. Penalty function method for imposing nonlinear multi freedom and multi node constraints in finite element analysis of frame systems[C]// IOP Conference Series: Materials Science and Engineering. Hanoi: IOP Publishing, 2020: 1621-1629.

[21]

Ahn

J G

, YangH

, Kim

J G

. Multipoint constraints with lagrange multiplier for system dynamics and its reduced-order modeling[J]. AIAA Journal, 2020, 58(1): 385-401.

[22]

Tsuta

, Yamaji

. Finite element analysis of contact problem theory and practice in finite element structural analysis[M]. Tokyo: University of Tokyo Press, 1973.

[23]

韩学川, 陶连金, 张宇, 等. 不同类型地震波作用下一体化车站结构地震响应分析[J]. 科学技术与工程, 2021, 21(4): 1506-1514.

Han

Xuechuan

, Tao

Lianjin

, Zhang

, et al. Seismic response of integrated station structures under different types of seismic waves[J]. Science Technology and Engineering, 2021, 21(4): 1506-1514.

[24]

Peng

X Y

, Yu

P C

, Zhu

, et al. Proposal of a coupled DDA-SPH method incorporating a new contact algorithm for soil-structure interaction simulations in geotechnical hazards[J]. Computers and Geotechnics, 2023, 164. DOI: 10.1016/J.COMPGEO.2023.105849.

[25]

钟万勰. 计算结构力学微机程序设计[M]. 北京: 水利电力出版社, 1986.

Zhong

Waixie

. Computer programming for computational structural mechanics[M]. Beijing: Water Conservancy and Electric Power Press, 1986.

[26]

唐松花, 罗迎社, 周筑宝. 线弹性动力学中的最小势能原理(含最小余能原理)[J]. 动力学与控制学报, 2005, 3(1): 36-40.

Tang

Songhua

, Luo

Yingshe

, Zhou

Zhubao

. Minimum potential energy principle in linear elastic dynamics(including minimum complementary energy principle)[J]. Journal of Dynamics and Control, 2005, 3(1): 36-40.

2025年第25卷第10期

PDF下载

346

140

引用本文

BibTeX

文章信息

doi: 10.12404/j.issn.1671-1815.2403703

接收时间：2024-05-19
首发时间：2025-07-09
出版时间：2025-04-08

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-05-19
修回日期：2025-01-13

基金

中国博士后科学基金(2024M750780)

中建六局科技研发计划(CSCEC6B-2024-Z-7)

作者信息

¹ 河南大学建筑工程学院, 开封 475004

² 同济大学桥梁工程系, 上海 200092

通讯作者:

* 李运华(1980—),男,汉族,河南项城人,硕士,高级实验师。研究方向:材料学。E-mail:24561591@qq.com。

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/kxjsygc/CN/10.12404/j.issn.1671-1815.2403703

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

类型	主要形式	处理方式
Dirichlet 边界^[17]	固定约束	直接消去法
指定位移	对角元素改1,其他元素置0
地震动作用	地震作用的施加 (一致激励、非一致激励、人工边界地震动输入)	本质上是一种运动着的边界约束条件,直接代入计算
主从约束	单级约束, (刚性约束、柔性约束)	将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]
多级约束^[18],
时变系数约束
非线性约束	大位移大变形问题,	转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]
曲线或曲面滑动支承,
接触碰撞问题

类型

主要形式

处理方式

Dirichlet
边界^[17]

固定约束

直接消去法

指定位移

对角元素改1,其他元素置0

地震动
作用

地震作用的施加
(一致激励、非一致激励、
人工边界地震动输入)

本质上是一种运动着的边界约束条件,直接代入计算

主从
约束

单级约束,
(刚性约束、柔性约束)

将从节点特性转换到主节点上^[18-19],通常使用罚函数法^[20]、拉格朗日乘子法^[21]

多级约束^[18],

时变系数约束

非线性
约束

大位移大变形问题,

转换为单级主从约束方程;使用罚函数法、拉格朗日乘子法处理等^[22-24]

曲线或曲面滑动支承,

接触碰撞问题

模型	N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
1	1	1	24	23	1	—	—	—
2	2	1	24	23	1	12	12	1
		2	24	22	2

模型

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

—

N	单元	Y_i	D_i	G_i	n_i₍_i_-1)	$n i (i - 1) D$	$n i (i - 1) G$
3	1	3	3	0	4 4	3 3	1 1
2	6	5	1
3	8	6	2

单元

Y_i

D_i

G_i

n_i₍_i_-1)

n i (i - 1) D

n i (i - 1) G

4
4

3
3

1
1

第一阶段:
1.初始化 1-1.初始化K_II块、v_alue_K_II[ $K I I n u m$ ]、r_ow_K_II[ $K I I n u m$ ]、 c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber]; 1-2.临时变量N_um初始化为0。(记录元素位置)。
2.遍历E₀(1,1,N_E) 2-1.获取当前单元E_LEM(E₀)_i; 2-2.获取当前单元独立自由度个数D_i及相应K_II子块; 2-3.判断单元是否循环完毕,若循环完毕,集成结束;若循环没有完毕,进入操作3。
3.遍历K_II子块 3-1.获取子块矩阵元素v_alue、r_ow、c_ol; 3-2.判断r_ow是否小于c_ol,若小于进行3-3,若大于返回3-1; 3-3.计算o_ff=U_BW(i-1)+j-i+1,并判断o_ff是否存在于B_and_K_II[o_ff,N_umber]中,若存在,根据o_ff索引出N_umber,将v_alue_K_II[N_umber]进行数值叠加;若不存在,存储相应数据:v_alue_K_II[N_um]、r_ow_K_II[N_um]=r_ow、c_ol.height_K_II[c_ol]++、 B_and_K_II[o_ff,N_umber](默认第一次遍历不存在,直接存储即可)。
4.单元循环完毕,返回操作操作2

第一阶段:

1.初始化
1-1.初始化K_II块、v_alue_K_II[

K I I n u m

]、r_ow_K_II[

K I I n u m

]、
c_ol.height_K_II[D_i]、B_and_K_II[o_ff,N_umber];
1-2.临时变量N_um初始化为0。(记录元素位置)。

4.单元循环完毕,返回操作操作2

第二阶段:
1.初始化 1-1.临时变量N初始化为0。(记录循环次数);初始化 $R o w_K I I$ [ $K I I n u m$ ]、 $V a l u e_K I I$ [ $K I I n u m$ ]、 $C o l . o f f_K I I$ [D_i+1]。
2.遍历 $B a n d_K I I$ [o_ff,N_umber] 2-1.获取当前o_ff对应的编号N_umber; 2-2.获取 $r o w_K I I$ [N_umber]、 $v a l u e_K I I$ [N_umber]; 2-3. $R o w_K I I$ [N]= $r o w_K I I$ [N_umber]、 $V a l u e_K I I$ [N]= $v a l u e_K I I$ [N_umber];判断N是否小于D_i+1,若小于,将 $c o l . h e i g h t_K I I$ [D_i]的前N个数值进行叠加并赋值给数组 $C o l . o f f_K I I$ [N+1]。 2-4.变量N++,返回操作2-1。

第二阶段:

1.初始化
1-1.临时变量N初始化为0。(记录循环次数);初始化

R o w_K I I

[

K I I n u m

]、

V a l u e_K I I

[

K I I n u m

]、

C o l . o f f_K I I

[D_i+1]。

2.遍历

B a n d_K I I

[o_ff,N_umber]
2-1.获取当前o_ff对应的编号N_umber;
2-2.获取

r o w_K I I

[N_umber]、

v a l u e_K I I

[N_umber];
2-3.

R o w_K I I

[N]=

r o w_K I I

[N_umber]、

V a l u e_K I I

[N]=

v a l u e_K I I

[N_umber];判断N是否小于D_i+1,若小于,将

c o l . h e i g h t_K I I

[D_i]的前N个数值进行叠加并赋值给数组

C o l . o f f_K I I

[N+1]。
2-4.变量N++,返回操作2-1。

位置	OpenSees/ cm	CSC存储/ cm	变带宽存储/cm	误差1/ %	误差2/ %
0#墩顶	1.079 30	1.079 00	1.079 00	0.027 0	0.027 0
1#墩顶	1.079 10	1.079 07	1.079 07	0.002 7	0.002 7
2#墩顶	1.081 44	1.081 41	1.081 41	0.002 7	0.002 7
3#墩顶	1.084 60	1.084 57	1.084 57	0.002 7	0.002 7
4#墩顶	1.102 53	1.102 48	1.102 48	0.004 5	0.004 5

位置

OpenSees/
cm

CSC存储/
cm

变带宽存
储/cm

误差1/
%

误差2/
%

0#墩顶

1.079 30

1.079 00

0.027 0

1#墩顶

1.079 10

1.079 07

0.002 7

2#墩顶

1.081 44

1.081 41

0.002 7

3#墩顶

1.084 60

1.084 57

0.002 7

4#墩顶

1.102 53

1.102 48

0.004 5