有限元固结分析中若干问题探讨_805603882

第二届全国土力学教学研讨会1有限元固结分析中若干问题探讨介玉新1, 2（1 清华大学水利水电工程系，北京，100084； 2 水沙科学与水利水电工程国家重点实验室(清华大学)，北京 100084）摘 要：在用基于比奥固结理论的有限元法进行固结分析时，涉及平衡方程和连续方程的有限元格式构造，荷载、位移和水头、流量等边界条件的处理，以及静孔压和超静孔压的计算等问题，本文对这些问题进行澄清，重点探讨有限元格式中以水头为未知数和以孔压为未知数两者之间的关系。

对计算中所涉及的土的本构模型，尤其邓肯-张模型中关于泊松比的限制条件也进行了初步分析。

关键词：固结；有限元；比奥固结理论；邓肯-张模型中图分类号：O302, O211 文献标识码：A 文章编号：xxxxxxxxxxxxxxxxxxxx作者简介：介玉新(1970－ )，男，博士，副教授，主要从事土工数值分析、基础工程、土工合成材料等方面的教学和科研工作。

E-mail: jieyx@On analysis of consolidation by finite element methodJie Y uxin 1, 2(1. Department of Hydraulic Engineering, Tsinghua University, Beijing, 100084;2. State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing, 100084)Abstract : Consolidation analysis with finite element method (FEM) involves constructing the finite element schemes ofequilibrium equation and continuity equation, and processing such boundary conditions as force, displacement, water head and discharge. Total pore water pressure and excess pore water pressure are also puzzling. In this paper, the author investigates these problems, discusses the relationship between the finite element schemes with water head as unknown and with pore water pressure as unknown. For the constitutive model concerned, the author also probes into the constraints on Poisson ’s ratio.Key words : consolidation; finite element method (FEM); Biot ’s consolidation theory; Duncan-Chang ’s model0 引 言比奥固结方程是土力学数值计算的基本方程。

有限元法是求解比奥固结方程的强有力的数学工具。

比奥固结方程可以看作是平衡方程和连续方程的联立，以平面应变为例，可以写为：1a ）或1b ）其中p 为孔隙水压力，按弹性力学符号规定，压应力为负，与常规土力学规定的符号相反；h h w γ=，h 为________________________ 基金项目：国家重点基础研究发展计划资助(2007CB714102) 收稿日期：2水头。

在渗透系数不随坐标位置变化时，上式中的连续方程也可以写为01s v 2222=∂∂+∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-t S t y k x k wy x w γεγ （2a ） 或01s v 2222=∂∂-∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂tp S t y p k x p k w y x w γεγ （2b ） 利用物理方程和几何方程，可以把式（1）进一步表示成位移分量和孔压（或水头）表示的形式[1~5]，这是位移法有限元格式构造的基础。

求解上述偏微分方程组，除了物理方程、几何方程外，还需要与平衡方程有关的荷载和位移边界条件，以及与渗流有关的水头（孔压）和流量边界条件。

太沙基-伦杜立克理论也可以进行固结分析，但它假定主应力之和在固结过程中不随时间变化，显然是不严密的。

它的优点是能够把孔压计算和应力应变计算分离开来，减少计算工作量。

这在计算机技术发展早期还是有一定意义的，但在今天已没有太大意义。

1 求解比奥固结方程的有限元格式有限元格式可以用里兹法推导，也可以用伽辽金法推导。

