计算机内插在土工分析中的应用

1计算机内插在土工分析中的应用周逢君（深圳市工勘岩土工程有限公司）摘要本文对岩土工程的土工分析中常用的有限元法等进行了介绍，然后再对计算机内插方法，如拉格朗日内插法、不规则三角网（TIN）内插、标准格网内插法进行了算法研究并编写计算机程序，最后根据实际工程项目，应用软地基处理的监测点进行内插计算，对软地基的沉降和回弹问题进行分析。关键词土工分析、内插方法、软地基沉降1．土工计算机分析的几种常用方法从上世纪20年代起，数值分析方法开始逐步应用于土工分析，特别是计算机技术的发展，使得数值计算在土工分析中的应用更为方便和快捷，计算机软硬件的发展使得很多理论得以在实际工程中应用，有些数学领域的方法借助计算机技术广泛应用于土工分析中。土工分析中常用的数值分析方法主要有有限差分法、有限元法、离散单元法、拉格朗日（Lagrange）元法和流形元法等方法，这些方法都涉及较复杂的数学问题，应用于工程实践则需研究计算机算法，将这些数学模型用计算机能运行的程序进行解算，因此研究计算机算法是实现工程应用的必经之途径，因此本文除对土工分析的数值分析方法作了简单介绍以外，笔者还根据多年的工作实践，编写了计算机程序对软地基的沉降和回弹进行了计算和分析。在求解变分问题中，物理量一般可表示为待定函数及其导数的积分式（如傅里叶函数和傅里叶积分），对其积分区域Ω，可仿照差分法的离散化方法（按一定方式划分网格并取结点），把Ω划分为有限个子区域Ωi，待定函数及其导数在子区域的结点上的值是待求的，这就是有限元法的基本理论，在各个Ωi内解析的插值函数（一般采用多项式插值，多项式系数一般由结点参数来确定，最常用的方法就是Lagrange插值法）来逼近。只要选取适当的插值幂次和结点参数的范围，一般就能满足对插值函数对精度和插值光滑性的要求，如后面还要讲到的2标准格网内插采用的线性内插和五次幂内插法等均被广泛应用。总之，有限元法是以变分原理为基础吸取差分格式的思想而发展起来的一种比较好的数值解法，它把求解无限自由度的待定函数归结为求解有限单元的待定值问题，在数学上讲，具有按一定分布形式的结点及其一定类型的结点参数值的子区域Ωi称为单元。简而言之，有限元法就是将分析域离散成有限个只在结点相连接的子域，即有限元，然后在各单元中采用低阶多项式插值，建立单元特性矩阵，求得各单元的解，再利用能量变分原理集合形成总特性矩阵，最后结合初始条件和边界条件求解。有限元法已广泛应用于地基基础工程、地下工程、堤坝工程、基坑工程及各类土工建筑物的性状分析。除此以外，土工分析还经常用到有限差分法、离散元法、拉格朗日（Lagrange）元法、非连续变形分析法、流形元法等方法。不论采用何种方法，在求解单元插值时，Lagrange插值法都是最常用和最简单的插值方法，也是目前计算机编程插值计算最常用的方法。2．常用的内插法及计算机编程在土工分析中的应用2．1拉格朗日内插法不论采用何种方法在求解单元插值时，拉格朗日内插法都是最常用和最简单的插值方法，也是目前计算机编程插值计算最常用的方法。拉格朗日内插法的基本原理是采用一个n次多项式（一元或二元）去拟合一个函数f(x)或f(x,y)，为计算方便和限制计算误差，一般n不大于5（阶数太大，除增加计算工作量外，有时还会出现不收敛情况和出现很大的计算误差），为保证内插的精度，必要时也采用分段插值，但由于离散点的特性，很容易在分段处其导数不连续，内插曲线不光滑，因此需要在插值结点的小区间内用低幂次插值，以保证内插曲线（或曲面）的光滑。在土工分析中，最常用的方法是在现场采集数量有限的离散样本点，根据样本点求解拉格朗日内插多项式的系数，从而获得内插函数，求解任意点的解，在数学上则表示为：给定n个互异结点xi及对应的函数值yi=f(xi),(i=1,2,···n)，写成多项式为：3)1(121nnxcxccy对大多数问题是求解该多项式的系数ci，将n个离散点的（1）式写成矩阵方程：)2(1111212222112112121nnnnnnnnxxxxxxxxxcccyyy求解的系数Ci的方程式为：)3(11121112122221121121nnnnnnnnyyyxxxxxxxxxccc对于二元五次拉格朗日多项式来说（1）式可写为：)4(),(501kkiiikkiyxcyxz2．2计算机编程内插算法分析对于（4）式，要求解的系数Cki个数为：212)25)(15(n(个)实际工程应用中，往往样本点（m个观测值）多于所要求解的变量Cn(mn），这样就需要对观测值z1、z2、···zm用最小二乘法进行高斯约化求解系数Cn：首先求观测值与理论差的平方和：mjmjkkijijikjkjjjjzyxczyxz1150022)5()()),((再求σ对系数Ckj的偏导数，并令其为0：4mjkkikjjijjkkjkjczyxcc1500)6(0)(2这样，线性代数方程组有21个方程和21个未知系数Ckj，求出Ckj后就确定了二元五次内插的表达式。拉格朗日内插方法一般用于曲线和曲面拟合，在岩土工程中广泛应用于边坡稳定性分析中的滑裂面曲线拟合、基坑回弹分析，地下巷道底臌分析及软地基处理的沉降和回弹分析等。3．不规则三角网（TIN）及标准格网内插及计算机编程在土工分析中的应用3．1TIN的基本理论不规则三角网（TriangulationIrregularNetwork，简写为TIN）内插法是近几年，特别是计算机编程广泛应用于工程实践以来最为普遍的内插方法，它最早由俄国数学家Delaunay于1934年提出的三角形最小内角最大化准则，并证明了在没有四点或四点以上共圆条件下的平面离散点存在一种三角化方式，使得连成的三角网中的三角形满足三角形外接圆不包含第四个及以上的点的条件，该三角形最接近等边三角形，所以也通常称这类三角形为Delaunay三角形。由离散点组成的三角网内插法一般用于空间内插，如表面插值等。其内插方式如图1，内插点P满足共面条件：P1（x1,y1,z1）)7(0cBzAyx·P(x,y,z)P2(x2,y2,z2)P3(x3,y3,z3)图1：三角形内插共面条件A、B、C三个待求系数由三角形三个顶点的已知值求得，求出该三角形的三个系数后，该三角形内任一点（x,y值已知）内插的z值都可以用（7）式求得。图1和（7）式为线性内插，它在计算上是最简单的内插方法，线性内插只5能保证在边界处是连续的，但不能保证光滑，如图2。在边界处不光滑在边界处光滑`法线方向法线方向法线方向法线方向图2：线性内插在边界处的非光滑性图3：高幂次内插在边界处的光滑特性高幂次多项式内插（以二元五次内插为例）的数学公式见（4）式，它类似于线性内插，认为表面模型是连续和光滑的（其二阶导数不为0），从数学上来说，由于五次多项式在内插时考虑到相邻三角形的几何特性，则垂直于每一个三角形的法线是连续变化的，相邻三角形边线上的法线无断裂，内插的表面是光滑的。3．2TIN和标准格网内插在土工分析中的应用笔者应用以上各类内插理论结合计算机编程应用于岩土工程实际，取得了比较好的效果，在此以深圳软地基处理工程为例作一详细介绍。我们承担了深圳西部地区一填海项目的软基处理工程，软基场地面积约1.5万平方米，全部为填海软地基。工程完成后按规范要求进行了4次地基沉降监测，监测期限为12个月，沉降监测时由于场地陆续进行了后续施工，所以布设的监测点有限及部分监测点在施工中受到破坏，最终只有37个监测点有完整的监测数据，监测点的分布比较均匀，其中最大沉降11mm，最大回弹5mm，监测的仪器和精度符合规范要求。笔者将该37个监测点作为样本点，采用TIN和五幂次内插法进行内插计算，以2×2m为单元网格，解算场地内沉降模型。监测点坐标及沉降量见表1。TIN解算软件界面及构成的Delaunay见图4，五次幂曲面内插见图5。内插6计算完成后生成了每个单元网格的三维坐标值数据文件，其中z值代表沉降量，正值代表回弹，负值代表沉降，为便于直观分析，笔者还在ArcGIS下采用二次开发的方法将插值数据导入到ArcScene下用三维图形表示，见图6，为方便观察，其中z值放大400倍表示。表1：软地基沉降观测数据表点号X坐标Y坐标沉降量点号X坐标Y坐标沉降量点号X坐标Y坐标沉降量197034.8214819.260.0021397037.3114850.09-0.0012597086.7514858.360.004297029.1814837.27-0.0011497018.5314913.27-0.0032697081.1614877.070397023.4414855.770.0021597054.4614824.90-0.0022797075.8914896.830.001497016.7114877.110.0011697048.5614846.06-0.0032897069.0314914.460.001597043.5014865.96-0.0051797074.2914830.990.0052997062.8814933.25-0.007697036.9214886.30-0.0051897068.5414853.130.0033097056.9914951.51-0.009797030.9614905.49-0.0081997062.5014872.2803197112.8514843.180897025.5614925.40-0.012097057.2614891.4503297106.6914860.290.002997020.4314942.00-0.0072197050.9614911.15-0.0033397101.1614877.9301097009.9814898.690.0022297045.5214929.76-0.0063497095.7114900.340.0011197004.0314917.340.0012397040.2314947.18-0.0113597089.4714920.180.0021296998.8014934.3902497094.0814836.850.0023697084.1114938.880.00103797079.6014958.02-0.004图4：Delaunay三角网生成图5：五次曲面内插7图6：软地基沉降分析三维图形最后为便于分析，在计算机软件下将沉降数据绘成等值线，见图7。图7：软地基沉降分析等值线参考文献1．龚晓南土工计算机分析中国建筑工业出版社20002．何光渝VisualBasic常用数值算法集科学出版社20023．李胜乐、陆远忠、车时MapInfo地理信息系统二次开发实例电子工业出版社20044．EsriArc/Info帮助文件