编程.docx

编程.docx

《编程.docx》由会员分享，可在线阅读，更多相关《编程.docx（12页珍藏版）》请在冰豆网上搜索。

编程

有限单元法的渗流计算

一、基本原理和方法

有限单元法是把连续体或研究区域离散化为有限个单元体的集合体来进行研究的。

引

用变分原理或伽勒金法，对研究问题建立模式，推导近似解产生——组方程，最后归结为求

解大量联立方程式的计算。

因此，可以把有限单元法概括为划分单元体来模拟实物进行物

理量分布上的近似．以电子计算机为工具在矩阵分析和近似计算的基础上进行所需精度的

数值计算。

（一）支配方程与边界条件

一船以测压管水头h的分布为研究对象时．作为二向渗流问题考虑（x-y

Y垂直剖面平面），渗流的基本微分方程式及其边界条件如下（见图354）：

（3-120）

初始条件：

水头边界

流量边界

上式为考虑土体压缩的非稳定渗流运动方程，式中的．称为单位贮水量

（尺度1／l）．即单位体积的饱和土体，在下降1个单位水头时，由于土体压缩（ρga）和

水的膨胀（ρgnβ）所释放出来的贮存水量。

这里的α和β分别为土和水的压缩系数;ρg=

为水的容重。

当不考虑压缩性时

则式（3—120）变为:

（3-121）

若渗透系数Ky＝Kx＝K，式（3-121）即变为拉普拉斯方程。

需要指出式（3—121）也就是

稳定渗流的微分方程；

渗流场常遇到两种边界条件。

第一类边界Γ1为已知边界水头值，第二类边界Γ2为已

知或计算出的边界流量值，不透水边界屑第二类的特例，自由面（浸润线）边界，在稳定渗流时

自由面上各点又必须满足水头恒等于垂直坐标高度的一类边界条件：

即

对于上述渗流微分方程的解答，根据变分原理为求下面的泛函取极小值

（3-124）

上式右端的未项为第二类边界积分，经过取泛函极小值后在计算中即自动达成第二类边界条件。

至于第一类边界条件，则在计算中直接赋给已知的边界水头值。

（二）有限单元分析与计算公式

1．应用有限单元法求解，首先要把研究的渗流场划分成有很多个单元，常用者为三角形单元。

其中任意一个三角形单元有i,j,m三个结点．并设单元内水头函数h的分布为下面最简单的线性变化模式:

（3-125）

上式中的三个常数a1、a2、a3可由3个结点i、j、m用的坐标及其相应水头hi、hj、hm来表示，即依次代入上式写出三个方程式可求出常数。

；再代回上式就可得出此典型单元内

任意点的水头，写成矩阵形式为（右上角e指单元而言）

（3-126）

N只是坐标X，Y的函数，称为形状函数，其值为

（3-127）

式中的系数a、b、c和三角形单元的面积Δ都是坐标x,y的函数

（3-128）

其它系数按照i、j、m的次序轮换排列。

三角形单元的面积为:

（3-129）

为了不使面积得出负值，三个角点i、j、m的次序应按逆时针方向填写。

其次求单元内水头函数的导数为：

（3-130）

由式（3—126）．求其对时间的导数，得

（3-131）

将式（3—126）、式（3—130）、式（3—131）以及式（3—122）的各值代入式（3-124）的泛函，并分别对单元三个结点水头求泛函的微商，经过推演可得

（3-132）

这里应当注意．如果是流场的内部单元，则仅有前两项的和；只有在自由面才有最后一项。

式中的系数矩阵为

（3-133）

其中在自由面上代表有效排水孔隙率的矩阵[P]．规定j、m为边界上单元的相邻两点，如图3—54所示把由式（3—132）求得各单元泛函的微商相加．并使其等于零求极小值，就得到整个流场的泛函对各结点水头求导数的方程组。

即

（3-134）

式中，n为结点的总数，

表示对所有的单元求和。

应当指出两点：

（1）因为对特定的某—特定的结点i只有其相连接的结点值出现，而且系数中也只有相邻

的单元才有贡献，因此上式对所有单元求和，实际上就相当于对环绕结点i的各单元求和。

（2）在划分单元的流场中、对于未知结点都存在上式的方程；但是对于任何一个已

知结点，是不能给予变分的，因而就不存在上式的条件。

因此，上式线性方程的数目就等于未知结点的数日n。

至于已知结点值，将作为自由项而存在于各方程中。

将上面汇总的方程组写成矩阵形式，为

（3-135）

式中{F}为形成总矩阵时、由已知结点值所得的常数项的列向量。