全量形式的有限元格式为⎪⎭⎪⎬⎫=---=-e h q e e p e e s e eT c e hF e e c e e R h K h K K R K K }{}{][}{][}{][}{{][}{][ δδ（3a ）或⎪⎭⎪⎬⎫=--=+e p q e e p e e s e eT c epF eec eeR p K p K K R p K K }{}{][}{][}{][}{}{][}{][ δδ （3b ）其中[]A B D B K T e d ]][[][⎰⎰= []A N MB K Tec ]d [}{][⎰⎰=A B k B K s T s e s d ]][[][][⎰⎰=A N N SK T ws e p d ][][][⎰⎰=γ ⎰⎰⎰⎰⎰-+=A y M B A F N s F N R w T T Te h F d }{][d }{][d }{][}{γs N R n Te h q d ][}{⎰=υ⎰⎰⎰+=A F N s F N R T Te p F d }{][d }{][}{A kB s N R w Ts n Te pq d 0][][d ][}{⎰⎰⎰⎭⎬⎫⎩⎨⎧--=γυ{}{}TM 011={}{}T y xF FF =为面力，{}{}Ty xF F F =为体力。

对四节点四边形等参元来说，⎥⎦⎤⎢⎣⎡=432143210000][N N N N N N N N N []4321][N N N N N =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=y N yN y N y N x N x N x N x N B s 43214321][其他符号与一般应力应变分析和渗流分析的符号含义相同。

在时段t n-1～t n 之间，增量形式的有限元格式为()⎪⎭⎪⎬⎫+∆=∆+∆-∆-∆=∆-∆-en e s e h q e e p e s e eT c e hF e e c e e h K R t h K K t K R h K K 1}{][}{}){][][(}{][}{}{][}{][θδδ （4a ）或()⎪⎭⎪⎬⎫+∆=∆+∆-∆∆=∆+∆-e n e s e p q e e p e s e eT c e pF e e c e e p K R t p K K t K R p K K 1}{][}{}){][][(}{][}{}{][}{][θδδ （4b ）其中10<<θ。

5.0<θ时计算不稳定，θ宜取为0.5～1之间，通常可取3/2=θ。

另外，上式中⎰⎰⎰+==AF N s F N R R TT e p F e h F d }{][d }{][}{}{∆∆∆∆ （4c ）在以水头为未知数的有限元格式中，荷载项eh F R }{中的A y M B w Td }{][⎰⎰γ在相减中约去。

注意虽然式（4）中eK ][与式（3）中的相应矩阵形式相同，但此时eK ][指的是切线矩阵，式（3）中e K ][对应的是割线矩阵。

3由于式（4b ）中ep q R }{是以全量形式出现的，所以A k B w T s d 0][][⎰⎰⎭⎬⎫⎩⎨⎧γ项仍然保留在式（4b ）连续方程的右边项中。

在用式（4b ）进行计算时，忽略这一项会在流量边界和y k 不相等的单元体相邻边上产生虚假流量。

多出这一项的原因是在流量边界上水头对y 的一次微分变为孔压表达时需要考虑位置水头的影响。

流量边界n y x m l υυυ=+，用孔压表示为m k l y p k l x p k y n y x w +=⎪⎪⎭⎫ ⎝⎛⋅∂∂+⋅∂∂υγ1 考虑m k y 后，ep q R }{可写为s m k N R y n T e p q d )(][}{⎰+-=υ利用微积分的有关知识，可以把关于m k y 的边界积分变换为面积分，即⎰⎰⎰∂∂=A yN k s m k N iyy i d d 进而可以证明⎭⎬⎫⎩⎨⎧=∂∂w T s T y k B N y k γ0][][][ Sandhu 和Wilson 曾采用组合单元进行固结分析[6]。

位移采用6节点二次形函数，孔压采用3节点一次形函数。

这样位移精度高于孔压精度，而应力精度则可以与孔压精度相当。

但作者认为这种做法会增加编程和网格剖分的难度，必要性不是很大。

因为既然不考虑固结的固体力学有限元计算中应力的精度能够保证，在固结分析中也就没有必要为了孔压的关系提高应力计算的精度。

2 孔压增量与水头增量的关系（4a ）和式（4b ）其实是完全等价的，差别只在一个正负号上，这就意味着应当有()()e n e s e p q en e s e h q p K R t h K R t 11}{][}{}{][}{--+∆-=+∆即ep q e h q e n e s e n e s R R p K h K }{}{}{][}{][11--=+-- （5a ）A kB y A B k B w Ts ew s Ts d 0][][}{]d ][[][⎰⎰⎰⎰⎭⎬⎫⎩⎨⎧=γγ（5b ） 其中{}T w w w w e w y y y y y 4321}{γγγγγ=将⎭⎬⎫⎩⎨⎧w γ0写为 ()()⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y y x y w w γγ 并利用∑==41i i i y N y 进行运算，不难证明式（5b ）确实是成立的。

3 超静孔压与总孔压计算的关系一些读者常常在超静孔压和静孔压的计算中困惑。

从式（4a ）和式（4b ）可以看出，固结计算进行的先决条件是计算初时刻连续方程的右边项{}0h或{}0p 应首先赋值，否则无法进行下一步计算。

对于计算初时刻存在一稳定渗流场的情况，如果计算初时刻给{}0h或{}0p 赋以渗流场的水头或孔压，则是按总的水头或孔压进行计算；如果以0赋值，则是按超静孔压进行计算。