基于ABAQUS二次开发的索单元
总第 236 期 交 通 科 技 Serial No. 236
2009 年第 5 期 Transportation Science & Technology No. 5 Oct. 2009
DOI 1013963/ j1issn1167127570120091051004
收稿日期 :2009205203
基于 ABAQUS 二次开发的索单元
李 昊 简方梁 曹一山
(同济大学土木工程学院 上海 200092)
摘 要 基于大型通用有限元分析软件 ABAQUS 计算平台 ,利用该软件提供的用户单元子程序
U EL 接口 ,以悬链线单元为基础 ,推导了迭代分析公式 ,向 ABAQUS 中添入悬链线单元的刚度矩
阵及等效荷载列阵 ,扩充 ABAQUS 单元库 ,使悬链线索单元用于 ABAQUS 有限元分析 ,并采用
Fortran 语言开发了接口程序 ,研究结果表明 ,所研制的接口程序开发思路正确 ,计算精度满足要求。
关键词 U EL 二次开发 索单元 非线性分析
通用有限元软件 ABAQUS 具有完备的前后
处理功能 ,丰富的单元库 ,大量的材料模型库以及
强大的求解器 ,在工程实践和理论研究中得到了
广泛的应用。但是作为通用软件 ,不免在某些专
业领域内有所欠缺 ,因此 ABAQUS 提供了大量
的用户子程序接口 ,它允许用户结合专业问题定
义符合自己问题的模型。本文以悬链线索单元为
基础 ,结合 ABAQUS 的用户自定义单元接口 ,编
制 Fortran 计算程序 ,模拟大跨度悬索结构悬索、
斜拉桥拉索进行有限元分析。相对于研制和开发
大型专业软件 ,该方法可以节约大量时间 ,为今后
分析软件的开发提供一条新的途径。
1 ABAQUS 用户自定义单元子程序
ABAQUS 具有良好的开放性 ,为了满足不
同专业分析的需要 ,ABAQUS 提供了用户子程
序接口生成非标准分析程序来满足用户的需要 ,
它允许用户通过子程序以代码的形势来扩展主程
序的功能 ,并给用户提供了强大而又灵活的用户
子程序接口和应用程序接口 ,这些用户子程序接
口使得用户解决一些问题时有很大的灵活性 ,同
时大大地扩充了 ABAQUS 的功能。
ABAQUS 的 U EL 是 ABAQUS 为用户提供
的自定义单元工具接口 ,它可以使用户根据自己
的需要 ,向 ABAQUS 程序内部嵌入线性和非线
性有限单元 ,值得强调的是 ,ABAQUS 程序允许
在一个 U EL 中实现多个用户单元 ,同时 ,这些单
元适用于 ABAQUS 的大部分求解程序。
ABAQUS 单元用户子程序主要包括以下几
部分 :
(1) ABAQUS 约定的子程序题名说明 :
SUBROU TIN E U EL (RHS ,AMA TNX , ?,
PERIOD) 。
(2) ABAQUS 定义的参数声明。
(3) 用户自定义的局部变量声明。
(4) 根据用 户自身 需求所 定义的 程序代
码段。
ABAQUS 主程序调用 U EL 子程序 ,在增量
步开始时 ,通过接口将数据由主程序进入 U EL ,
将变量值赋予 U EL 的相应变量 ,在增量步结束
后 ,变量将会被更新 ,并通过 U EL 的接口传入
ABAQUS 主程序[122 ] 。
在 U EL 中 ,L FLA GS 数组是 ABAQUS 的
默认数组 ,它的值决定了用户单元子程序需定义
的变量 ,在一般的非线性分析步骤中 ,在当前增量
步的开始 ,包含有这些定义变量值的数组就被传
递到 U EL 中 ,除非程序中要求 U EL 不做更新 ,
否则在增量步结束时 ,这些定义的变量值应当被
更新。本 文 根 据 悬 链 线 索 单 元 公 式 , 编 制 了
ABAQUS 用 户 单 元 子 程 序 , 该 单 元 可 与
ABAQUS 其他单元一起使用 ,适用于悬索桥、斜
拉桥拉索的精细分析。
2 悬链线索单元
2. 1 基本假定
(1) 索在弹性阶段工作。
(2) 大位移小应变。
(3) 索是理想柔性的 ,只承受拉力而不能承
受弯矩和压力。
(4) 假设索的几何形状为悬链线。
2. 2 位移模式的确定
索单元 局 部 坐 标 系 o2x y z 和 整 体 坐 标 系
O2X Y Z如图 1 所示。
图 1 索单元的坐标系
索的两端点为 A 、B , A B 连线所在直线为 x
轴 , u、v、w 为索截面沿 x 、y 、z 轴的结点位移 , u1 、
v1 、w1 和 u2 、v2 、w2 分别为 A 、B 两点的节点位
移 ,取索单元的位移模式如下[3 ] :
u = (1 - x
L ) u1 + x
L u2
u = (1 - x
L ) v1 + x
L v2
u = (1 - x
L ) w1 + x
L w 2 - Δz
(1)
局部坐标系中 ,索在均布荷载 q 作用下的悬
链线方程为[4 ] :
z = - 1
2β coshβ- cosh (2βx
L - β) (2)
式中 :L 为索的弦长;β= ql
2 H ,其中 : l 为索的水平
长度 ; q 为均布荷载集度 ; H 为索的水平张力。推
导得
Δz = D
2β2 [coshβ- cosh (2βx
L - β) -
βsinhβ+ (2 x
L - 1) βsinh (2βx
L - β) ]ΔL
(3)
式中 : D = β(β- 2sinhβ+ sinhβcoshβ)
2 (βcoshβ- sinhβ) (4)
ΔL = ( u2 + u1 ) +
( u2 - u1 ) 2 + ( v2 - v1 ) 2 + ( w2 - w1 ) 2
2L (5)
因此 ,索单元的位移模式为 :
u = Nue (6)
式中 : ue 为节点位移列阵;
N =
1 - x
L 0 0 x
L 0 0
0 1 - x
L 0 0 x
L 0
Φ 0 1 - x
L Φ 0 x
L
(7)
其中 : Φ = D
2β2 (coshβ- cosh (2βx
L - β) -
βsinhβ+ (2 x
L - 1)βsinh (2βx
L - β) ) (8)
2. 3 单元刚度矩阵
采用上述位移模式 ,得索单元的切线刚度矩
阵 :
Ke = K0 + Kσ (9)
式中 : K0 为弹性刚度阵; Kσ 为几何刚度阵。
K0 = EA
L
α1 0 0 - α1 0 0
0 0 0 0 0 0
0 0 α2 0 0 - α2
- α1 0 0 α1 0 0
0 0 0 0 0 0
0 0 - α2 0 0 α2
(10)
Kσ = N
L
α3 0 0 - α3 0 0
0 1 0 0 - 1 0
0 0 α2 0 0 - 1
- α3 0 0 α3 0 0
0 - 1 0 0 1 0
0 0 - 1 0 0 1
(11)
式中 :α1 = 1 + D sinh 2β
4β2 - cosh 2β
2β +
D2 ( - 1
24 + sinh 4β
32β - cosh 4β
64β2 + sinh 4β
256β3 ) (12)
α2 = sinh 2β
4β - 1
2 (13)
α3 = D2 ( 1
6 + sinh 2β- 2βcosh 2β+ 2β2 sinh 2β
8β3 )
(14)
要得到刚度矩阵的表达式 ,必须先求得β,β
满足约束方程 :
f (β) = c2 + lsinhβ
β
2 - s0 - ql2
4 EAβ
cothβ
β sinh2 β+ 2 βc
l
2
+ 1 = 0
采用牛顿迭代法对β进行求解 ,迭代方程为 :
βn+1 = βn - f (β)
f′(β) (15)
112009 年第 5 期 李 昊等 :基于 ABAQUS 二次开发的索单元
由 matlab 得出 f′(β) 表达式为 :
f′(β) =
l2 coshβsinhβ
β 2 - l2 - sinh2 β
β 3
c2 + lsinhβ
β
2
+ ql2
4 EAβ2 1 +
cothβ 2c2β2
l2 + sinh2 β
β
- ql2
4 EAβ
cothβ(4c2β
l2 + 2cothβsinhβ)
β -
cothβ(2c2β
l2 + sinh2 β)
β2 -
csch2 β(2c2β2
l2 + sinh2 β)
β
(16)
2. 4 质量矩阵
在对索进行动力分析时 ,需要建立质量矩阵。
质量矩阵采用集中质量矩阵形式 :
M =
1
2 ρA S 0 0 0 0 0 0
0 1
2 ρA S 0 0 0 0 0
0 0 1
2 ρA S 0 0 0 0
0 0 0 1
2 ρA S 0 0 0
0 0 0 0 1
2 ρA S 0 0
0 0 0 0 0 1
2 ρA S 0
(17)
3 算例
例 1 由 5 个单元组成的两端铰接的索杆结
构 ,高 5 m ,长 10 m ,弹性模量 E = 200 GPa ,截面
积 A = 1. 963 ×10 - 5 m2 ,密度ρ = 7 800 kg/ m3 ,
初张力 F = 1 000 kN 。6 个节点号码依次为 1～
6 ,5 个单元号码依次为 1～5 ,如图 2 所示 ,计算自
由振动的频率和周期 ,结果见表 1 所列。
图 2 由 5 个单元组成的索杆结构(单位 :m)
表 1 用户索单元和梁单元计算频率比较 cycles/ s
振动模态 梁单元 索单元 振动模态 梁单元 索单元
1 112. 48 112. 42 7 294. 48 295. 52
2 112. 48 112. 42 8 294. 48 346. 65
3 213. 95 213. 84 9 294. 48 346. 65
4 222. 55 213. 84 10 346. 19 476. 15
5 222. 75 251. 38 11 422. 37 654. 11
6 222. 75 295. 53 12 423. 69 769. 66
计算结果表明 ,从 11 阶模态开始 ,由于梁的
弯曲刚度影响 ,梁单元的振动周期降低 ,而频率高
于索单元[5 ] 。
例2 图3为索在自重作用下的线型 ,弹性模量
E = 190 GPa ,截面积 A = 0. 85 cm2 ,自重集度 q =
3. 16 N/ m ,计算其在一集中荷载 P = 8 kN 作用下 ,
B 点相对于恒载状况的位移。相关研究见表 2。
图 3 例 2 计算模型 (单位 :cm)
表 2 索在集中力作用下的位移 cm
作用点位移 文献[5 ] 文献[6] 文献[7] ABAQUS 结果
竖向位移 - 17. 951 - 18. 226 - 18. 456 - 18. 463
水平位移 - 2. 772 - 2. 804 - 2. 819 - 2. 825
4 结语
本文依据 ABAQUS 提供的用户单元子程序
接口 ,推导了悬链线索单元迭代分析公式 ,把悬链
线索单元模型加入 ABAQUS ,扩大了 ABAQUS
在拉索有限元分析中的应用范围 ,并为有限元程
序的开发提供了一个新的并且非常实用的途径。
由于开发过程是基于 ABAQUS 平台 ,所设计的
开发工作量要比独立开发有限元程序少得多 ,有
利于提高开发周期。
参考文献
[ 1 ] Sorensen lnc , ABAQUS/ Standard User Manual
[ M]. Hibbitt , Karlsson , Sorensen Ine. 2002.
[2 ] 叶志才 ,徐 磊 ,王 朝. 基于 ABAQUS 的三维锚
杆单元开发[J ]. 三峡大学学报 ,2008 ,30 (5) :29232.
[3 ] 杨孟刚 ,陈政清. 两节点曲线索元精细分析的非线
性有限元法[J ]. 工程力学 ,2003 ,20 (1) :42247.
[4 ] 张其林 ,预应力结构非线性分析的索单元理论[J ].
工程力学. 1993 ,10 (4) :932101.
[5 ] 庄 茁 ,张 帆 ,岑 松 ,等. ABAQUS 非线性有限
元分析与实例[ M ]. 北京 :科学出版社 ,2005.
[6 ] Knudson W C. Static and dynamic analysis of cable
net structure[D ]. Berkeleg ,Califomia University of
California ,19711
[7 ] 罗喜恒 ,肖汝诚 ,项海帆. 基于精确解析解的索单元
[J ]. 同济大学学报 ,2005 ,30 (4) 4452450.
21 李 昊等 :基于 ABAQUS 二次开发的索单元 2009 年第 5 期
总第 236 期 交 通 科 技 Serial No. 236
2009 年第 5 期 Transportation Science & Technology No. 5 Oct. 2009
DOI 1013963/ j1issn1167127570120091051005
收稿日期 :2009205205
自平衡试验等效转换的附加压缩量修正
马远刚1 ,2 陈开利2 杨春和1
(1. 中国科学院武汉岩土力学研究所 武汉 430071 ; 2. 中铁大桥局集团武汉桥梁科学研究院有限公司 武汉 430034)
摘 要 提出一种新的自平衡试验数据简化等效转换方法 ,它考虑了自平衡加载方式下底层荷载
箱以上桩体发生的压缩量及其转换后桩顶加载方式下的压缩量 ,对等效转换的桩顶位移进行附加
压缩量修正 ,并给出了单层及双层荷载箱的计算公式 ;算例对比分析表明 ,这种方法比其他他简化
转换方法确定的极限承载更准确 ,而且有更高的等效转换曲线相似性。
关键词 自平衡 等效转换 压缩量 极限承载力 相关系数
自平衡加载设备、位移测试手段等在工程实
践中不断完善 ,特别是试验数据桩顶等效转换方
法 ,尤其受到众多学者的关注 ,从理论及实用方面
进行了深入研究。自平衡试验数据等效转换的目
的是要便于与传统静载试验曲线对比 ,并结合桩
顶位移合理确定试桩极限承载力 ,因此等效转换
方法要从其受力机理入手保证等效转换曲线的精
度。归纳而言可分为简化转换法和精确转换
法[123 ] ,由于简化转换法只需要测试荷载箱向上、
向下位移和桩顶位移 ,而无需在桩身预先埋设应
力测试元件 ,现场施工简便 ,便于推广应用。
目前应用的简化转化方法 ,虽然在刚性体假
设的基础上考虑了上段桩的压缩量 ,但压缩量的
计算方式存在一定的误差 ,从而影响基桩承载力
测试的精度。本文从基桩等效转换过程中受力和
变形进行分析 ,提出一种新的简化转换方法。
1 改进简化转换法
将自平衡试验实测分段荷载2位移曲线等效
转 换 为 桩 顶 加 载 方 式 下 的 单 一 荷 载2位 移 曲
线[126] ,就是以底层荷载箱以下桩体的实测荷载2
位移曲线为基础 ,考虑底层荷载箱以上桩体的有
效阻力及附加压缩量 ,叠加等效转换为桩顶加载
方式下的荷载2位移曲线。
1. 1 计算假定及初步转换
对于单层(或双层)荷载箱而言 ,荷载箱将整桩
划分为 2(或 3) 个桩段 (包括桩端) ,简化转换法就
是将这些桩段视为若干个单元 ,并做出以下合理而
偏于保守的假定 : ①对于单层荷载箱 (图 1 中
的 A荷载箱) ,当荷载箱置于桩端 (或附近) 时 ,自
Cable Element Based on ABAQUS
L i Hao , J ian Fangliang , Cao Yishan
(1. School of Civil Engineering Tongji University ,Shanghai 200092 ,China)
Abstract : Based On the user2element subroutine interface provided by ABAQUS ,a catenary cable ele2
ment is established. The catenary cable element' s stiffness matrix and equivalent load vector are in2
duced and corresponding interface program is developed , in which Fortran language is adopted.
ABAQUS element library are expanded. The feasibility and validity of the interface subroutine has al2
so been proved.
Key words : U EL ; redevelopment ; cable element ; nonlinear analysis
基于ABAQUS二次开发的索单元.pdf
