第11章程序设计 11.1前言 从本书介绍的内容可以发现,随着近代数值分析技术的进展,边坡稳定分析的极限平衡 分析方法已经发展到十分成熟的程度。但是,随之而出现的问题是,这一分析步骤必须通过 计算机程序才能得以实现。而现状是,尽管计算机的硬件已经获得了飞速的发展,尽管土力 学的理论和实现这些理论的工具和手段均已十分先进,但岩土工程许多应用软件仍然处于起 步状态,远不能满足工程实际的要求,且大部分岩土工程师对边坡稳定分析计算的实际应用 水平还不很高。使用近代科学技术的成果来提高岩土工程设计的效率和水平,是工程师的强 烈愿望,也是保证建筑物的安全和减少投资非常重要手段。当前岩土工程软件的落后状况, 是有多方面因素造成的。 首先,开发岩土工程软件需要投入大量的人力和物力。研究一个新的数值分析方法并编 制相应的科研程序是一回事,将其变为一个能够在工程设计中广泛应用的软件是另一回事 其源代码的数量起码增加十倍。在 Windows的界面上编制一个直观、方便的具有处理图形 对话和纠错功能的软件,不仅是源代码数量上的增加,而且意味着对编程人员素质的新要求, 而这一要求是由所支付的高额报酬体现在程序价格中的 其次,岩土工程软件的市场极为有限,对于这种专业性的软件,单靠市场回报来实现其 自身的发展是不现实的,而这一问题并未得到共识。众所周知,设计人员在图纸上多画一根 预应力锚索,多画一根抗滑桩则是十万甚至百万元的投资。但对于开发这些有关软件的研究 经费(尽管无法与施工费用相比)却难以得到批准。 再次,我们的软件管理还远未达到规范化、市场化的水平。盗版现象普遍存在,因此, 营造一个保护软件的法律和道德约束环境,是软件作为一个产业兴旺发达的最基本的条件 工程应用软件开发的另一个重要问题是确认其可靠性和准确性。工程软件应用的对象是 关系到人民生命财产安全的建筑物,因此,绝非儿戏。对程序进行鉴定和考核自然是一个必 须屡行的手续。但笔者认为对一个包含成千上万行指命的、又需要面临复杂多变的实际情况 的工程应用软件,单靠鉴定和考核是不可能真正发现其全部问题的。这样的程序只有在长期 的、大量的实际应用过程中,才能变得成熟,成为可以信赖的设计工具。 作者希望通过第11章和第12章的内容,推动边坡稳定计算机程序的普及和推广。本章 将介绍编制一个边坡稳定分析程序的核心步骤和源程序。将这几个子程序串联起来,就是一 个计算任意几何形状边坡、包含第2章和第3章介绍的各种分析方法的程序。本章还将讨论 与程序推广、使用相关的问题。在第12章,作者提供STAB程序的一个简化版。使用这 版本,读者可以进行本书介绍的部分稳定分析工作。需要说明的是,本书所列的所有源程序
第11章 程 序 设 计 11. 1 前言 从本书介绍的内容可以发现 随着近代数值分析技术的进展 边坡稳定分析的极限平衡 分析方法已经发展到十分成熟的程度 但是 随之而出现的问题是 这一分析步骤必须通过 计算机程序才能得以实现 而现状是 尽管计算机的硬件已经获得了飞速的发展 尽管土力 学的理论和实现这些理论的工具和手段均已十分先进 但岩土工程许多应用软件仍然处于起 步状态 远不能满足工程实际的要求 且大部分岩土工程师对边坡稳定分析计算的实际应用 水平还不很高 使用近代科学技术的成果来提高岩土工程设计的效率和水平 是工程师的强 烈愿望 也是保证建筑物的安全和减少投资非常重要手段 当前岩土工程软件的落后状况 是有多方面因素造成的 首先 开发岩土工程软件需要投入大量的人力和物力 研究一个新的数值分析方法并编 制相应的科研程序是一回事 将其变为一个能够在工程设计中广泛应用的软件是另一回事 其源代码的数量起码增加十倍 在 Windows 的界面上编制一个直观 方便的具有处理图形 对话和纠错功能的软件 不仅是源代码数量上的增加 而且意味着对编程人员素质的新要求 而这一要求是由所支付的高额报酬体现在程序价格中的 其次 岩土工程软件的市场极为有限 对于这种专业性的软件 单靠市场回报来实现其 自身的发展是不现实的 而这一问题并未得到共识 众所周知 设计人员在图纸上多画一根 预应力锚索 多画一根抗滑桩则是十万甚至百万元的投资 但对于开发这些有关软件的研究 经费(尽管无法与施工费用相比)却难以得到批准 再次 我们的软件管理还远未达到规范化 市场化的水平 盗版现象普遍存在 因此 营造一个保护软件的法律和道德约束环境 是软件作为一个产业兴旺发达的最基本的条件 工程应用软件开发的另一个重要问题是确认其可靠性和准确性 工程软件应用的对象是 关系到人民生命财产安全的建筑物 因此 绝非儿戏 对程序进行鉴定和考核自然是一个必 须屡行的手续 但笔者认为对一个包含成千上万行指命的 又需要面临复杂多变的实际情况 的工程应用软件 单靠鉴定和考核是不可能真正发现其全部问题的 这样的程序只有在长期 的 大量的实际应用过程中 才能变得成熟 成为可以信赖的设计工具 作者希望通过第 11 章和第 12 章的内容 推动边坡稳定计算机程序的普及和推广 本章 将介绍编制一个边坡稳定分析程序的核心步骤和源程序 将这几个子程序串联起来 就是一 个计算任意几何形状边坡 包含第 2 章和第 3 章介绍的各种分析方法的程序 本章还将讨论 与程序推广 使用相关的问题 在第 12 章 作者提供 STAB 程序的一个简化版 使用这一 版本 读者可以进行本书介绍的部分稳定分析工作 需要说明的是 本书所列的所有源程序
336土质边坡稳定分析一原理·方法·程序 只供科研和教学使用,作者不对其承担法律责任,他人也不得将其应用于商业目的 11.2几何图形的识别和分析 11.2.1边坡剖面图形的数字信息 通常,坝体并非均一,而是由若干个物理力学指标不同的土质区域组成。图11.1所示 坝体由6个土质不同的区域组成,控制点17个,边界线22条。图中滑弧上阴影所示的土条 中心线和3条边界线相交,将该土条分成了三段,其高度分别为H、H和H,并分别位于 I区、区和ⅥI区。所有这些,在图解法中,依靠人的直观操作是很容易判断和量测得的 编制程序的关键是如何让计算机代替人判断土条的中心线究竟与哪几个线条相交,各段高度 分别为多少。我们的办法是,让计算机对所有的边界线循环一遍,逐个判断每一条线是否与 土条中心线相交。边界线用所压土层号的数组来表示。输入这些数据,如交上,则记下该交 点的X、Y坐标值及该线段所压的土层号。这样,根据相邻两交点的Y坐标值的差值,即可 算得H、H和H3。根据土层编号,即可查得相应的物理力学指标,随即可算得土条重量、 地震力、孔隙水压力等。这一对几何图形进行处理的方法,充分利用了计算机可进行大量重 复的判断,忌讳过多的逻辑判断的特点,为编制一个效率高、鲁棒性( Robust)强的边坡稳 定分析程序创造了条件。 10 图11.1由6个土质不闻的区域组成坝体 11.2.2构造圆弧滑裂面的步骤 设某一坝坡坡面为具有n个折点的折线,如图11.2所示(其中起点1和尾点n在极远处)。 建立坐标系统,X轴水平,指向滑坡方向为正,Y坐标轴垂直向下为正。某一半径为R0的滑 弧圆心O至坝坡各折点距离分别为R,R2,…,Rn。该圆弧自左至右在连接、计1两点的线 段进入坝体(上交点),在连接八、+1线段逸出坝体(下交点)。 不难证明,某一半径为R0的圆弧与线段A1A2有一个且只有一个交点的必要充分条件是: R2≤R0≤R1,其中R1、R2为圆心至线段端点的距离。如果R同时大于R1、R2,则圆弧与A1A2 没有交点;如果R0同时小于R1、R2,则圆弧与A1A2不是没有交点,就是有两个交点见图
336 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 只供科研和教学使用 作者不对其承担法律责任 他人也不得将其应用于商业目的 11. 2 几何图形的识别和分析 11. 2. 1 边坡剖面图形的数字信息 通常 坝体并非均一 而是由若干个物理力学指标不同的土质区域组成 图 11.1 所示 坝体由 6 个土质不同的区域组成 控制点 17 个 边界线 22 条 图中滑弧上阴影所示的土条 中心线和 3 条边界线相交 将该土条分成了三段 其高度分别为 H1 H2和 H3 并分别位于 I 区 V 区和 VI 区 所有这些 在图解法中 依靠人的直观操作是很容易判断和量测得的 编制程序的关键是如何让计算机代替人判断土条的中心线究竟与哪几个线条相交 各段高度 分别为多少 我们的办法是 让计算机对所有的边界线循环一遍 逐个判断每一条线是否与 土条中心线相交 边界线用所压土层号的数组来表示 输入这些数据 如交上 则记下该交 点的 X Y 坐标值及该线段所压的土层号 这样 根据相邻两交点的 Y 坐标值的差值 即可 算得 H1 H2和 H3 根据土层编号 即可查得相应的物理力学指标 随即可算得土条重量 地震力 孔隙水压力等 这一对几何图形进行处理的方法 充分利用了计算机可进行大量重 复的判断 忌讳过多的逻辑判断的特点 为编制一个效率高 鲁棒性 (Robust) 强的边坡稳 定分析程序创造了条件 图 11. 1 由 6 个土质不同的区域组成坝体 11. 2. 2 构造圆弧滑裂面的步骤 设某一坝坡坡面为具有n个折点的折线 如图11.2所示(其中起点1和尾点n在极远处) 建立坐标系统 X 轴水平 指向滑坡方向为正 Y 坐标轴垂直向下为正 某一半径为 R0的滑 弧圆心 O′ 至坝坡各折点距离分别为 R1, R2, …, Rn 该圆弧自左至右在连接 i i+1 两点的线 段进入坝体 上交点 在连接 j j+1 线段逸出坝体 下交点 不难证明 某一半径为R0的圆弧与线段A1A2有一个且只有一个交点的必要充分条件是 R2 R0 R1 其中 R1 R2为圆心至线段端点的距离 如果 R0同时大于 R1 R2 则圆弧与 A1A2 没有交点 如果 R0 同时小于 R1 R2 则圆弧与 A1A2 不是没有交点 就是有两个交点,见图
第11章程序设计337 1.3。 则根据上述原理可知,i、j应满足下面条件: 图112國弧滑裂面与边坡相 图11.3具有n个折点的边坡线与圆弧滑裂面相交的关系 交的关系 R1+1≤R,+sT+(+s)Ro-(, -T)'I (11.6) T (11.7) 此联立方程组有两个根,计算下交点时取“+号,计算上交点时取‘-’号 11.2.3构造任意形状滑裂面的步骤 如果要将若干点用光滑的曲线相连,可以采用样条函数予以实现,在STAB程序中,采 用了第二种边界条件的三次样条函数。 第二种边界条件三次样条函数曲线两端点的二阶导数值为已知值,即y12=y(x) yn2=y(xn)。在构造滑裂面时取y2=y(x1)=0,yn2=y(xn)=0,由此确定的样条称为
第 11 章 程序设计 337 11.3 则根据上述原理可知 i j 应满足下面条件 图 11. 3 具有 n 个折点的边坡线与圆弧滑裂面相交的关系 图 11. 2 圆弧滑裂面与边坡相 交的关系 (11.1) Ri+1 ≤ R0 < Ri (11.2) R j ≤ R0 < R j+1 找到了上 下交点所在的线段号 i j 之后 就可以根据该线段两端点的坐标值和圆弧的 圆心坐标值 x0 y0和半径值 R 算得其交点的坐标值 对于点 m (m = i, j) 计算公式为 (11.3) 2 2 0 2 0 (x − x ) + ( y − y ) = R (11.4) ( y − ym )/(x − xm ) = s (11.5) ( )/( ) s = ym+1 − ym xm+1 − xm 联立以上方程组 求解得 [ (1 ) ( ) ] (1 ) 1 2 0 0 2 2 x0 sT s R sx T s x + ± + − − + = (11.6) (11.7) 0 T = sxm − ym + y 此联立方程组有两个根 计算下交点时取 + 号 计算上交点时取 − 号 11. 2. 3 构造任意形状滑裂面的步骤 如果要将若干点用光滑的曲线相连 可以采用样条函数予以实现 在 STAB 程序中 采 用了第二种边界条件的三次样条函数 第二种边界条件三次样条函数曲线两端点的二阶导数值为已知值 即 ( ) 12 1 y = y′′ x ( ) n2 n y = y′′ x 在构造滑裂面时取 y12 = y′′(x1 ) = 0 yn2 = y′′(xn ) = 0 由此确定的样条称为
土质边坡定分析一原理·方法程序 自然样条。根据此边界条件,计算n个节点上的一阶导数值m,j=12,…,n。相应的计算 公式如下 a1=-0.5 y(x2)-y(x1) 2(x2-x1) (11.9) (11.10) B=3(1-a,y-y-)/h-1+a(y1+1-y)/h b=B-(-a,)12+(-,-1J=2.3…,n-1 mn=[3(vn-yn-)/hn1-bn-14(2+an1) (11.15) (11.16) 于是对于x≤x<x/H,利用数值导数m,(j=12,…,n),计算插值节点处函数值的公式 如下 -x)2-2(xm-x)2p+2(x-x)2-2(x-x)Py +b(x1+1-x) h,(x-x1)2 (x-x 如果插值区间(x1,xH1)被等距划分为段,则位于第s(0<s<段的被插值节点 x=x1+(s/)h处的函数值的计算公式如下 hy,(K3+hs, (E (11.18) 其中 hg=(x-x)/h,=s/ (11.19) K1=2y(x)-y(x1+1)+(m2+m+)h (11.20) 5y(x+)-y(x)-(2m+m/+1) 3-m K4=y(i)
338 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 自然样条 根据此边界条件 计算 n 个节点上的一阶导数值m j n j , = 1,2,L, 相应的计算 公式如下 0.5 (11.8) a1 = − 2( ) ( ) ( ) 3 2 1 2 1 1 x x y x y x b − − = × (11.9) (11.10) 1,2, , 1 h j = x j+1 − x j j = L n − ( ) α j = h j−1 hj−1 + h j (11.11) 3[(1 )( ) ( ) ] j j j j 1 j 1 j j 1 j j = − y − y h + y − y h β α − − α + (11.12) [2 (1 ) ] j = − j + − j j−1 a α α a (11.13) = [ − (1− ) 1] [2 + (1− ) 1] = 2,3, , −1 − − b b a j j n j j β j α j j α L (11.14) [3( ) ]/(2 ) n = n − n−1 hn−1 − bn−1 + an−1 m y y (11.15) (11.16) 1, 2, ,1 m j = a jm j+1 + b j j = n − n − L 于是对于 , 利用数值导数 ( j ) j ≤ < j+1 x x x m j = 1,2,L,n 计算插值节点处函数值的公式 如下 1 3 3 2 2 3 1 3 2 1 2 1 3 3 2 2 3 1 3 2 1 2 ( ) 2 ( ) 1 ( ) 1 ( ) 1 ( ) 2 ( ) 3 ( ) 2 ( ) 3 ( ) + + + + + + − − − − + − − − + − − − = − − − j j j j j j j f j j j f j j j j j j j j j j x x m h x x h x x m h h x x h h x x y h x x h x x y h x x h y x (11.17) 如果插值区间 (x j , x j+1 ) 被等距划分为 l 段 则位于第 s (0<s<l) 段的被插值节点 j j x = x + (s l)h 处的函数值的计算公式如下 (11.18) ( ) [ ( )] κ4 sj κ3 sj κ2 sjκ1 y x = + h + h + h 其中 h x x h s l sj = ( − j) j = (11.19) (11.20) j j j j j 2[ y(x ) y(x )] (m m )h κ1 = − +1 + + +1 (11.21) j j j j j 3[ y(x ) y(x )] (2m m )h κ2 = +1 − − + +1 (11.22) m j h j = ⋅ κ 3 ( ) (11.23) 4 j κ = y x
第11章程序设计339 11.3边坡稳定分析的源程序 11.3.1概述 本章介绍边坡稳定分析程序STAB的核心步骤和程序。程序 CIRCLE. FOR具体实施 1122节的计算步骤,对给定的边坡剖面信息,构造一个圆弧滑裂面。程序 SPLINE FOR具 体实施112、3节的计算步骤,对给定的滑面控制点信息构造一个任意形状滑裂面。程序 PROFILE FOR具体实施11.2.1节的计算步骤,计算土条重量和物理力学参数。程序SFOR 计算安全系数,包含了第2章和第3章介绍的各种分析方法的源程序,即极限平衡法的通用 条分法、瑞典法、毕肖普法、陆军工程师团法和简化法。将上述几个子程序串联起来,就是 一个计算复杂几何形状边坡稳定的安全系数的程序 在介绍这些程序时,将配以图114所示的一个例子。这一例子是将在114节中介绍的 澳大利亚 ACADS程序考核的EX1(c)。该例已在第2章和第4章中多次使用。读者在移植这 些程序时,可通过这些考核题进行考核 图11.4滑裂面和条分 11.3.2构造圆弧滑裂面的程序 CIRCLE. FOR 1.说明 本子程序对给定的边坡剖面信息,构造一个圆弧滑裂面(图114)。通过子程序SEAR找 到圆弧滑裂面与边坡的上、下交点所在线段。子程序DIVI计算上、下交点的坐标,并划分 土条。 2.源程序 REAL YTENSION, DS, CX, CY, XN(99), YN(99), X(80), Y(80), ALF(80) INTEGER INI, LSL, N, LNUM(80), IC(80,3), TWRS, L, NN, NLOW, NUPP OPEN(S, FILE=, STATUS=UNKNOWN) OPEN(6, FILE=,STATUS=UNKNOWN ILOW=0 NUPP=0 READGS, )IWRS, YTENSION, DS, CX,CY IWR5=0,1分别相应顶部有、无拉力缝; YTENSION=拉力缝底y坐标 DS. CX.CY=滑弧深度,圆心坐标 READ(5,*)N土条总数 READ(,*)INl!线段总数 DO 3021=L.INI READ5,*)(C(D)J=1,3)线段信息 302 CONTINUE
第 11 章 程序设计 339 11. 3 边坡稳定分析的源程序 11. 3. 1 概述 本章介绍边坡稳定分析程序 STAB 的核心步骤和程序 程序 CIRCLE.FOR 具体实施 11.2.2 节的计算步骤 对给定的边坡剖面信息 构造一个圆弧滑裂面 程序 SPLINE.FOR 具 体实施 11.2.3 节的计算步骤 对给定的滑面控制点信息构造一个任意形状滑裂面 程序 PROFILE.FOR 具体实施 11.2.1 节的计算步骤 计算土条重量和物理力学参数 程序 S.FOR 计算安全系数 包含了第 2 章和第 3 章介绍的各种分析方法的源程序 即极限平衡法的通用 条分法 瑞典法 毕肖普法 陆军工程师团法和简化法 将上述几个子程序串联起来 就是 一个计算复杂几何形状边坡稳定的安全系数的程序 在介绍这些程序时 将配以图 11.4 所示的一个例子 这一例子是将在 11.4 节中介绍的 澳大利亚 ACADS 程序考核的 EX1(c) 该例已在第 2 章和第 4 章中多次使用 读者在移植这 些程序时 可通过这些考核题进行考核 图 11. 4 滑裂面和条分 11. 3. 2 构造圆弧滑裂面的程序 CIRCLE.FOR 1. 说明 本子程序对给定的边坡剖面信息 构造一个圆弧滑裂面(图 11.4) 通过子程序 SEAR 找 到圆弧滑裂面与边坡的上 下交点所在线段 子程序 DIVI 计算上 下交点的坐标 并划分 土条 2. 源程序 REAL YTENSION, DS, CX, CY, XN(99), YN(99), X(80), Y(80), ALF(80) INTEGER IN1, LSL, N, LNUM(80), IC(80,3), IWR5, L, NN, NLOW, NUPP OPEN(5,FILE='',STATUS='UNKNOWN') OPEN(6,FILE='',STATUS='UNKNOWN') NLOW=0 NUPP=0 READ(5,*)IWR5,YTENSION,DS,CX,CY C IWR5=0,1 分别相应顶部有 无拉力缝 YTENSION=拉力缝底 Y 坐标 C DS,CX,CY=滑弧深度 圆心坐标 READ(5,*)N!土条总数 READ(5,*)IN1!线段总数 DO 302 I=1,IN1 READ(5,*)(IC(I,J),J=1,3)!线段信息 302 CONTINUE
340土质边坡德定分析一原理·方法程序 READ(5,*)LSL外边坡线总数 READ(S,*)LNUM()l=1,LSL)外边坡线编号 READ(5,NN!点总数 READ,*)L,XN(I,YN(D)1=1,NN)点坐标 703 FORMAT(T25,***) WRITE(6, 704) ORMAT(IOX, THEABSCISSAVALUESOFTHEUPPERANDLOWER') CALL SEARCLNUM,C,DS,CXCY,XN,YN,LSL, NLOW NUPP)!寻找上、下交点所在线段 CALL DIVI(NIC,IWRS,DS,CXCY,XN,YN!确定上、下交点坐标,条分 S, YTENSION, NLOW, NUPP X, Y, ALI 711 FORMAT(TS, NO, T17, X,T32,Y, T47,ALF) DO309I=1 WRITE(6,710)L, X),Y(),ALFO 710 FORMAT(1X,15,3F15.6) 309 CONTINUE 12 CALLCLOSEFL END SUBROUTINE CLOSEFL CLOSE(5) CLOSE(6) RETURN SUBROUTINE DIVI(N, IC, TWRS, DS, CX, CY,XN,YN S, YTENSION, NLOW,NUPP,X,Y,ALF) INTEGERNS, N, K, IC(80,3), KI, KTEMP, TWRS,J,J1,JTEMPNLOW, NUPP REAL RO, DS, CY,XN(80), YN(80),SK, T, XL, CX, XUl, YTENSION REAL YUI,SJ, XU, ROO, YU, X2(80), Y2(80), X(80), Y(80), ALF(80) S(E1,F1E2,F2)=(F2-F1)(E2-E1) NS=N-1 RO=DS-CY KI=IC(NLOW, 2) IF(XN(K).GTXN(KDGOTO190 KTEMP=K KI=KTEMP 190 IF(ABS(XN(K)-XN(KI).LT0.00005 ) GOTO23 SK=S(XN(K),YN(K), XN(KD,YN(KD) T=SK*XN(K)-YN(K)+CY XL=(CX+SK*T+SQRT(1+SK*2)RO*2-(SK*CX-T)*2))(1+SK**2) !按式(16)计算下交点的横坐标 GOTO28 23 XLFXNK 28 CONTINUE IF(WRS NE.O)THEN XUI=CX-SQRT(RORO-(YTENSION-CY*2) YUI=CY+SQRT(RO**2-ABS(XUl-CX)**2) J=IC(NUPP, D) IF(XNO).GT.XNJDGOTO191 JTEMP- JI=JTEMP
340 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 READ(5,*)LSL!外边坡线总数 READ(5,*)(LNUM(I),I=1,LSL)!外边坡线编号 READ(5,*)NN!点总数 READ(5,*)(L,XN(I),YN(I),I=1,NN)!点坐标 WRITE(6,703) 703 FORMAT(T25,'******************'//) WRITE(6,704) 704 FORMAT(10X,'THEABSCISSAVALUESOFTHEUPPERANDLOWER') CALL SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP)!寻找上 下交点所在线段 CALL DIVI(N,IC,IWR5,DS,CX,CY,XN,YN 确定上 下交点坐标 条分 $,YTENSION,NLOW,NUPP,X,Y,ALF) WRITE(6,711) 711 FORMAT(T5,'NO.',T17,'X',T32,'Y',T47,'ALF') DO 309 I=1,N WRITE(6,710)I,X(I),Y(I),ALF(I) 710 FORMAT(1X,I5,3F15.6) 309 CONTINUE 12 CALLCLOSEFL END SUBROUTINE CLOSEFL CLOSE(5) CLOSE(6) RETURN END SUBROUTINE DIVI(N,IC,IWR5,DS,CX,CY,XN,YN $,YTENSION, NLOW,NUPP,X,Y,ALF) INTEGERN S,N,K,IC(80,3),K1,KTEMP,IWR5,J,J1,JTEMP,NLOW,NUPP REAL RO,DS,CY,XN(80),YN(80),SK,T,XL,CX,XU1,YTENSION REAL YU1,SJ,XU,ROO,YU,X2(80),Y2(80),X(80),Y(80),ALF(80) S(E1,F1,E2,F2)=(F2-F1)/(E2-E1) NS=N-1 RO=DS-CY K=IC(NLOW,1) K1=IC(NLOW,2) IF(XN(K).GT.XN(K1))GOTO190 KTEMP=K K=K1 K1=KTEMP 190 IF(ABS(XN(K)-XN(K1)).LT.0.00005)GOTO23 SK=S(XN(K),YN(K),XN(K1),YN(K1)) T=SK*XN(K)-YN(K)+CY XL=(CX+SK*T+SQRT((1+SK**2)*RO**2-(SK*CX-T)**2))/(1+SK**2) !按式(11.6)计算下交点的横坐标 GOTO28 23 XL=XN(K) 28 CONTINUE IF(IWR5.NE.0)THEN XU1=CX-SQRT(RO*RO-(YTENSION-CY)**2) YU1=CY+SQRT(RO**2-ABS(XU1-CX)**2) ENDIF J=IC(NUPP,1) J1=IC(NUPP,2) IF(XN(J).GT.XN(J1))GOTO191 JTEMP=J J=J1 J1=JTEMP
第11章程序设计341 191 F(ABS(XNO)-XNOI)LT0.00005)GOTO82 J=S(XN),YN),XNI),YNJD) T=SJ"XNOFYNOHCY XU=Cx+SJ*T-SQRI(1+SJ**2)RO**2-(SJCX-)*2))(1+SJ**2) XU上交点的横坐标 ROO=RO**2-ABS(XU-CX)**2 IF( ROO.LT0 RETURN无交点,需修改RO,重新计算 YU=CY+SQRT(RO*#2-ABS(XU-CX) 2) IF( CYGTYN( J)AND CYGT.YN(J) RETURN凸弧否定,重算 IF(WRSNE.0 AND YU. LTYUI)XU=XUI XU=XNO) DO541=1,NS划分土条求各点坐标 X2(=(-)+(XLXU(NS-1)+XU!土条宽度B=(XL-XU(NS-1) Y2(FCY+SQRT(RO*2-ABS(X2()-CX)**2) 41 CONTINUE IF(ABS(Y2(NS-Y2(1)). LT0.OlTHEN WRITE(6, THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE RETURN ENDIF N=NS+1 Y2(NFY2(NS 计算土条中点坐标 DO544=2.NS X(=(X2(1)+X2(1-1)2 Y(=Y2(+Y2(-1)y ALF(I=ATAN(Y2(1)Y2(-1))(X2(1}X2(-1)) 4 CONTINUE ALF(IFALF(2) ALF(NALF(NS) X(NX2(N) X(1)=x2(1) Y(1)=Y2(1) RETURN END SUBROUTINE SEAR(LNUM,IC, DS, CX, CY,XN,YN, LSL, NLOW, NUPP) REAL R, DS, CY, R2, CX, XN(80),YN(80),Rl, D, XXL, SLO INTEGER I, I1, K, LNUM(80),J1, J2, IC(80,3), IN, NLOW,NUPP LSL I=1NLOW=下交点所在直线段的线段号 I1=0NUPP=上交点所在直线段的线段号 12K=LNUMD JI=IC(K, 1) J2=IC(K, 2) R= DS-CY!R=圆弧半 R2=SQRT(ABS(CX-XN(2)**2+ABS(CY-YNO2)* *2) RI=SQRT(ABS(CX-XNO1)**2+ABS(CY-YNJD))**2) R,R2为圆心至该线段两端点的距离 IF(REQ R2AND II.NEO)GOTOll IF(RGTR1 ANDR GTR2)GOTO11!无交点 IF(RLTR1 ANDRLTR2GOTO19有两个或零个交点 GOT020 19 CALL DD(XNJI),YND), XN(2), YN(J2), CX, CY, D, IN) XXL=XN(2)-XNO1
第 11 章 程序设计 341 191 F(ABS(XN(J)-XN(J1)).LT.0.00005)GOTO82 SJ=S(XN(J),YN(J),XN(J1),YN(J1)) T=SJ*XN(J)-YN(J)+CY XU=(CX+SJ*T-SQRT((1+SJ**2)*RO**2-(SJ*CX-T)**2))/(1+SJ**2) !XU 上交点的横坐标 ROO=RO**2-ABS(XU-CX)**2 IF(ROO.LT.0)RETURN!无交点,需修改 RO,重新计算 YU=CY+SQRT(RO**2-ABS(XU-CX)**2) IF(CY.GT.YN(J).AND.CY.GT.YN(J1))RETURN!凸弧,否定,重算 IF(IWR5.NE.0.AND.YU.LT.YU1)XU=XU1 GOTO20 82 XU=XN(J) 20 X2(1)=XU DO541I=1,NS!划分土条,求各点坐标 X2(I)=(I-1)*(XL-XU)/(NS-1)+XU!土条宽度 B=(XL-XU)/(NS-1) Y2(I)=CY+SQRT(RO**2-ABS(X2(I)-CX)**2) 541 CONTINUE IF(ABS(Y2(NS)-Y2(1)).LT.0.01)THEN WRITE(6,*)'THE SLIP SURFACE INTERCEPTS A HORIZONTAL SURFACE' RETURN ENDIF N=NS+1 X2(N)=X2(NS) Y2(N)=Y2(NS) C 计算土条中点坐标 DO544I=2,NS X(I)=(X2(I)+X2(I-1))/2. Y(I)=(Y2(I)+Y2(I-1))/2. ALF(I)=ATAN((Y2(I)-Y2(I-1))/(X2(I)-X2(I-1))) 544 CONTINUE ALF(1)=ALF(2) ALF(N)=ALF(NS) X(N)=X2(N) Y(N)=Y2(N) X(1)=X2(1) Y(1)=Y2(1) RETURN END SUBROUTINE SEAR(LNUM,IC,DS,CX,CY,XN,YN,LSL,NLOW,NUPP) REAL R,DS,CY,R2,CX,XN(80),YN(80),R1,D,XXL,SLO INTEGER I,I1,K,LNUM(80),J1,J2,IC(80,3),IN,NLOW,NUPP,LSL I=1 !NLOW=下交点所在直线段的线段号 I1=0 !NUPP=上交点所在直线段的线段号 12K=LNUM(I) J1=IC(K,1) J2=IC(K,2) R=DS-CY!R=圆弧半径 R2=SQRT(ABS(CX-XN(J2))**2+ABS(CY-YN(J2))**2) R1=SQRT(ABS(CX-XN(J1))**2+ABS(CY-YN(J1))**2) !R1,R2 为圆心至该线段两端点的距离 IF(R.EQ.R2.AND.I1.NE.0)GOTO11 IF(R.GT.R1.AND.R.GT.R2)GOTO11 !无交点 IF(R.LT.R1.AND.R.LT.R2)GOTO19 !有两个或零个交点 GOTO20 19 CALL DD(XN(J1),YN(J1),XN(J2),YN(J2),CX,CY,D,IN) XXL=XN(J2)-XN(J1)
342土质边坡德定分析一原理·方法程序 IF(ABS(XXL).LT0.0001THEN SLO=L0E-6 ELSE SLO=(YN(J2YN(J1)XXL!斜率 ENDIF IF(NEQ 1 OR SLO. LEO)GOTOll C第一条外边坡线为水平同时圆弧与其有两个交点不接受继续搜索 IF( R. LT D ORREQ. D)GOTO1无交点 NLOW=LNUM NUPPENLOW RETURN 0 IF(I NE.OGOTO14 NLOW=LNUM 11=11+I 11 CONTINUE IF(EQ 1 AND R.LTRIAND R GT R2)RETURN l=+1 IF(GT.LSLGOTO14 GOTO12 14 NUPP=LNUM 15 CONTINUE RETURN SUBROUTINE DD(X1, Y1, X2, Y2, XC,YC, D, IN) REAL XI, Y1, X2, Y2, XC,YC, D REAL XXL, AK, AKF X,YXT,YT NTEGER IN XXLEX2-XI IF(ABS(XXL). LT0.001)XXL=1.0E-6 IF(ABS((Y2-Y1)XXL). LT0.OOIRETURN AK=(Y2-YI/XXL AKF=l/AK X=(AK*Xl-Y1+AKF*XC+YC)(AK+AKF) Y=AK*XAK*X1+Y1!从圆心作垂直与线段的线获与线段的交点 X1OX-X2 YT=(YY1)*(Y-Y2) IF(XT.LT0OR YT. LT.OGOTO12 IN=11交点不在线段内 RETURN 12D=SQRT(X-XC)**2+(YYO)*2)!线段与圆心的垂直距离 RETURN END 3.例题 数据文件 CIRCLE DAT 00. 24.235-35.164-42.046 21 2 3
342 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 IF(ABS(XXL).LT.0.0001)THEN SLO=1.0E-6 ELSE SLO=(YN(J2)-YN(J1))/XXL !斜率 ENDIF IF(IN.EQ.1.OR.SLO.LE.0)GOTO11 C 第一条外边坡线为水平,同时圆弧与其有两个交点,不接受,继续搜索 IF(R.LT.D.OR.R.EQ.D)GOTO11 !无交点 NLOW=LNUM(I) NUPP=NLOW RETURN 20 IF(I1.NE.0)GOTO14 NLOW=LNUM(I) I1=I1+1 11 CONTINUE IF(I.EQ.1.AND.R.LT.R1.AND.R.GT.R2)RETURN I=I+1 IF(I.GT.LSL)GOTO14 GOTO12 14 NUPP=LNUM(I) 15 CONTINUE RETURN END SUBROUTINE DD(X1,Y1,X2,Y2,XC,YC,D,IN) REAL X1,Y1,X2,Y2,XC,YC,D REAL XXL,AK,AKF,X,Y,XT,YT INTEGER IN D=9999. IN=0 XXL=X2-X1 IF(ABS(XXL).LT.0.001)XXL=1.0E-6 IF(ABS((Y2-Y1)/XXL).LT.0.001)RETURN AK=(Y2-Y1)/XXL AKF=1./AK X=(AK*X1-Y1+AKF*XC+YC)/(AK+AKF) Y=AK*X-AK*X1+Y1 !从圆心作垂直与线段的线,获与线段的交点 XT=(X-X1)*(X-X2) YT=(Y-Y1)*(Y-Y2) IF(XT.LT.0.OR.YT.LT.0)GOTO12 IN=1 !交点不在线段内 RETURN 12 D=SQRT((X-XC)**2+(Y-YC)**2) !线段与圆心的垂直距离 RETURN END 3. 例题 数据文件 CIRCLE.DAT 0 0. -24.235 -35.164 -42.046 21 9 1 2 3 2 3 1 3 4 1 2 5 3 5 6 2 6 7 2
第11章程序设计343 000-25000 70.000-35.000 -40.000-27.000 54.000-31.000 70.000-31.000 52.000-24.000 10-70.000-24.000 计算成果 CIRCLEC: THE ABSCISSA VALUES OF THE UPPER AND LOWEI NO 51.522050 1.094647 50.955680 3-49822940-31.978910 -48.690200 30491210 863838 -47.557460 29278430 770531 46.424730 28265840 685123 7-45.29199027411070 605365 41.893780 25.566700 387709 11-40.761050 -25.147810 319849 1239628310-24.813490 13-38495570 24.558860 188260 14-37.362830-24.380470 123833 3509736024244130 003743 33.964630-24.284500 067425 1832831890 -131382 4.584810 -.19588 30.566420 4848610 261233 30.000050 5.000020 261233 11.3.3构造任意形状滑裂面的程序 SPLINE FOR 1.说明 本程序(刘德贵,1980)用来计算已知第二种边界条件(曲线两端点处的二阶导数值) 三次样条函数的一阶导数值,并利用样条函数对一元函数进行插值,进而构筑包含直线段的 复合滑裂面。对于滑裂面上控制点区间为直线段时,则直接线性插值,而对于非直线段区间 则运用上述介绍的方法进行三次样条插值。滑裂面被NS1个点分为NS1-1段。从上交点向 下交点编号为1,2,…NS1-1。程序用LS代表此NS1-1段中为直线段的线段总数。在数组 LSA(I)(=1,2,…LS)中存入这些直线段的编号 LS,滑裂面上NS1-1条线段中为直线段的总数 SA(LS),滑裂面上线段为直线段的编号;
第 11 章 程序设计 343 7 8 2 5 9 3 9 10 3 3 1 2 3 10 1 .000 -25.000 2 -30.000 -25.000 3 -50.000 -35.000 4 -70.000 -35.000 5 -40.000 -27.000 6 -50.000 -29.000 7 -54.000 -31.000 8 -70.000 -31.000 9 -52.000 -24.000 10 -70.000 -24.000 计算成果 CIRCLE.C THE ABSCISSA VALUES OF THE UPPER AND LOWER NO. X Y ALF 1 -51.522050 -35.000000 1.094647 2 -50.955680 -33.901810 1.094647 3 -49.822940 -31.978910 .969015 4 -48.690200 -30.491210 .863838 5 -47.557460 -29.278430 .770531 6 -46.424730 -28.265840 .685123 7 -45.291990 -27.411070 .605365 8 -44.159260 -26.687410 .529825 9 -43.026520 -26.076840 .457513 10 -41.893780 -25.566700 .387709 11 -40.761050 -25.147810 .319849 12 -39.628310 -24.813490 .253490 13 -38.495570 -24.558860 .188260 14 -37.362830 -24.380470 .123833 15 -36.230100 -24.275990 .059923 16 -35.097360 -24.244130 -.003743 17 -33.964630 -24.284500 -.067425 18 -32.831890 -24.397580 -.131382 19 -31.699150 -24.584810 -.195887 20 -30.566420 -24.848610 -.261233 21 -30.000050 -25.000020 -.261233 11. 3. 3 构造任意形状滑裂面的程序SPLINE.FOR 1. 说明 本程序 刘德贵 1980 用来计算已知第二种边界条件(曲线两端点处的二阶导数值) 三次样条函数的一阶导数值 并利用样条函数对一元函数进行插值 进而构筑包含直线段的 复合滑裂面 对于滑裂面上控制点区间为直线段时 则直接线性插值 而对于非直线段区间 则运用上述介绍的方法进行三次样条插值 滑裂面被 NS1 个点分为 NS1−1 段 从上交点向 下交点编号为 1, 2, , NS1−1 程序用 LS 代表此 NS1−1 段中为直线段的线段总数 在数组 LSA(I)(I=1, 2, , LS)中存入这些直线段的编号 LS, 滑裂面上 NS1−1 条线段中为直线段的总数 LSA(LS) 滑裂面上线段为直线段的编号
344土质边坡稳定分析一原理·方法·程序 AI(14),滑裂面上控制点处四个与样条插值有关的数,K,K2,K3,K4; A(),滑裂面上控制点处的一阶导数值m B(D),一维数组,工作单元 2.源程序 SUBROUTINE SPI(LS, LSA, X2, Y2) DIMENSION KO2(50),X1(50,Y1(50),AI(50,4)LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1, KQ2, X1, Y1/A13/DD, IWRl, IWR2, IWR3 COMMON/WALL/TWALL, G WALL, HMW,EWALL, ETA INTEGER4 OPTION(6),OP1(6) COMMON/OPP/OPTION. OPI OPEN(S, FILE STATUS=UNKNOWN OPEN(6, FILE STATUS=UNKNOWN READ(S NSI DO=LNSI READO, )KQ2), Xl,YI) ENDDO 如果滑裂面上控制点总数NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的xy值 LSAFO ENDDO IF(LS.NEO)READ(, (LSA I=l, LS) sA(S+IFNS C 滑裂面上第NS个控制点对应的土条侧面边界线的编号 Xl(①),Yl(),滑裂面上第I个控制点的x,y值 N1=0 DO 100IL=LLS+ LI=NI+1 N2=NI-1 IF(LSEQ0GOTO83!滑裂面上没有直线段 IF(NIGE NS1)GOTO 40 AI(Nl,1=0.!计算直线段左端点处K1,K2,K,K4 AI(Nl,3)=Y1(N1+1}Y1(N1) AI(N1, 4FYI(ND) 40IF(N1L1EQ0GOTO100!滑裂面上没有直线段 83B(Ll=1.5°(Y1(L+1Y1(L1)/(X1(L1+1)-Xl(L1) A(L1)=0.5!插值区间左端点处的一阶导数值 CN=3.0*(Y1(Nl)Y1(N2)(X1(Nl}X1(N2) IF(N1-L1EQ.1)GOTO55!插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2 N2 HI=XIO)(-1) H=XIO+1)-XI) AF=H1/(H+H1)!计算α值 BT=3.*(1.-AF)°Y1(J-Y1(J-1)MHI+AF(Y1(J+1}Y1)/H) 计算β值 (JAF(2.+(1.-AF)*A(J-)!计算插值区间内插值节点处的一阶导数值m 50B(J=(BT-(1.-AF)*B(J-1)/(2+(1.-AF)A(J-1) 5A(NI=(CN-B(N2)/(2.+A(N2) A(N2=A(N2)*A(N2+1)+B(N2)!计算插值区间右端点处的一阶导数值mn
344 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 AI(I,4) 滑裂面上控制点处四个与样条插值有关的数 κ1, κ2, κ3, κ4 A(I) 滑裂面上控制点处的一阶导数值m j B(I) 一维数组 工作单元 2. 源程序 C SUBROUTINE SP1(LS,LSA,X2,Y2) DIMENSION KQ2(50),X1(50),Y1(50),AI(50,4),LSA(60),A(50),B(50) DIMENSION X2(80),Y2(80) COMMON /A23/NS1,KQ2,X1,Y1/A13/IDD,IWR1,IWR2,IWR3 COMMON/WALL/IWALL,GWALL,HMW,EWALL,ETA INTEGER*4 OPTION(6),OP1(6) COMMON/OPP/OPTION,OP1 OPEN(5,FILE=' ',STATUS='UNKNOWN') OPEN(6,FILE=' ',STATUS='UNKNOWN') READ(5,*)NS1 DO I=1,NS1 READ(5,*)KQ2(I),X1(I),Y1(I) ENDDO C 如果滑裂面上控制点总数 NS1=2,则按直线段处理并计算滑裂面与土条侧面边界线交点的 x,y 值; READ(5,*)LS DO I=1,NS1-1 LSA(I)=0 ENDDO IF(LS.NE.0)READ(5,*)(LSA(I),I=1,LS) LSA(LS+1)=NS1 C 滑裂面上第 NS1 个控制点对应的土条侧面边界线的编号 C X1(I) Y1(I) 滑裂面上第 I 个控制点的 x y 值 N1=0 DO 100 IL=1,LS+1 L1=N1+1 N1=LSA(IL) N2=N1-1 IF(LS.EQ.0)GOTO 83 ! 滑裂面上没有直线段 IF(N1.GE.NS1) GOTO 40 AI(N1,1)=0. ! 计算直线段左端点处κ1,κ2,κ3,κ4 AI(N1,2)=0. AI(N1,3)=Y1(N1+1)-Y1(N1) AI(N1,4)=Y1(N1) 40 IF(N1-L1.EQ.0)GOTO 100 ! 滑裂面上没有直线段 83 B(L1)=1.5*(Y1(L1+1)-Y1(L1))/(X1(L1+1)-X1(L1)) A(L1)=-0.5 ! 插值区间左端点处的一阶导数值 CN=3.0*(Y1(N1)-Y1(N2))/(X1(N1)-X1(N2)) IF(N1-L1.EQ.1)GOTO 55 ! 插值区间只具有左右两个端点 L2=L1+1 DO 50 J=L2,N2 H1=X1(J)-X1(J-1) H=X1(J+1)-X1(J) AF=H1/(H+H1) ! 计算α值 BT=3.*((1.-AF)*(Y1(J)-Y1(J-1))/H1+AF*(Y1(J+1)-Y1(J))/H) !计算β值 A(J)=-AF/(2.+(1.-AF)*A(J-1)) ! 计算插值区间内插值节点处的一阶导数值 mj 50 B(J)=(BT-(1.-AF)*B(J-1))/(2.+(1.-AF)*A(J-1)) 55 A(N1)=(CN-B(N2))/(2.+A(N2)) 60 A(N2)=A(N2)*A(N2+1)+B(N2) ! 计算插值区间右端点处的一阶导数值 mn