用差分表示时间项时，则得：

（3-136）

此式就是对于一般非稳定渗流问题最后要求解的线性代数方程组。

其中总的系数短阵和常

数的列向量中的典型元素，均为对各单元的求和，即

（3-137）

这里Kij，Sij，Pij为总的系数矩阵中第i行第j列的元素；

为各单元中相应

于总坐标编号第i行第j列的元素。

各单元的系数矩阵见式（3—133）。

同样，常数项

也是相应于总坐标编号的。

式中求和项上面的m为单元总数。

k为自由面边界上的单元数（参看图3—54）。

就是说，除与自由面上结点水头有关的矩阵项Pij是对自由面边界单元求和外,其余的为所有单元求和。

由式（3—136）可知，由前一时刻t-Δt的结点水头分布，就可计算该时刻t的结点水头分布；因此只需要知道初始边界条件时的稳定流场并作为初始值就可计算以后边界水位改变条件下各时段流场分布。

而且在计算时，只是重复了自由面上边界单元各结点水头h及h·它们是同一未知数．对解得的自由面边界水头h，还必须用式（3—123）h·=y核算加以修正。

’

稳定渗流是上述计算的一种特例．即矩阵[S]、[P]项都等于零，则求解方程组变为

[K]{h}+{F}=0（3-129）

上面均为与未知结点数n相同数目的线性代数方程式，[K]为n*n的方阵，{h}为未知结

点水头列向量，{F}为常数列向量。

在求得流场的水头分布后，可计算通过任意断面的渗流量。

计算渗流量的断面一船按

单元两条边的中点连线而取为折线断面，即所谓中线法；也可以按单元的边选取，但前者

有更高的精度。

关于中线法计算公式及其推演如下。

图3—56表示经过某一流区的一个断面，其两端与边界Ψ=0和Ψ=q，q为该断面的总渗流量。

该断面是由一排单元中线相连接的折线。

假如通过中点连线的流量为Δq，则

q=Δq（3-140）

取图中ijm作为一个典型三角形单元进行分析，并令中线平行于jm边时．则有

（3-141）

式中υn表示垂直于单元中线a1，a2的流速。

根据达西定律，υx，υy表示主渗透方向的流速时，引用式（3—130）后，得

（3-142）

将上式代入式（3—141）可得典型单元中线a1a2的通过流量为

（3=143）

按照式（3—140）累加各面上各单元即得过流断面的渗流量q。

应当注意，用式（3—143）计算时是具有方向的，对计算断面上的单元角点编号与υn的

正方向均须作好规定，再累加。

二．示例：

如图3—57所示为一个矩形土坝，在上下游固定水位条件下的稳定渗流，划分少数单元用手算来具体说明计算的过程。

[解]稳流场的单元划分如图3—57所示，先假定一自由面．结点编号先编未知结点l

一8，再编已知结点9—13；自由面上的结点和单元可以编成连号也可以不编成连号，但

前者可使形成总矩阵的带宽为最小。

三角形单元和结点的数据如表3—5。

数据列好后，即可计算各三角形单元的矩阵。

例如单元①代入式（3-128）。

按照三结点i，j，m次序轮换算出系数a，b，c

代入式（3—129）计算单元面积：

代入式（3—133）计算单元的渗流矩阵

当Kx=Ky=K时，为

渗流矩阵是的右上角数码代表单元编号，右下角为其结点编号。

同样可计算其余各单单元矩阵如下：

将各单元矩阵中相应编号的系数迭加，或代入式（3—137）．即形成总的渗流矩阵[K]。

例如总矩阵中第1行第l列的系数，第1行第3列的系数和第2行第5列的系数等，应为

其次计算方程组中的常数项。

对于本例来说，只有在计算各单元时涉及到的已知结点，它与水头列向量中已知水头值的乘积所得的常数，例如已知结点9—13的相应水头值为h9=h10=1,h11=h12=h13=6,其中第l行系数与水头相乘后的常数项就有0x6，9x6。

这样将所有相应的常数项迭加．就可得出方程组中相当式（3—139）的列向量{F}：

，例如方程组第一式中的常数项应为

因此求得式（3—139）的矩阵形式为

解上式的方程组就得出各结点的水头脑。

在电子数字汁算机上按照常用方法其得结果

如表3书，并用甘油模型试验结果的自由面位置相比较．可知甚为一致。

从而推知，虽然划

分的流场单元数不多，但只要安排合理，也可得到高精度的计算结果

程序说明

本程序由三个部分构成：

zjbsj.for，SHJ.DAT和zjb.for。

