引入接触单元模拟桩土共同作用

表 1 地基土物理力学性质指标统计表 Table 1 Statistic table of foundation soil’s indexes of physico-mechanical properties
层号
土层 名称
土层厚度 层底标高
/m
/m
ω
γ
%
/kN·m-3
主要物理力学性质指标平均值
α1-2
2 计算模型
由于桩土的几何形状、边界条件和载荷都对称
于桩的轴线，因此，将桩土相互作用三维实体模型 简化为二维轴对称模型，其计算简图如图 1 所示。
接触是一种非线性行为，本文通过采用参变量 变分原理及基于此原理的相应的参变量二次规划法 [3]，对接触单元刚度矩阵进行推导。
接触单元 桩体
土体
图 1 桩体-接触单元-土体分析模式 Fig.1 Analysis mode of pile-contact element-soil
Abstract: The contact element model is introduced to simulate the pile-soil’ s reciprocity,then the paper derived the contact element stiffness matrices by adopting parameter variable’s variational principle and parameter variable’s second programming method.With comparison of computing result and test in situation ,it shows that the contact element model is suited to simulating the pile-soil’ s shear force transmission and relative displacement. Keywords: contact element; parameter variable’s variational principle; the pile-soil’ s reciprocity
2.1 接触单元 接触单元是 Okamoto 与 Nakazawa 提出的一种
收稿日期：2003-10-30
修改稿收到日期：2003-12-27
基金项目：上海高校优秀青年教师后备人选资金资助(03YQHB162)
作者简介：齐良锋，女，1973 年生，博士，讲师，从事结构工程中高层建筑上部结构与基础地基的共同作用。
第1期
齐良锋等：引入接触单元模拟桩土共同作用
129
min. ∏[λ(•)] = 1 uˆΤ kuˆ − uˆ(Φλ + Pˆ)
(13)
2
Cuˆ −Uλ − d +υ = 0⎫
υΤ λ = 0,υ ,λ
≥
0
⎬ ⎭
(14)
式中 k, Φ , Pˆ, C, U , d 的矩阵的组装与普通有限单
元法相同。
λkυk = 0 υk , λk ≥ 0 k =1，2，3
(7)
将 fk 作一阶 Taylor 展开，可以将状态方程写成另一 种表达形式：
f
0 k
+ ωkεc
−
mk λ
+υk
=
0；
λk ⋅υk = 0 υk 百度文库 λk ≥ 0 k =1，2，3
(8)
其中
f
0 k
为起始时
fk
值，而
mk
=
⎡ ⎢ ⎣
∂f k ∂Pc
摘 要：在桩与桩侧土界面引入接触单元模型，并采用参变量变分原理及基于此原理的参变量二次规划法，对接触单元刚度
矩阵进行了推导。模型计算结果与实测比较表明，接触单元能较好地模拟桩与土之间的剪力传递和相对位移。
关 键 词：接触单元；参变量变分原理；桩土共同作用
中图分类号：O 343.3
文献标识码：A
The computing model of adopting contact element to simulate the pile-soil’ s reciprocity
—
—
—
—
⑥ 粉质粘土 11.80 369.30
24.3
20.1
0.690
0.31
0.17
10.1
⑦ 粉质粘土 6.60
362.70
24.0
20.0
0.688
0.26
0.14
15.0
b
中砂
2.50
360.20
—
—
—
—
—
—
⑧ 粉质粘土 8.70
351.50
23.5
20.3
0.665
0.23
0.14
13.4
128
岩土力学
2005 年
模拟接触的单元，它根据接触条件，把接触点对（在
接触面上坐标相同的节点）的位移和接触力，以单
元的形式进行表示，可以直接向刚度矩阵中组装，
形成的总刚可以进行向接触面上的凝聚，得到在接
触点的经过缩聚的刚度矩阵。
接触面上一点的接触力可分解为法向接触力
Pn 和切向接触力 Pτ ，即 Pc = {Pτ ,Pn }Τ ，接触条件为：
−
u
(2) n
+δ
*
= =
∆uτ ∆u n
⎫
+
δ
*
⎬ ⎭
(2)
式中
u(a) n
，
u(a) τ
分别为桩土域相对接触点局部坐标
系下的切向、法向位移；δ * 为桩土间的初始间隙（如
打入桩存在过盈配合， δ * < 0 ），也可以表示为：
εc
=
ε
e c
+
ε
p c
，ε
e c
为弹性相对位移，即发生接触但未
达到滑动时的相对位移；ε
(12)
此问题可以通过参数二次规划法求解，具体解 法参见文献[3]。
将接触物体进行单元划分，在接触边界用接触 单元划分，接触点对的主位移与通常的节点位移一 起编排在总体位移向量 uˆ 中，而接触点对的相对位 移编排在总体位移向量的最前端，离散化后平面弹性 接触问题有限元二次规划求解方程的提法归纳如下：
g1 g2
= Pτ + c0 = −Pτ + c0
⎪⎫ ⎬
(3)
g3 = Pn ⎪⎭
c0 为任意常数，则滑动相对位移可表示为：
∑ ε
p c
=
3
λk
k =1
∂g k ∂Pc
，式中 ∂gk ∂Pc
体现相对滑动位移方向，
λk 代表滑动量大小。位移与接触力的关系条件可表 示为：
β (ε n ) =Pnε=n [1E−n βsi(gεnn()ε n )]/ 2⎭⎬⎫
第 26 卷第 1 期 2005 年 1 月
文章编号：1000－7598－(2005) 01―0127―04
岩土力学 Rock and Soil Mechanics
Vol.26 No.1 Jan. 2005
引入接触单元模拟桩土共同作用
齐良锋 1 ，简 浩 2 ，唐丽云 3
(1. 上海应用技术学院 土木与建筑工程系，上海 200235；2.上海岩土工程勘察设计研究院有限公司，上海 200002； 3. 西安科技大学 建筑工程系，陕西 西安 710054）
⎬
(19)
∑8
u = N iui
i=1
⎪ ⎪ ⎪
∑ w =
8 i=1
N i wi
⎪ ⎪⎭
单元刚度矩阵表达式：
[k] e= ∫∫A[B]Τ[D] [B]dA
(20)
2.3 材料特性
在桩土体系中，桩土的弹性模量相差很大，有
2 个量级左右，桩受荷后，一般处于弹性阶段，而
其周边土容易达到塑性状态。因此，在桩土相互作
f
f1
2
==−PPτ τ++µµPPn≤n ≤00⎪⎬⎫
(1)
f3 = Pn≤0 ⎪⎭
µ 为桩土间的摩擦系数，在接触力空间的几何图形
见图 2 所示。当 fk < 0 时，接触力状态点处于 fk = 0 所围区域之内，虽发生接触，但未滑动。当 fk = 0 时， 状态点处于区域边界，桩土作接触滑动或脱开
p c
为滑动相对位移，定义
滑动相对位移
ε
p c
如下：当状态点落在
f1
=
0 时，对
应滑动量为 λ1 ， λ1 > 0 ；当状态点落在 f2 = 0 时，
对应滑动量为 −λ2 ，λ2 > 0 ；当状态点落在 f3 = 0 时，
物体间相互滑脱，对应滑脱量为 λ3 ， λ3 > 0 ；对应
fk ，k =1, 2, 3 定义滑动势函数：
⑨ 粉质粘土 13.40 338.10
24.2
20.0
0.682
0.32
0.14
12.5
(4)
Pτ
=
⎩⎨⎧−
Eτ ετ µPn sign(ε
n
)
(5)
Pτ < −µPn Pτ = −µPn
(6)
其中， En , Eτ → ∞ ，为惩罚因子。 引入原约束松弛变量υk ,令： fk + υk = 0 , k =1，2，3 则接触问题的状态方程为：
f k(∆uτ , ∆un ,λk ) + υk = 0 ;
（ f3 = 0 时）。
εs
f1=0
Pr
b
εs
Pn
b f3=0
f2=0 εs
图 2 滑动条件示意图 Fig.2 Abridged general view of sliding condition
接触相对位移 εc = {ε τ ,ε n}Τ 为
εn
ετ = uτ(1) − uτ(2)
=
u (1) n
⎤Τ ∂g
⎥ ⎦
Dc ( ∂Pc
)Τ
(9)
Dc
=
⎡ Eτ ⎢⎣ 0
0⎤ E n ⎥⎦
(10)
参变量最小势能原理：在所有满足应变-位移关
系式和几何边界条件式的可能位移场中，真实解使
总势能泛函 ∏[λ(•)] 在接触系统状态方程的控制下
取总体最小值：
∏[λ(•)]
=
∫Ω
1ε 2
Τ DεdΩ
−
⎡ ⎢⎣
∫Ω
e
IL
/MPa
E1-2 /MPa
③
黄土
9.20
394.10
28.4
18.7
0.848
0.81
0.11
17.8
④
古黄土
4.70
389.40
27.6
19.6
0.765
0.58
0.20
9.4
⑤ 粉质粘土 5.90
383.50
24.2
20.0
0.692
0.41
0.13
12.1
a
中粗砂
2.40
381.10
—
—
me
=U
e
=
⎢− ⎢⎣
Eτ 0
Eτ 0
µEn En
⎥ ⎥⎦
(18)
2.2 桩、土单元模式 桩、土部分采用平面 8 节点 4 边形轴对称等参
单元，圆柱坐标为整体坐标，r 表示径向坐标，z 表 示轴向坐标，坐标变换和位移模式如下：
8
r = ∑ Niri
⎫ ⎪
i=1
⎪
∑8
z = N i zi
i=1
⎪ ⎪⎪
有限元分析平面点接触单元的刚度矩阵和约束
矩阵为
kce
=
N
e c
Τ
Dce
N
e c
=
⎡Eτ ⎢⎣ 0
0⎤ En ⎥⎦
(15)
Φ
c e
=
N
e c
Τ
ReΤ
=
⎡Eτ ⎢⎣ 0
− Eτ 0
0⎤ En ⎥⎦
(16)
⎡ Eτ µEn ⎤
Cec
=
ωe
N
e c
=
⎢− ⎢⎣
Eτ 0
µEn En
⎥ ⎥⎦
(17)
⎡ Eτ − Eτ µEn ⎤
bΤudΩ
+
∫P Sp
Τ
ud s
⎤ ⎥⎦
+
∫Sc
(
1ε 2
c
Τ
Dcε
c
−
λΤ Rε c )dS
(11)
式中
R
=
(
∂g ∂Pc
)Dc
为常数矩阵；位移
u
为自变量；λ
为不直接参加变分的参变量。因此，接触问题的求
解化为：min. ∏[λ(•)]
s.t.
f (uc , υΤλ =
λ) + 0,υ,
υ λ
=0⎫ ≥ 0⎭⎬
1 引言
随着计算机的发展，有限元在桩土共同作用分 析中得到了愈来愈多的应用，其一个突出的优点就 是通过对土体区域的划分，解决土体分布的层状性 问题。同济大学沈伟跃、赵锡宏[1]应用有限元法对 桩土进行了分析，得出了一些满意的结果，但这种 分析方法认为，桩侧和桩周土之间没有相对位移。 大量工程实践表明，桩侧和桩周土之间不仅能传递 摩阻力，而且它们之间还有较大的相对位移。文献 [2]在有限元分析的基础上，在桩与土之间引入摩擦 单元（节理单元），本文在桩与土之间引入另一种单 元--接触单元以模拟桩与土之间的剪力传递和相对 位移。
QI Liang-feng1 , Jian Hao2 , TANG Li-yun3
(1. Department of Civil Eng., Shanghai Institute of Technology, Shanghai 200235, China ; 2. Shanghai Geotechnical Engineering Investigation & Design Institute Co.Ltd. ,Shanghai 200002, China ; 3. Department of Civil Eng, Xi’an University of Science and Technology, Xi’an 710054, China)
用分析中，桩采用弹性材料，土采用理想弹塑性材
料，采用 Drucker-Prager 屈服准则，使用相关流动
准则，其屈服面并不随着土的逐渐屈服而改变，因
而没有强化准则。
3 实例分析
陕西信息大厦主楼试桩桩长 83 m，地质资料和 物理力学性质指标见表 1，用本文有限元计算模型 进行计算，其桩侧摩阻力分布如图 3，实测试桩截 面轴力传递[4]如图 4 所示