粘弹性固体的精细积分有限元算法
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
在流变问题的有限元解析中一般采用时步2粘性初应变法来描述1213也就是把粘性应变看成初应变对每一时步叠加上相应的初应变对于线性流变问题只要连续?断地求解线性弹性方程即可能得到应?应变等随时间的变化情况
第21卷第1期 2004 年 2 月
计算力学学报 Ch inese Journa l of Com puta tiona l M echan ics
(2. 6)
A =
0
I
0
0
0
…0
I
00
0
0
0
…I
- an I - an- 1 I - an- 2 I … - a1I (2. 7)
B = {Β1 I Β2 I … Βn I }T
(2. 8)
C = [I 0 0 … 0] D = Β0 I = b0 I
(2. 9)
x = Ρ - Β0Ε= Ρ - E 0Ε 则 (2. 11) 可写为
关键词: 粘弹性固体本构方程; 状态空间方程; 精细积分有限元 中图分类号: O 241. 8 文献标识码: A
1 引 言
一般而言, 粘弹性固体本构方程数学表达式分 为积分型和微分型, 其数值求解, 主要是时域上的 离散计算[ 123 ]。对积分型的粘弹性体本构方程而言, 若直接采用数值积分, 则为求算由整个应力历史所 描述的积分项, 一般必须保存每个离散时刻的应力 等大量信息, 如果在空间也做大规模的数值离散, 将给计算带来很大困难。对于微分型的粘弹性体本 构方程, 可以通过有限差分而离散为多步递推公 式, 但如微分方程的阶数较高, 则也会出现信息贮 存量过大的问题。 对于某些积分型的本构方程, 可 将其化为阶数不高的微分方程再行求解, 以避免直 接求解积分方程带来的困难。由于粘弹性问题的复 杂性, 求解方法的研究就显得很重要, 许多年来, 人 们一直不断地对它进行探索、改进和完善[429], 以求 更准确、简捷、有效的进行实际问题的数值分析。
式中 Ta 0 = Ta , 因此
[ I + TaN ] = [ I + TaN - 1 ]2 =
[ I + TaN - 2 ]4 = … = [ I + Ta ]2N = T (Σ) (3. 6)
上三式构成了精细积分的基础, 主要应用 2N 算法 的高效性。然而一个关键步是单位矩阵不能参与中 间计算, 以免损失数字精度。
表达式[ I + TaN ] 由 [ I + TaN - 1 ] 计算; [ I + TaN - 1 ] 则由[ I + TaN - 2 ] 计算等。一般的迭代式为
Ta i = 2 × Ta i- 1 + Ta i- 1 × Ta i- 1
( i = 1, 2, …, N ) (3. 7)
将迭代结果与单位矩阵相加形成 T (Σ)。 可以认为非齐次项在时间步 ( tk , tk+ 1) 内是线
2 粘弹性固体本构方程的状态空间 表示及其与微、积分型本构 方程间的转换
不失一般性, 普遍线性粘弹性体的应力、应变 速率关系为[122 ]
Ρ(n) + a1Ρ(n- 1) + … + anΡ = b0Ε(n) + b1Ε(n- 1) + … + bnΕ
选择状态变量[12 ]
(2. 1)
Β0 = b0
Β1 = b1 - a1Β0
Β2 = b2 - a1Β1 - a2Β0
(2. 3)
……
Βn = bn - a1Βn- 1 - … - an- 1Β1 - an Β0
110
计算力学学报
第 21 卷
则
x= Ax + B Ε
(2. 4)
式中
Ρ = Cx + D Ε
(2. 5)
x = {x 1 x 2 … x n}T
收稿日期: 2001209225; 修改稿收到日期: 20022012101 作者简介: 刘靖华3 (19782) , 男, 博士生;
钟万勰 (1934) , 男, 中科院院士 1
步划分为 2N 个精细步长, 其精度可以达到计算机 的精度, 具有理论简单, 编程计算容易, 计算精度高 等优点。最后, 本文由有限元初应变法出发, 形成一 套求解粘弹性问题的精细积分有限元体系。
将特解 vp ( t) 代入动力系统的一般解中得
v ( tk+ 1) = T (Σ) × [ v ( tk ) + H - 1 (Ε0 + H - 1Ε1) ] -
H - 1 [ Ε0 + H - 1Ε1 + Ε1Σ]
(3. 11)
式中 Σ = tk+ 1 - tk 1 文献[ 10 ] 通过数值算例的计算结果和简要的
记方程 (3. 1) 的特解为 vp ( t) , 并暂时假设已经 求得, 当注意到 Σ= 0 时的初始条件时, 方程的一般 解为
v ( t) = T (Σ) [ vh ( tk ) - vp ( tk ) ] + vp ( t) (3. 3)
对于高精度精细积分法, 关键是计算矩阵 T (Σ) , 它 关系到这种方法的效率与特性, 因此精细积分法将 精细划分时间步, 即
Ε= Εe + Εv
(4. 1)
弹性应变 Εe 与应力 Ρ 呈线性关系, 即物理方程 为
Ρ = D Εe = D (Ε- Εv )
其中 D 为弹性矩阵。 几何方程为
(4. 2)
Ε= B ∆e 由式 (4. 3) 可得单元节点力为
(4. 3)
∫ F e = B TD (Ε- Εv ) dv = v ∫B TD (B ∆e - Εv ) dv = v
性的, 即
Ε= Ε0 + Ε1 ( t - tk )
(3. 8)
式中 Ε0, Ε1 分别为与时间无关的向量。动力方程特 解 vp 满足下式
vp = H vp + Ε0 + Ε1 ( t - tk )
(3. 9)
很容易验算特解 vp 可由下式表示
vp = H - 1{- Ε0 + [H - 1 + I ( t - tk ) Ε1 ]} (3. 10)
以图 1 三参数粘弹性固体模型为例来具体说明上 述方程。该模型的微分表示为
Ρ+
E0
+ Γ
E 1Ρ
=
E 0Ε+
E
0E
Γ
1
Ε
(2. 10)
式中 E 0, E 1, Γ 为材料常数。式 (2. 10) 写为更为一 般的表达式为
式中
Ρ+ a1Ρ = b0Ε+ b1Ε
(2. 11)
a1 =
E0
+ Γ
E 1,
∑ F v =
F
e v
,
为粘性应变所产生的总的附加节点荷载
(4. 11)
下面来说明平衡方程 (4. 8) 求解的步骤: (1) 加载瞬间的弹性分析 在 施加载荷的瞬间, 粘性应变尚未发展, 只有
弹性应变。即在初始时刻时 t = 0 时, Εv = 01 于是, 整个系统在施加荷载瞬间的平衡方程
(2. 13)
x= -
E0
+ Γ
E 1x
-
E
2 0
Γ
Ε
Ρ = x + E 0Ε
写成积分形式则为
(2. 14)
Ρ ( t) = exp -
(E 0
+ Γ
E 1)
t
[ Ρ (0) - E 0Ε(0) ] +
∫t exp 0
(E 0
+ Γ
E 1)
(t
-
Σ)
(-
E
2 0
Γ
)
Ε(Σ) dΣ+ E 0Ε( t)
V o l. 2 1 , N o. 1 Feb ruary 2004
文章编号: 100724708 (2004) 0120109206
粘弹性固体的精细积分有限元算法
刘靖华3 1, 范益群2, 钟万勰3
(1. 上海交通大学 建筑工程与力学学院, 上海 200240; 2. 上海隧道工程与轨道交通设计研究院, 上海 200070; 3. 大连理工大学 工程力学系, 辽宁 大连, 116024)
摘 要: 粘弹性固体本构方程的数学表达式分为微分型和积分型两种, 其数值求解主要是时域上离散计算。文中 从微分型表达式出发导出其状态空间方程的数学表达式, 通过严格推导论证了它与微、积分型表达式的等价性; 引入状态空间方程, 从而利用精细积分格式来求解粘弹性固体本构方程; 给出了粘弹性固体本构方程的精细积 分有限元算法, 为求解粘弹性固体本构方程的数值解提供了一个新的途径, 具有计算简便, 求解精度高等优点。
(2. 15)
上式即为粘弹性体本构方程的积分表示。
3 粘弹性体本构方程的精细积分格式
式 (2. 14) 状态空间方程可表示为一般的形式 为
x= H x + Ε Ρ = x + E 0Ε
(3. 1)
对于上述方程, 钟万勰教授提出了一种高精 度的精细积分格式, 本文与此引入[10]。
动力系统 (3. 1) 的通解是满足于方程右端项 为零的齐次方程的解, 即
Hale Waihona Puke b0=E 0, b1 =
E 0E 1
Γ
令
(2. 12)
Β0 = b0 = E 0
Β1 = b1 -
a1Β0 =
E 0E 1
Γ
-
E0
+ Γ
E 1 Β0
=
-
E
2 0
Γ
vh = exp (H Σ) C = T (Σ) C (3. 2)
式中 Σ∈ [ tk , tk+ 1 ] 是积分步长, C 是由初始条件决 定的向量。
=
B TD Εvdv
v
(4. 6)
112
计算力学学报
第 21 卷
而 Ke 为单元刚度矩阵:
∫ Ke = B TDB dv v
(4. 7)
所以由式 (4. 5) 可写出总体平衡方程为
式中
K∆= F + F v
(4. 8)
∑ K = Ke, 为总体刚度矩阵
(4. 9)
∑ F = F e, 为总的节点荷载 (4. 10)
在流变问题的有限元解析中, 一般采用“时步 2 粘性初应变法”来描述[12, 13 ], 也就是把粘性应变 看成初应变, 对每一时步叠加上相应的初应变, 对 于线性流变问题, 只要连续不断地求解线性弹性方 程, 即可能得到应力、应变等随时间的变化情况。下 面来说明一下初应变法的思路与步骤。
总的应变 Ε分为弹性应变 Εe 和粘性应变 Εv , 即
误差分析, 表明了精细积分的高度精确性, 其数值 结果实际上是计算机的精确解, 可以比拟于精确解 的数值结果。
4 有限元法解析的思路与步骤
工程上的流变问题在理论上可以采用L ap lace 变换来求解, 但对于实际问题, 由于拉氏反演的困 难, 还难以用这类解析解来求出计算结果。而采用 有限元数值分析方法, 对求解工程上的复杂情况是 非常有效的。
[ I + H ∃ t + (H ∃ t) 2 2 + (H ∃ t) 3 6 + (H ∃ t) 4 24 ]m =
[ I + Ta ]m
(3. 4)
注意到下式, 即
[ I + Ta i ] = [ I + Ta i- 1 ]2 =
[ I + 2 × Ta i- 1 + Ta i- 1 × Ta i- 1 ] i = 1, 2, …, N (3. 5)
∃t = Σm
式中 m = 2N , 一般取N = 20, 则m = 1048576, 对 于 ∃ t 则是相当的小了, 以至于达到计算机的有效 位数。
第 1 期
刘靖华: 粘弹性固体的精细积分有限元算法
111
应用 T aylo r 展开法, 将矩阵 T (Σ) 展开至四阶 为
T (Σ) = [ exp (H ∃ t) ]m =
本文从粘弹性体本构关系的微分型表达式出 发, 导出其状态空间方程的数学表示。 状态方程的 求解很多, 可以改写成单一变量的微分方程, 然后 按微分方程的正规解法或微分方程的 L ap lace 变 换解法解出, 也可采用矩阵指数函数求解状态方 程, 还可利用状态转移矩阵求得。但以上解法, 计算 繁琐, 精度不是很高。 本文采用钟万勰提出的精细 积分格式[10] 来求解粘弹性本构方程。 精细积分利 用2N 类 的 算 法 , 相 当 于 在 每 一 时 间 步 内 再 进 一
∫ ∫ B TDB dv ∆e - B TD Εvdv =
v
v
K e∆e -
F
e v
(4. 4)
因而, 单元平衡方程为
Ke∆e =
Fe +
F
e v
(4. 5)
因此, 如果把粘性应变 Εv 当作初应变, 式 (4. 4) 中
的
F
e v
即表示由于粘性应变这一初应变而产生的单
元节点荷载, 即
∫ F
e v
K∆0 = F
(4. 12)
由上式解得施加荷载瞬时的位移场 ∆0, 将 ∆0 代入几何方程
Ε0 = B ∆0
(4. 13)
可求得应变场 Ε0, 将 Ε0 代入物理方程
Ρ0 = D Ε0
(4. 14)
可求得初始时刻的应力场 Ρ01 此时, 粘性应变 Εv0 = 01
x1 = Ρ - Β0Ε
x2 = Ρ- Β0Ε- Β1Ε= (x1 - Β1Ε)
x3 = Ρ- Β0Ε- Β1Ε- Β2Ε= (x2 - Β2Ε) ……
x n = Ρ - (n- 1) (xn- 1 -
Β0Ε(n- 1) Βn- 1Ε)
…-
Βn- 2Ε-
Βn- 1Ε= (2. 2)
式中, 待定系数 Β0, Β1, Β2, …Βn 的计算公式为
第21卷第1期 2004 年 2 月
计算力学学报 Ch inese Journa l of Com puta tiona l M echan ics
(2. 6)
A =
0
I
0
0
0
…0
I
00
0
0
0
…I
- an I - an- 1 I - an- 2 I … - a1I (2. 7)
B = {Β1 I Β2 I … Βn I }T
(2. 8)
C = [I 0 0 … 0] D = Β0 I = b0 I
(2. 9)
x = Ρ - Β0Ε= Ρ - E 0Ε 则 (2. 11) 可写为
关键词: 粘弹性固体本构方程; 状态空间方程; 精细积分有限元 中图分类号: O 241. 8 文献标识码: A
1 引 言
一般而言, 粘弹性固体本构方程数学表达式分 为积分型和微分型, 其数值求解, 主要是时域上的 离散计算[ 123 ]。对积分型的粘弹性体本构方程而言, 若直接采用数值积分, 则为求算由整个应力历史所 描述的积分项, 一般必须保存每个离散时刻的应力 等大量信息, 如果在空间也做大规模的数值离散, 将给计算带来很大困难。对于微分型的粘弹性体本 构方程, 可以通过有限差分而离散为多步递推公 式, 但如微分方程的阶数较高, 则也会出现信息贮 存量过大的问题。 对于某些积分型的本构方程, 可 将其化为阶数不高的微分方程再行求解, 以避免直 接求解积分方程带来的困难。由于粘弹性问题的复 杂性, 求解方法的研究就显得很重要, 许多年来, 人 们一直不断地对它进行探索、改进和完善[429], 以求 更准确、简捷、有效的进行实际问题的数值分析。
式中 Ta 0 = Ta , 因此
[ I + TaN ] = [ I + TaN - 1 ]2 =
[ I + TaN - 2 ]4 = … = [ I + Ta ]2N = T (Σ) (3. 6)
上三式构成了精细积分的基础, 主要应用 2N 算法 的高效性。然而一个关键步是单位矩阵不能参与中 间计算, 以免损失数字精度。
表达式[ I + TaN ] 由 [ I + TaN - 1 ] 计算; [ I + TaN - 1 ] 则由[ I + TaN - 2 ] 计算等。一般的迭代式为
Ta i = 2 × Ta i- 1 + Ta i- 1 × Ta i- 1
( i = 1, 2, …, N ) (3. 7)
将迭代结果与单位矩阵相加形成 T (Σ)。 可以认为非齐次项在时间步 ( tk , tk+ 1) 内是线
2 粘弹性固体本构方程的状态空间 表示及其与微、积分型本构 方程间的转换
不失一般性, 普遍线性粘弹性体的应力、应变 速率关系为[122 ]
Ρ(n) + a1Ρ(n- 1) + … + anΡ = b0Ε(n) + b1Ε(n- 1) + … + bnΕ
选择状态变量[12 ]
(2. 1)
Β0 = b0
Β1 = b1 - a1Β0
Β2 = b2 - a1Β1 - a2Β0
(2. 3)
……
Βn = bn - a1Βn- 1 - … - an- 1Β1 - an Β0
110
计算力学学报
第 21 卷
则
x= Ax + B Ε
(2. 4)
式中
Ρ = Cx + D Ε
(2. 5)
x = {x 1 x 2 … x n}T
收稿日期: 2001209225; 修改稿收到日期: 20022012101 作者简介: 刘靖华3 (19782) , 男, 博士生;
钟万勰 (1934) , 男, 中科院院士 1
步划分为 2N 个精细步长, 其精度可以达到计算机 的精度, 具有理论简单, 编程计算容易, 计算精度高 等优点。最后, 本文由有限元初应变法出发, 形成一 套求解粘弹性问题的精细积分有限元体系。
将特解 vp ( t) 代入动力系统的一般解中得
v ( tk+ 1) = T (Σ) × [ v ( tk ) + H - 1 (Ε0 + H - 1Ε1) ] -
H - 1 [ Ε0 + H - 1Ε1 + Ε1Σ]
(3. 11)
式中 Σ = tk+ 1 - tk 1 文献[ 10 ] 通过数值算例的计算结果和简要的
记方程 (3. 1) 的特解为 vp ( t) , 并暂时假设已经 求得, 当注意到 Σ= 0 时的初始条件时, 方程的一般 解为
v ( t) = T (Σ) [ vh ( tk ) - vp ( tk ) ] + vp ( t) (3. 3)
对于高精度精细积分法, 关键是计算矩阵 T (Σ) , 它 关系到这种方法的效率与特性, 因此精细积分法将 精细划分时间步, 即
Ε= Εe + Εv
(4. 1)
弹性应变 Εe 与应力 Ρ 呈线性关系, 即物理方程 为
Ρ = D Εe = D (Ε- Εv )
其中 D 为弹性矩阵。 几何方程为
(4. 2)
Ε= B ∆e 由式 (4. 3) 可得单元节点力为
(4. 3)
∫ F e = B TD (Ε- Εv ) dv = v ∫B TD (B ∆e - Εv ) dv = v
性的, 即
Ε= Ε0 + Ε1 ( t - tk )
(3. 8)
式中 Ε0, Ε1 分别为与时间无关的向量。动力方程特 解 vp 满足下式
vp = H vp + Ε0 + Ε1 ( t - tk )
(3. 9)
很容易验算特解 vp 可由下式表示
vp = H - 1{- Ε0 + [H - 1 + I ( t - tk ) Ε1 ]} (3. 10)
以图 1 三参数粘弹性固体模型为例来具体说明上 述方程。该模型的微分表示为
Ρ+
E0
+ Γ
E 1Ρ
=
E 0Ε+
E
0E
Γ
1
Ε
(2. 10)
式中 E 0, E 1, Γ 为材料常数。式 (2. 10) 写为更为一 般的表达式为
式中
Ρ+ a1Ρ = b0Ε+ b1Ε
(2. 11)
a1 =
E0
+ Γ
E 1,
∑ F v =
F
e v
,
为粘性应变所产生的总的附加节点荷载
(4. 11)
下面来说明平衡方程 (4. 8) 求解的步骤: (1) 加载瞬间的弹性分析 在 施加载荷的瞬间, 粘性应变尚未发展, 只有
弹性应变。即在初始时刻时 t = 0 时, Εv = 01 于是, 整个系统在施加荷载瞬间的平衡方程
(2. 13)
x= -
E0
+ Γ
E 1x
-
E
2 0
Γ
Ε
Ρ = x + E 0Ε
写成积分形式则为
(2. 14)
Ρ ( t) = exp -
(E 0
+ Γ
E 1)
t
[ Ρ (0) - E 0Ε(0) ] +
∫t exp 0
(E 0
+ Γ
E 1)
(t
-
Σ)
(-
E
2 0
Γ
)
Ε(Σ) dΣ+ E 0Ε( t)
V o l. 2 1 , N o. 1 Feb ruary 2004
文章编号: 100724708 (2004) 0120109206
粘弹性固体的精细积分有限元算法
刘靖华3 1, 范益群2, 钟万勰3
(1. 上海交通大学 建筑工程与力学学院, 上海 200240; 2. 上海隧道工程与轨道交通设计研究院, 上海 200070; 3. 大连理工大学 工程力学系, 辽宁 大连, 116024)
摘 要: 粘弹性固体本构方程的数学表达式分为微分型和积分型两种, 其数值求解主要是时域上离散计算。文中 从微分型表达式出发导出其状态空间方程的数学表达式, 通过严格推导论证了它与微、积分型表达式的等价性; 引入状态空间方程, 从而利用精细积分格式来求解粘弹性固体本构方程; 给出了粘弹性固体本构方程的精细积 分有限元算法, 为求解粘弹性固体本构方程的数值解提供了一个新的途径, 具有计算简便, 求解精度高等优点。
(2. 15)
上式即为粘弹性体本构方程的积分表示。
3 粘弹性体本构方程的精细积分格式
式 (2. 14) 状态空间方程可表示为一般的形式 为
x= H x + Ε Ρ = x + E 0Ε
(3. 1)
对于上述方程, 钟万勰教授提出了一种高精 度的精细积分格式, 本文与此引入[10]。
动力系统 (3. 1) 的通解是满足于方程右端项 为零的齐次方程的解, 即
Hale Waihona Puke b0=E 0, b1 =
E 0E 1
Γ
令
(2. 12)
Β0 = b0 = E 0
Β1 = b1 -
a1Β0 =
E 0E 1
Γ
-
E0
+ Γ
E 1 Β0
=
-
E
2 0
Γ
vh = exp (H Σ) C = T (Σ) C (3. 2)
式中 Σ∈ [ tk , tk+ 1 ] 是积分步长, C 是由初始条件决 定的向量。
=
B TD Εvdv
v
(4. 6)
112
计算力学学报
第 21 卷
而 Ke 为单元刚度矩阵:
∫ Ke = B TDB dv v
(4. 7)
所以由式 (4. 5) 可写出总体平衡方程为
式中
K∆= F + F v
(4. 8)
∑ K = Ke, 为总体刚度矩阵
(4. 9)
∑ F = F e, 为总的节点荷载 (4. 10)
在流变问题的有限元解析中, 一般采用“时步 2 粘性初应变法”来描述[12, 13 ], 也就是把粘性应变 看成初应变, 对每一时步叠加上相应的初应变, 对 于线性流变问题, 只要连续不断地求解线性弹性方 程, 即可能得到应力、应变等随时间的变化情况。下 面来说明一下初应变法的思路与步骤。
总的应变 Ε分为弹性应变 Εe 和粘性应变 Εv , 即
误差分析, 表明了精细积分的高度精确性, 其数值 结果实际上是计算机的精确解, 可以比拟于精确解 的数值结果。
4 有限元法解析的思路与步骤
工程上的流变问题在理论上可以采用L ap lace 变换来求解, 但对于实际问题, 由于拉氏反演的困 难, 还难以用这类解析解来求出计算结果。而采用 有限元数值分析方法, 对求解工程上的复杂情况是 非常有效的。
[ I + H ∃ t + (H ∃ t) 2 2 + (H ∃ t) 3 6 + (H ∃ t) 4 24 ]m =
[ I + Ta ]m
(3. 4)
注意到下式, 即
[ I + Ta i ] = [ I + Ta i- 1 ]2 =
[ I + 2 × Ta i- 1 + Ta i- 1 × Ta i- 1 ] i = 1, 2, …, N (3. 5)
∃t = Σm
式中 m = 2N , 一般取N = 20, 则m = 1048576, 对 于 ∃ t 则是相当的小了, 以至于达到计算机的有效 位数。
第 1 期
刘靖华: 粘弹性固体的精细积分有限元算法
111
应用 T aylo r 展开法, 将矩阵 T (Σ) 展开至四阶 为
T (Σ) = [ exp (H ∃ t) ]m =
本文从粘弹性体本构关系的微分型表达式出 发, 导出其状态空间方程的数学表示。 状态方程的 求解很多, 可以改写成单一变量的微分方程, 然后 按微分方程的正规解法或微分方程的 L ap lace 变 换解法解出, 也可采用矩阵指数函数求解状态方 程, 还可利用状态转移矩阵求得。但以上解法, 计算 繁琐, 精度不是很高。 本文采用钟万勰提出的精细 积分格式[10] 来求解粘弹性本构方程。 精细积分利 用2N 类 的 算 法 , 相 当 于 在 每 一 时 间 步 内 再 进 一
∫ ∫ B TDB dv ∆e - B TD Εvdv =
v
v
K e∆e -
F
e v
(4. 4)
因而, 单元平衡方程为
Ke∆e =
Fe +
F
e v
(4. 5)
因此, 如果把粘性应变 Εv 当作初应变, 式 (4. 4) 中
的
F
e v
即表示由于粘性应变这一初应变而产生的单
元节点荷载, 即
∫ F
e v
K∆0 = F
(4. 12)
由上式解得施加荷载瞬时的位移场 ∆0, 将 ∆0 代入几何方程
Ε0 = B ∆0
(4. 13)
可求得应变场 Ε0, 将 Ε0 代入物理方程
Ρ0 = D Ε0
(4. 14)
可求得初始时刻的应力场 Ρ01 此时, 粘性应变 Εv0 = 01
x1 = Ρ - Β0Ε
x2 = Ρ- Β0Ε- Β1Ε= (x1 - Β1Ε)
x3 = Ρ- Β0Ε- Β1Ε- Β2Ε= (x2 - Β2Ε) ……
x n = Ρ - (n- 1) (xn- 1 -
Β0Ε(n- 1) Βn- 1Ε)
…-
Βn- 2Ε-
Βn- 1Ε= (2. 2)
式中, 待定系数 Β0, Β1, Β2, …Βn 的计算公式为