非线性无网格伽辽金法的实现
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
司建辉, 李九红, 简 政
( 西安理工大学 水利水电学院 , 陕西 西安 710048)
摘
要 : 提出了应用无网格伽辽金法计算非 线性混凝土问 题的基本方 法 . 无 网格伽辽金 法 ( EFGM ) 是近些年
发展 起来的一种数值算法 , 它采用移动的最小二乘法构造 形函数 , 从能量泛函的弱变 分形式中得 到控制方 程 , 该法 只需节点信息 , 不需将节点连成单元 . 在积分 网格 中 , 取高斯 点的 本构关 系随 应力 变化来 反映 混凝土 的非 线性性 质 . 文中混凝土的本构模型选用 OT T OSEN 本构 , 屈服准则选 用修正的莫 尔库仑强 度准则 , 通过算例 分析 , 验证了 程序的可靠性及应用无网格伽辽金法解决非线性混凝土问题的可行性 , 表明该方法在混凝 土材料领域 有着广阔的 应用前景 . 关 键 词 : 无网格伽辽金法 ; 本构模型 ; 移动最小 二乘法 ; 非线性 文献标识码 : A 中图分类号 : O 241
m
Lx ^ u( x ) =
j= 1
p j ( x ) aj ( x ^ ) = p T ( x ) a( x ^ ) ( 3)
1
无网格伽辽金法
无网格伽辽金方法与有限元方法最大的不同在
最小移动二乘法要求近似函数与真实函数在各 已知点上的取 值之离 差的 加权 2 范 数最 小. 取权 wi(x ^ ) 为 i 节点的权函数在 x ^ 点的取值 , 则要求:
2
上式中: E 0 Ec D
5
算例分析
为了验证理论及程序的可靠性 , 以如图 1 所示
对下降段影响很大; Ef
受集中力作用的悬臂梁为例 , 进行分析 .
以下公式计算: Ef =
E0 J2 1 式中 : A = E c , x = ( f c ) f 3
图 1 受集中力作 用的悬臂梁
j= 1 n - 1 ji
( 11) ( 12) 所以可得 : ( 13a ) ( 13b ) (x) = ( (x) = L
x
ni ( x ) = A (x) =
i= 1
(x) , y( x), N( x )
( x))T u
( 21b) ( 22) ( 23a)
wi(x) p(x i)pT( xi) wi(x) p(x i), , w n( x ) p ( x n )
p j ( x ) A- 1 ( x ^ ) B( x ^)
( 9)
d( x )
( 20)
现在使 x ^ 成为域内任一点, 即 x ^ 与 x 不加区别 , 则 可得到: Gu( x ) = L x u( x ) 即全域近似可以写成 :
n
( 10)
( 21a)
Gu ( x ) =
i= 1 m
ni ( x ) u i p j ( x ) A ( x ) B( x )
[ 1~ 4]
划分单元和单元联结信息, 而后者则在单元基础上 对节点进行插值. 下面简述一下无网格伽辽金方法 的实施过程. 移动最小二乘法 移动最小二乘法是用加权最小二乘法来近似场 函数的一种 方法. 设场函数为 u ( x ) , 其中 x = ( x , 1. 1 y ) 为场点 , 给定 u( x ) 在 n 个节点上的值: u( x i ) = u i , i = 1 , 2, , n ( 1) 则使用移动最小二乘法可以确定 u( x ) 的一个近似 函数 :
Ni ( x ) = ( 8)
ji
ni ( x ^ , x ) ui
i= 1 m
N ( x ) 称为形函数矩阵, n i ( x ) 可由式 ( 12 ) 计算 . 域内任一点应变为 : (x ) = L 其中 L 为偏导数矩阵: x L= 0 y 0 y x
xy
其中
ni ( x ^ , x) =
j= 1
- 1 T
T
为拉氏乘子, H , H 是一维和零维的 sobolev 空间 . ( 6a ) ( 6b ) ( 6c ) ( 7) 其中 : d( x ) = ( u( x ) , v( x ) ) T N( x ) = [ N 1 ( x ) , N 2 ( x ) , , N n( x ) ] , i = 1, 2, ,n ni ( x ) 0 0 ni ( x ) ( 19a) ( 19b) ( 19c) 根据 ML S 可得域内任意点近似位移为: d( x ) = N ( x ) u ( 18)
广, 如果将材料常数即弹性模量 E 和泊松比
为常数, 而是确定为随应力状态而变化的参数 , 则这 种关系就变为非线性关系了 . 本文针对混凝土材料 的特点 , 采用割线形式的本构关系, 选用 Ot t osen 本 构模型. Ot t osen 建 议的本构模型, 关键是 要明确 3 个 条件 , 详见文献 [ 5] 所述. 本文选取的屈服准则为修正的莫尔库仑强度准 则, 采用的非线性指标为双向应力状态下非线性指 标, 弹性模量及泊松比的计算详见文献 [ 5] 所述, 即 : E s = 1 E0 - ( 1 E0 - E c) 2 2 [ 1 1 2 2 E 0 - ( E 0 - E f ) ] + E c [ D( 1 - ) - 1 ] 2 2 混凝土初始弹性模量; 混凝土应力达到 f c 时的割线模量; 系数, 对 曲线上升段影响不大, 而 非线性指标. 三轴压 缩破坏割线模 量, 取值可 按 Ec 1 + 4( A - 1) x 0. J 2 / f c) f 是 当 <
B( x ^) = 所以
, w n( x ^ ) p ( x n)
u = ( u 1 , u2 ,
a( x ^ ) = A (x ^ ) B( x ^ )u L x^ u ( x ) = p ( x ) A ( x ^ ) B( x ^ )u =
n T - 1
将式 ( 7) 带入式 ( 3) 可得:
第 S2 期
司建辉 等 : 非线性无网格伽辽金法的实现
T T T
t
47 bd T
u
最小 , 所以 J = 2 wi(x ^ ) p ( x i ) p T ( x i ) a( x ^ ) - ui = 0 a i= 1 即 式中 A( x ^ ) 为m 维向量,
n n
(
s
): d -
td -
T
u
n
于插值函数的 构造上 , 前 者通过 移动 最小 二乘 法
收稿日期 : 2005 09 10 作者简介 : 司建辉 ( 1976 ) , 男 , 硕士 , 讲师 , 现从事结构分析研究 .
J a( x ^) =
i= 1
wi(x ^ ) p ( x i ) a( x ^ ) - ui
T
2
( 4)
E mail: s jhf r@ xaut . edu. cn
第 51 卷 第 S2 期 2005 年 12 月
武汉大学 学报( 理学版 ) J. W uhan U niv. ( N at. Sci. Ed. )
V o l. 51 N o. S2 Dec. 2005, 046~ 050
Hale Waihona Puke Baidu
文章编号 : 1671 8836( 2005) S2 0046 05
非线性无网格伽辽金法的实现
1. 2
无网格伽辽金法控制方程 考察弹性力学中的二维边值问题:
ij , j ij
( x ) = D B( x ) (x),
y
( 24) ( 25)
(x),
xy
(x))
+ fi = 0
在 在 在
t u
内 上 上
i
( 16a ) ( 16b ) ( 16c )
其中 , D 为弹性矩阵 . 从而可以得到刚度矩阵 : K= B
式中 a( x ) 为 m 维系数向量 , p T ( x ) 为 m 维基向量. 为了保证收敛性, p ( x ) 应取完全多项式, 本文选用 线性基, 即: p ( x ) = ( 1, x , y ) . 为了求出 Gu( x ) , 先固 定一点 x ^ (x ^ 称为 估值 点) , 求 x ^ 微小邻域内 u ( x ) 的局部近似函数 :
迭代法中比较简单的一种, 又称直接迭代法 , 迭代过 程如文献 [ 6] 所述 .
4
应用无网格伽辽金法解决非线性混凝土 问题的计算流程
2
混凝土的本构模型
由于材料的线弹性关系非常简单, 并且应用较
应用无网格伽辽金法解决非线性混凝土问题的 计算流程如下 : ① 根据研究物体形状生成节点并计算节点影 响半径; ② 各网格所包含的节点号 ; ③ 定高斯点的位置及积分权重 ; ④ 高斯点上逐个计算 ; ⑤ 检索高斯点的影响节点 ; ⑥ 根据应力修改各点 E , ; ⑦ 检索该高斯点对整体平衡方程的贡献并集 入整体平衡方程; ⑧ 解平衡方程; ⑨ 求解应力应变等所需变量; 根据荷载等级 , 返回④ ; 循环所有荷载等级 , 计算最终应力应变等 所需变量 . 不取
x
, Bn ( x ) ]
0 ni ( x ) y ni ( x ) x u
T
( 23b)
ni, k =
j= 1
p j , k [ A- 1 B] j i + p j [ A , k B + A B, k ] ji
-1 - 1
( 14) ( 15)
1 - 1 A,k = - A A , kA - 1
t
为
( 27) N td
T
T
位移边界 . 与方程( 16 ) 对应的能量泛函的弱变分形式为 :
实施过程中, 在域
内布置无网格节 点, 然后
48
武汉大学 学报( 理学版 )
第 51 卷
划分积分网格积分网格是背景网格 , 它与无网格节 点无关, 仅用来完成 数值积分, 得到 K, F 后 , 即可 得到方程 ( 27) 的数值解 , 再由式 ( 22) 和 ( 24) 拟合出 域 内任意点的位移和应变, 进而求得应力.
m T
Gu( x ) =
j= 1
p j ( x ) aj ( x ) = p T ( x ) a( x )
( 2)
.
目前 , 关于无网格伽辽金法的研究主要集中在 线弹性材料上, 对于结构工程中主要面对的混凝土 材料所表现出的非线性本构模型则很少涉及 , 本文 主要介绍如何应用无网格伽辽金法结合非线性本构 模型解决混凝土材料的非线性问题 , 结合算例进行 分析 , 取得了较好的成果 .
( M L S) 对区域内互不相关的节点进行插值 , 不需要
0
引
言
移 动 最 小 二 乘 法 ( M oving L east Square ( M LS) ) 可以通过几个 互不相关节点上的数据 , 拟 合出一个函数, 该函数光滑性好且导数连续, 在求解 边值问题偏微分方程数值解时, 采用 M LS 构 造位 移函数, 这种位移函数的形成及区域积分可以脱离 单元的概念 . 这 一求 解过 程被 称 为 DEM ( Dif fuse Elem ent M ethod, 扩散单元法 ) , 部分学者对此作了 进一步改进, 使得 DEM 求解精度更高 , 应用更为广 泛, 这种改 进后的方 法称为无 网格伽辽 金法 ( Ele m ent F ree Galerkin Met hod( EF GM ) )
( u - u) d H 1, H0
T 0
d = 0, ( 17)
A( x ^ ) a( x ^ ) = B( x ^)u m 矩阵, B( x ^ )为 m
( 5) n 矩阵, u 为 n 式中
s
T
代表
1
的对称部分, u, 为试验函数,
A( x ^) =
i= 1
wi(x ^ ) p(xi) p ( xi) w i(x ^ )p (xi), un )
u = B( x )
B( x ) =
其中 w i ( x ) 为 i 节点的权函数在 x = ( x , y ) T 点的取 值, n i ( x ) 为 i 节点的形函数 x = ( x , y ) T 点的取值 . 由( 11) 式可得到形函数 ni ( x ) 关于坐标的偏导 数为 :
m
其中 B( x ) 为几何矩阵: B( x ) = [ B 1 ( x ) , B 2 ( x ) , ni ( x ) x Bi ( x ) = 0 ni ( x ) y 则任一点应力为: ( x) = D ( x) = (
T
n j - ti = 0
u i = ui 其中 量.
ij i
D Bd
( 26)
为应力张量分量, f 为体力分量 , t 为面力分
t
量, u i 为边界位移 分量, n j 为 边界法向 单位向量 分 为弹性力学问题计算域, 为应力边界 ,
u
弹性力学中的二维边值问题求解方程为: Ku = F 其中 F= N fd +
( 西安理工大学 水利水电学院 , 陕西 西安 710048)
摘
要 : 提出了应用无网格伽辽金法计算非 线性混凝土问 题的基本方 法 . 无 网格伽辽金 法 ( EFGM ) 是近些年
发展 起来的一种数值算法 , 它采用移动的最小二乘法构造 形函数 , 从能量泛函的弱变 分形式中得 到控制方 程 , 该法 只需节点信息 , 不需将节点连成单元 . 在积分 网格 中 , 取高斯 点的 本构关 系随 应力 变化来 反映 混凝土 的非 线性性 质 . 文中混凝土的本构模型选用 OT T OSEN 本构 , 屈服准则选 用修正的莫 尔库仑强 度准则 , 通过算例 分析 , 验证了 程序的可靠性及应用无网格伽辽金法解决非线性混凝土问题的可行性 , 表明该方法在混凝 土材料领域 有着广阔的 应用前景 . 关 键 词 : 无网格伽辽金法 ; 本构模型 ; 移动最小 二乘法 ; 非线性 文献标识码 : A 中图分类号 : O 241
m
Lx ^ u( x ) =
j= 1
p j ( x ) aj ( x ^ ) = p T ( x ) a( x ^ ) ( 3)
1
无网格伽辽金法
无网格伽辽金方法与有限元方法最大的不同在
最小移动二乘法要求近似函数与真实函数在各 已知点上的取 值之离 差的 加权 2 范 数最 小. 取权 wi(x ^ ) 为 i 节点的权函数在 x ^ 点的取值 , 则要求:
2
上式中: E 0 Ec D
5
算例分析
为了验证理论及程序的可靠性 , 以如图 1 所示
对下降段影响很大; Ef
受集中力作用的悬臂梁为例 , 进行分析 .
以下公式计算: Ef =
E0 J2 1 式中 : A = E c , x = ( f c ) f 3
图 1 受集中力作 用的悬臂梁
j= 1 n - 1 ji
( 11) ( 12) 所以可得 : ( 13a ) ( 13b ) (x) = ( (x) = L
x
ni ( x ) = A (x) =
i= 1
(x) , y( x), N( x )
( x))T u
( 21b) ( 22) ( 23a)
wi(x) p(x i)pT( xi) wi(x) p(x i), , w n( x ) p ( x n )
p j ( x ) A- 1 ( x ^ ) B( x ^)
( 9)
d( x )
( 20)
现在使 x ^ 成为域内任一点, 即 x ^ 与 x 不加区别 , 则 可得到: Gu( x ) = L x u( x ) 即全域近似可以写成 :
n
( 10)
( 21a)
Gu ( x ) =
i= 1 m
ni ( x ) u i p j ( x ) A ( x ) B( x )
[ 1~ 4]
划分单元和单元联结信息, 而后者则在单元基础上 对节点进行插值. 下面简述一下无网格伽辽金方法 的实施过程. 移动最小二乘法 移动最小二乘法是用加权最小二乘法来近似场 函数的一种 方法. 设场函数为 u ( x ) , 其中 x = ( x , 1. 1 y ) 为场点 , 给定 u( x ) 在 n 个节点上的值: u( x i ) = u i , i = 1 , 2, , n ( 1) 则使用移动最小二乘法可以确定 u( x ) 的一个近似 函数 :
Ni ( x ) = ( 8)
ji
ni ( x ^ , x ) ui
i= 1 m
N ( x ) 称为形函数矩阵, n i ( x ) 可由式 ( 12 ) 计算 . 域内任一点应变为 : (x ) = L 其中 L 为偏导数矩阵: x L= 0 y 0 y x
xy
其中
ni ( x ^ , x) =
j= 1
- 1 T
T
为拉氏乘子, H , H 是一维和零维的 sobolev 空间 . ( 6a ) ( 6b ) ( 6c ) ( 7) 其中 : d( x ) = ( u( x ) , v( x ) ) T N( x ) = [ N 1 ( x ) , N 2 ( x ) , , N n( x ) ] , i = 1, 2, ,n ni ( x ) 0 0 ni ( x ) ( 19a) ( 19b) ( 19c) 根据 ML S 可得域内任意点近似位移为: d( x ) = N ( x ) u ( 18)
广, 如果将材料常数即弹性模量 E 和泊松比
为常数, 而是确定为随应力状态而变化的参数 , 则这 种关系就变为非线性关系了 . 本文针对混凝土材料 的特点 , 采用割线形式的本构关系, 选用 Ot t osen 本 构模型. Ot t osen 建 议的本构模型, 关键是 要明确 3 个 条件 , 详见文献 [ 5] 所述. 本文选取的屈服准则为修正的莫尔库仑强度准 则, 采用的非线性指标为双向应力状态下非线性指 标, 弹性模量及泊松比的计算详见文献 [ 5] 所述, 即 : E s = 1 E0 - ( 1 E0 - E c) 2 2 [ 1 1 2 2 E 0 - ( E 0 - E f ) ] + E c [ D( 1 - ) - 1 ] 2 2 混凝土初始弹性模量; 混凝土应力达到 f c 时的割线模量; 系数, 对 曲线上升段影响不大, 而 非线性指标. 三轴压 缩破坏割线模 量, 取值可 按 Ec 1 + 4( A - 1) x 0. J 2 / f c) f 是 当 <
B( x ^) = 所以
, w n( x ^ ) p ( x n)
u = ( u 1 , u2 ,
a( x ^ ) = A (x ^ ) B( x ^ )u L x^ u ( x ) = p ( x ) A ( x ^ ) B( x ^ )u =
n T - 1
将式 ( 7) 带入式 ( 3) 可得:
第 S2 期
司建辉 等 : 非线性无网格伽辽金法的实现
T T T
t
47 bd T
u
最小 , 所以 J = 2 wi(x ^ ) p ( x i ) p T ( x i ) a( x ^ ) - ui = 0 a i= 1 即 式中 A( x ^ ) 为m 维向量,
n n
(
s
): d -
td -
T
u
n
于插值函数的 构造上 , 前 者通过 移动 最小 二乘 法
收稿日期 : 2005 09 10 作者简介 : 司建辉 ( 1976 ) , 男 , 硕士 , 讲师 , 现从事结构分析研究 .
J a( x ^) =
i= 1
wi(x ^ ) p ( x i ) a( x ^ ) - ui
T
2
( 4)
E mail: s jhf r@ xaut . edu. cn
第 51 卷 第 S2 期 2005 年 12 月
武汉大学 学报( 理学版 ) J. W uhan U niv. ( N at. Sci. Ed. )
V o l. 51 N o. S2 Dec. 2005, 046~ 050
Hale Waihona Puke Baidu
文章编号 : 1671 8836( 2005) S2 0046 05
非线性无网格伽辽金法的实现
1. 2
无网格伽辽金法控制方程 考察弹性力学中的二维边值问题:
ij , j ij
( x ) = D B( x ) (x),
y
( 24) ( 25)
(x),
xy
(x))
+ fi = 0
在 在 在
t u
内 上 上
i
( 16a ) ( 16b ) ( 16c )
其中 , D 为弹性矩阵 . 从而可以得到刚度矩阵 : K= B
式中 a( x ) 为 m 维系数向量 , p T ( x ) 为 m 维基向量. 为了保证收敛性, p ( x ) 应取完全多项式, 本文选用 线性基, 即: p ( x ) = ( 1, x , y ) . 为了求出 Gu( x ) , 先固 定一点 x ^ (x ^ 称为 估值 点) , 求 x ^ 微小邻域内 u ( x ) 的局部近似函数 :
迭代法中比较简单的一种, 又称直接迭代法 , 迭代过 程如文献 [ 6] 所述 .
4
应用无网格伽辽金法解决非线性混凝土 问题的计算流程
2
混凝土的本构模型
由于材料的线弹性关系非常简单, 并且应用较
应用无网格伽辽金法解决非线性混凝土问题的 计算流程如下 : ① 根据研究物体形状生成节点并计算节点影 响半径; ② 各网格所包含的节点号 ; ③ 定高斯点的位置及积分权重 ; ④ 高斯点上逐个计算 ; ⑤ 检索高斯点的影响节点 ; ⑥ 根据应力修改各点 E , ; ⑦ 检索该高斯点对整体平衡方程的贡献并集 入整体平衡方程; ⑧ 解平衡方程; ⑨ 求解应力应变等所需变量; 根据荷载等级 , 返回④ ; 循环所有荷载等级 , 计算最终应力应变等 所需变量 . 不取
x
, Bn ( x ) ]
0 ni ( x ) y ni ( x ) x u
T
( 23b)
ni, k =
j= 1
p j , k [ A- 1 B] j i + p j [ A , k B + A B, k ] ji
-1 - 1
( 14) ( 15)
1 - 1 A,k = - A A , kA - 1
t
为
( 27) N td
T
T
位移边界 . 与方程( 16 ) 对应的能量泛函的弱变分形式为 :
实施过程中, 在域
内布置无网格节 点, 然后
48
武汉大学 学报( 理学版 )
第 51 卷
划分积分网格积分网格是背景网格 , 它与无网格节 点无关, 仅用来完成 数值积分, 得到 K, F 后 , 即可 得到方程 ( 27) 的数值解 , 再由式 ( 22) 和 ( 24) 拟合出 域 内任意点的位移和应变, 进而求得应力.
m T
Gu( x ) =
j= 1
p j ( x ) aj ( x ) = p T ( x ) a( x )
( 2)
.
目前 , 关于无网格伽辽金法的研究主要集中在 线弹性材料上, 对于结构工程中主要面对的混凝土 材料所表现出的非线性本构模型则很少涉及 , 本文 主要介绍如何应用无网格伽辽金法结合非线性本构 模型解决混凝土材料的非线性问题 , 结合算例进行 分析 , 取得了较好的成果 .
( M L S) 对区域内互不相关的节点进行插值 , 不需要
0
引
言
移 动 最 小 二 乘 法 ( M oving L east Square ( M LS) ) 可以通过几个 互不相关节点上的数据 , 拟 合出一个函数, 该函数光滑性好且导数连续, 在求解 边值问题偏微分方程数值解时, 采用 M LS 构 造位 移函数, 这种位移函数的形成及区域积分可以脱离 单元的概念 . 这 一求 解过 程被 称 为 DEM ( Dif fuse Elem ent M ethod, 扩散单元法 ) , 部分学者对此作了 进一步改进, 使得 DEM 求解精度更高 , 应用更为广 泛, 这种改 进后的方 法称为无 网格伽辽 金法 ( Ele m ent F ree Galerkin Met hod( EF GM ) )
( u - u) d H 1, H0
T 0
d = 0, ( 17)
A( x ^ ) a( x ^ ) = B( x ^)u m 矩阵, B( x ^ )为 m
( 5) n 矩阵, u 为 n 式中
s
T
代表
1
的对称部分, u, 为试验函数,
A( x ^) =
i= 1
wi(x ^ ) p(xi) p ( xi) w i(x ^ )p (xi), un )
u = B( x )
B( x ) =
其中 w i ( x ) 为 i 节点的权函数在 x = ( x , y ) T 点的取 值, n i ( x ) 为 i 节点的形函数 x = ( x , y ) T 点的取值 . 由( 11) 式可得到形函数 ni ( x ) 关于坐标的偏导 数为 :
m
其中 B( x ) 为几何矩阵: B( x ) = [ B 1 ( x ) , B 2 ( x ) , ni ( x ) x Bi ( x ) = 0 ni ( x ) y 则任一点应力为: ( x) = D ( x) = (
T
n j - ti = 0
u i = ui 其中 量.
ij i
D Bd
( 26)
为应力张量分量, f 为体力分量 , t 为面力分
t
量, u i 为边界位移 分量, n j 为 边界法向 单位向量 分 为弹性力学问题计算域, 为应力边界 ,
u
弹性力学中的二维边值问题求解方程为: Ku = F 其中 F= N fd +