其中SHJ.DAT为数据文件，它是由zjbsj.for程序运行生成的。

zjb.for是稳定渗流的计算程序。

它在运行时需要调用SHJ.DAT数据文件。

针对不同的计算问题，需要对数据文件进行修改。

计算的结果输出在OUT.TET文档中。

一．输入数据

（1）NNNENFIXMK

NN-------结点总数

NE-------单元数

NFIX----已知水头的结点数

MK------渗透系数

（2）LOC（I.1）LOC（I.2）LOC（I.3）

第I个单元逆时针排列的结点号，即i，j，m按单元结点顺序填写

（3）CX（J）CY（J）

CX（J）:

第J个结点X方向的坐标

CY（J）:

第J个结点Y方向的坐标

按结点编号顺序填写

（4）IFIX（K）

第K个已知水头的结点编号

二．给定问题的已知水头的结点修改

（1）在程序中相应的部分要对已知水头的结点进行手动修改

（2）在进行程序的水头输出时需要对程序中的已知水头进行手动修改

三．输出数据（存在OUT.TXT）

（1）NNNENFIXMK

（2）单元结点按逆时针的编号

（3）结点坐标

（4）已知水头的结点编号

（5）各单元的标号及单元系数矩阵

（6）各结点进行高斯消元时的等效计算水头

（7）各结点计算出来的水头输出（包括已知的结点）

（8）各单元通过平行于ij边的中线的流量计算，速度输出

（9）截面的总流量输出

四．注意事项

（1）for程序所在的文件路径一定不能是中文，否则程序运行将会报错。

（2）zjbsj.for程序是用以生成名字为SHJ.DAT的数据文件的。

如果此文件已经生成运行zjbsj.for程序将会出错，显示错误为连接错误，此时因先删除SHJ.DAT数据文件或者将zjbsj.for程序中OPEN（5,FILE='SHJ.DAT',STATUS='NEW'）行的'NEW'改为'UNKNOWN'即可，数据文件（SHJ.DAT）用记事本打开。

（3）单元i，j，m按逆时针顺序填写。

（4）应当注意，用式（3—143）计算时是具有方向的，对计算断面上的单元角点编号与υn的

正方向均须作好规定，再累加。

这就要求在计算的时候要进行人为的控制，保证计算结果在求和时的正确性。

（5）要记的对已知的水头进行单独手动的输入和修改，且在程序中要在两处进行修改。

五．程序的局限性分析

（1）在数据文件已经建立的情况下，对于不同的计算问题，例如单元数目、结点数目、已知的水头大小及数目、等等条件改变时，需要在数据文件（SHJ.DAT）中对相应的部分进行修改。

（2）应当注意，用式（3—143）计算时是具有方向的，对计算断面上的单元角点编号与υn的

正方向均须作好规定，再累加。

程序是基于单元的中线流速进行推导的，这就要求在进行单元划分的时候，截面必须是由单元的中线连接而成，否则程序进无法计算。

（3）单元的编号最好是连续标定，尽量使半带宽最小，从而节省程序运行过程中的存储空间。

实际编号时很难找到使半带宽最小的方案。

且程序在形成总的系数矩阵时并未将已知水头的结点对应的行和列去掉，只是令其系数为一个大数，已知水头的结点的计算结果为0，然后在进行手动的修改，不是太方便。

（4）在计算稳定渗流的程序（zjb.for）中，调用的两个子程序有比较多的语句都是一样的，我们却未能将其统一起来。

六．总结

通过对《深基坑支护工程技术》的学习，我们学会了自己去看一些比较专业的课外书籍和资料。

特别在写论文的过程中，收获不少。

一方面培养了和同学分工合作的团队精神（teamwork）。

由于FORTRAN语言学的时间比较长，程序编写过程中遇到了不少问题。

在解决这些问题的过程中，我们复习了FORTRAN语言，同时也对有限元法有了更清楚的认识，初步学会了将有限元法应用与工程实践的计算。

由于经验不足，编写的程序过多的依赖于钱家欢，殷宗泽老先生编著的教材，其中的实例也是仿造坝体的有限元计算而编写，在应用于计算基坑开挖的时候是有局限性的，在此，我们也恳切的希望罗晓辉老师能批评和帮助改正。

最后，感谢罗晓辉老师给我们上了一门与众不同的课程。

七．参考书目

《土工原理与计算》---------钱家欢，殷宗泽

《深基坑支护工程技术》---------罗晓辉

《有限元法基础与程序设计》---------王元汉，李丽娟，李银平

《FORTRAN语言》-------潭浩强，田淑清