理学有限单元法课后习题全部答案王勖成.pdf

习题习题 1.2： 在用有限元法求解时，边界条件总是满足的，控制方程的不完全匹配，会产生误差。题中所 给出的近似函数： 23 0123 aa xa xa x=+，应该满足边界条件，对于情况（1） ，代入边 界条件可得 2 12 03 3 1 0, a La L aa L =，从而 333 2 12 23 ()() xxx a xax LLL =+ （1） 上式中的最后一项 3 3 x L 前面没有待定系数，这是由于使用了在 x=L 处=1 的强制边界条件。 从物理意义上说，相当于给定边界条件的解为齐次方程的通解加一个特解的缘故。将（1） 式代入教材（1.2.26）式，得到残量： 12 23 66 ( )( 6)(2)( ) xxx R xaaQ x LLL =+ 不同的求解方法，如配点法、子域法和伽辽金法，只是残量在某种意义上某个区域加权积分 为零。 配点法强制残量 R(x)在有限个点严格为零，点的个数取决于未知数个数，这里为 2，通常取 所选的点在域内均匀分布，则取 x=L/3 和 x=2L/3 处，R(x)=0,这样得到 2 ( )0,()0 33 LL RR=，从而可以解出待定系数 12 ,a a。带入（1）式可以得到。 配点法仅考虑了有限个点的局部特性，子域法则要求在有限个子域 i 内残量的积分 ( )0 i R x dx = 为零，子域的个数仍然取决于未知函数个数，通常选取各子域的并集为整个 待求区域，一般情况可以选择各子域大小相同，但对于某些局部变化较复杂的区域，可以缩 小 子 域 的 大 小 ， 使 得 子 域 分 布 更 合 理 。 例 如 取 子 域 为 12 |0/2, |/2xxLx LxL = =，则利用 12 ( )0,( )0R x dxR x dx = ， 可以求出待定系数 12 ,a a。 伽辽金法作为加权余量法的特殊形式，权函数选择为插值函数 12 ,N N， 这里 33 2 12 2 ( ),( ) xx N xxNxx LL =，这样，利用( ) ( )0,1,2 i N x R x dxi = 可以求出待 定系数 12 ,a a。 对于其余边界条件情况可依此类推。 练习题 1.4，注意近似函数要满足边界条件，从而可知截面及坐标系如图所示： ，很多同学把积分区域弄错了，也有不少同学计 算错误。这里，由于边界为零，采用泛函及其弱形式得到的积分结果是相同的。最终计算得 到：a1=4608/(134), a2=-512/(154), a3=-1536/(854)。 练习题 1.5，泛函的欧拉方程基本没太多问题，泛函为零得到边界条件： 23 23 0 0 L d wdwd w w dxdxdx = 1.5 如有一问题的泛函为 2 22 2 0 ( ) 22 L EId wkw wqw dx dx =+ ，其中 E, I, k 是常数，q 是给定函数，w 是未知函数，试导出原问题的微分方程和边界条件. 22 22 0 ( ) L d wd w wEIkw wq w dx dxdx =+ 2223 2223 00 0 234 234 0 00 ()() () L LL LL L d wdwd w dwd w dw EIdxEIEIdx dxdxdxdxdxdx d w dwd wd w EIEIwEIwdx dxdxdxdx = =+ 224 224 00 23 23 00 ( ) () LL LL d wd wd w wEIkw wq w dxEIkwqwdx dxdxdx d w dwd w EIEIw dxdxdx =+=+ + 微分方程： 4 4 0 d w EIkwq dx += 边界条件： 22 22 0 0 xx L d wd w dxdx = =， 33 33 0 0 xx L d wd w dxdx = = 分强制边界和自然边界。 补充题 试作加权余量发的最小二乘配点法，并给出所得到的求解方程系数矩阵的特点分析。 （最小二乘配点法思路是， 利用使求解域内所选各点处误差平方的总和为最少的条件， 去建 立求解试函数系数的方程。配点法是强迫余量误差在所选点上为 0，最小二乘配点法则是余 量在所选点上的误差，满足平方和最小。 ） 解：近似函数为( )( ) ii u xN x a= ，不失一般性 余量为：( )( )( )( ) )( ) ii R xA uf xA N x af x= 最小二乘配点法取权函数 () () jiik j wA N axx a = 其中j=1,.,n; k=1,.,m 且mn 加权余量要求0 j w Rd = 1 11 ( ) ) () ( ) )( ) ( ) () ( ) )( ) () () )() () ()() Ka-P T jiikii j T jkii m T jkikik k mm TT jiij kk w RdAN x axxA N x af x d a ANxxxA N x af x d ANxA N x af x ANA N aANf = = = = = = = ()写成矩阵形式 因此， 1 () () m T ijjiji k kANA Nk = = ， 系数矩阵对称，且无需积分。 复习题 1.7 自然边界条件强制边界条件的区别何在?为什么这样命名？对于一个给定的微分方程，如何 区分这两类边界条件？ 自然边界条件与强制边界条件，二者都是针对边值条件来说的。边值条件一般有三类边界条件。第一类：狄里克莱(Dirichlet) 条件；第二类，诺依曼(Neumann)条件；第三类，前两者的混合条件，也叫洛平(Robin)条件 在选择近似函数时， 已经事先满足的边界条件为强制边界条件。 而自然边界条件则是在将等 效形式化为弱形式时包含在边界积分场上的边界条件。 对于 2m 阶微分算子，含 0 到 m-1 阶导数的边界条件称为强制边界条件，近似函数应该事先 满足。含 m 到 2m-1 阶导数的边界条件称为自然边界条件，近似函数不必事先满足。对于给 定的微分方程，判断其阶次，再依据边界所含导数阶数可区分两类边界。 思考题 1.8 泛函在什么条件下有极值？了解泛函是否有极值的意义何在？ 2 00 = ，则3n = 方向最高达 4 阶，214n ，则3n = 因此积分为3 39 = 采用优化积分方案，由于多项式中完备性的阶次为 0，所以用 1 1 （2）方向最高达 3 阶，213n ，2n = 方向最高达 4 阶，214n ，则3n = 因此积分为2 36 = 采用优化积分方案，由于多项式中完备性的阶次为 0，所以用 1 1 5.1 节点 4 和 5 只属于单元（I）和（II） ，因此需要引入多点约束方程 对于单元（III）可知 1 2 2 3 1 (1)(1)(1) 4 1 (1)(1)(1) 4 1 (1)(1) 2 N N N =+ + = = 那么在 1-4-3-5-2 这条边上 3 00 1 (,) ii i N = = 将点 4 和 5 坐标带入上式可得 3 4123 1 3 5123 1 1313 ( 1, )() 2884 1133 ( 1,) 2884 ii i ii i N N = = =+ + = = + 5.2 (1)理论解 单位质量的离心力 2 fr= 则体积力为 2 Fr= 于是 C 点右侧部分的惯性力之和为 222 2 3 d 8 l C l FrA rAl = 其应力值 22 3 8 C l =， 同样的方法可求 22 1 2 O l = 其实可以求出全场应力解 22 2 1 ( ) 2 lr l = 进而得到全场应变解 2 32 0 1 () 23 r udrrl r E =+ O C 由应变可积分得到位移解 22 2 1 ( ) 2 lr l = (2)有限元分析 对于一维 3 节点杆单元 插值函数如下 1 2 2 3 1 (1) 2 1 1 (1) 2 N N N = = =+ 那么 123 N N N r = BLN =，且由(1) 4 l r=+可知 4 rrl = 所以 411 2 22l =+BLN = 弹性=DE（弹性模量） 雅可比矩阵 3 1 11 () 0()2 22244 i i i Nlll Jx = = += 计算刚度矩阵： 1 1 1 22 1 2 2 1 22 111213 212223 313233 111 ()2 () 224 16111 2 ()()2 () 4222 11 2 ()4 42 e J d El d l KKK KKK KKK = =+ + = KBDB 节点编号 1 2 3 4 5 其中 1 2 3 为单元 1，3 4 5 为单元 2 同理可求 2 e K 集成总刚为 1 2 3 14162 00 333 3216 00 33 28162 333 3216 33 14 3 E l K =， 下面求等效节点载荷 1 1T 1 1 22 1 1 2 2 2 2 23 1 (1) 2 (1) 1 44 1 (1) 2 0 12 24 f J d ll d P l P P l = + = + = PN f 同理可求 1 2T 1 1 22 1 2 2 3 2 2 4 2 25 1 (1) 2 (3) 1 44 1 (1) 2 24 4 12 f J d ll d l P l P P l = + = + = PN f 集成节点载荷向量为 1 2 2 2 3 4 5 0 1 12 1 12 1 4 1 12 P P Pl P P = P 把 1 0u =带入约束条件，利用对角线改 1 法，将求解方程=K AP改写为 2 2 2 3 4 5 100000 32161 00 0 3312 281621 33312 32161 334 114 123 u E ul l u u = 求解可得到 2 2 32 3 3 4 5 000 470.1223958 880.2291666 384 1170.3046875 1280.3333333 u ll u EE u u = ， 由 2 32 0 1 () 23 r udrrl r E =+ 可得理论解 * 2 2 3 3 4 5 0 31 0 256 11 48 37 128 1 3 u l u E u u = 由所求出的位移解可以计算得到应力值 1 2 2 2 3 0.520830 0.458333 0.395833 l = 3 2 2 4 5 0.395833 0.208333 0.020833 l = 总体应力磨平 设 3 * 1 ii i N = =，其中 * i 为待求的改进后的节点应力值 将其带入 2 * 1 ()0 e i e V CN dV = = 即 12 *(1)*(2) 11 ()()0,1,2,3,4,5 ii VV N dVN dVi EE += 单元应力磨平 5.3 证明题 证明： （1）利用 4 个高斯点外推，课本的（5.3.22）到（5.3.25）公式。 （2）利用 4 结点应力插值 （5.3.21）式为 * ()0 (1,2,3,4) e i V N dVi= 采用2 2高斯积分方案： 2 * ,1 ()( ,) e kijij V i j N dVH F = = * ( ,)()() ijmmkm FNx = m x为高斯积分点， * k N为 4 结点的插值函数。 m 为已 知的 4 个高斯点应力 由于2 2高斯积分的权系数1 ij H =，并且 4 * 1 () mkmk k Nx = =， ,mI II II IV=代表高斯点的标号，1,2,3,4k =代表结点的标号。 代入积分公式得 2444 * ,1111 ()() )()()0 kmkmkmkmkmkm i jkmk NxNxNxNx = = ， 由此，可得， 写成矩阵形式 4 * I 1 * 4 1111* II * 1 2222 * 4 * 3333 III * 1 4444 * IV 1 () ( )()()() () ( )()()() ( )()()() () ( )()()() () kmj k kmj k kmj k kmj k Nx NINIINIIINIV Nx NINIINIIINIV NINIINIIINIV Nx NINIINIIINIV Nx = = = = 4 0 0 = 0 0 由于插值函数矩阵的非奇异性，因此必然有 4 * I 1 4 * II 1 4 * III 1 4 * IV 1 () 0 () 0 = 0 () 0 () kmj k kmj k kmj k kmj k Nx Nx Nx Nx = = = = ，移项，得 * I 11111 * II 22222 * III 33333 * IV 44444 ( )()()() ( )()()() ( )()()() ( )()()() NINIINIIINIV NINIINIIINIV NINIINIIINIV NINIINIIINIV = 同公式（5.3.23） 5.7 利用旋转周期结构的原理和方法求解图 5.24 所示的受风载的塔， 风载在 22 区 域呈余弦分布。 5.9 将（5.6.19）式的被积函数展开，证明平行四边形的非协调元能通过分片试验。 证明：对于平行四边形单元，Jconst= 所以只需要 11 11B 0d d = 即可 22 22 0 1010 BLN=0 0101 2000 0002 0220 = = 200 002 002 020 T B = 代入计算，可得 11 11B 0d d = 第八章约束变分原理第八章约束变分原理 8.1 用直接代入法和拉格朗日乘子法求解函数用直接代入法和拉格朗日乘子法求解函数 22 465zxxyy=+在约束条件：在约束条件：238xy 下的极值问题。并检查采用拉格朗日乘子法时修正函数下的极值问题。并检查采用拉格朗日乘子法时修正函数 * z或或 * z是否仍保持极值性。是否仍保持极值性。 解： （1）直接代入：将238xy代入 22 465zxxyy=+得： 2 52464zyy=+， 2 2 102410 100 dz y dy dz dy = = ，则12/5y =时z取驻值，此时2/5x =。 min 176/5z= （2） *22 465(238)zxxyyxy=+ 驻值条件： * * * 8620 61030 2380 z xy x z xy y z xy =+= =+= =+= 2/5 12/5 44/5 x y = = = * 176/5z= 二次型矩阵为： 862 6103 230 ， 11 80a=， 1112 2122 86 440, 610 aa aa = 862 6103400 230 = with (student): A:=Doubleint(1/64*(xi+1)*(eta+1)*(2+xi+eta-xi*2-eta*2)*2,xi=-11,eta=-11); := Ad -1 1 d -1 1 () + 1 () + 1 () + + 222 2 64 value(A); 4 15 ） 11 1211 11 167 3150 e ab b MWN N d d = W （程序： :=Doubleint(1/64*(xi+1)*(eta+1)*(2+xi+eta-xi*2-eta*2)*b*eta*(xi+1)*(eta+1)*2*(eta-1),xi=-11,eta=-1 1); := Bd -1 1 d -1 1 () + 1 2( ) + 1 3( ) + + 222b () 1 64 ? value(B); 167b 3150 ） - 1 - 计算力学作业 第 14 章 力 92 方楠 991788 11 1311 11 167 3150 e ac a MWN N d d = W， 11 1412 11 613 3150 e aa MWN N d d = W 以下的程序基本上是类似的，就不一一列出了。 11 1512 11 167 3150 e ab a MWN N d d = W， 11 1612 11 2 315 e ac a MWN N d d = W 其他的数据可以仿照这样的方法进行求解，这里就不一一计算了. 1 1 12 12 1 1 000 000 1 0004 000 e N N MW N N = 其中：N 1 100 000 000 = 对于集中质量矩阵，相当于每个结点集中 1/4 的质量，即： 14.3 用中心差分法、集中质量矩阵求所示变截面杆在图 14.8 所示外载作用下的响应（初 始条件：u x） 。 ( ,0)( ,0)0u x = 解：设 E 为材料的弹性模量。采用轴力杆单元，两个结点。易知将质量集中在结点，得到： 集中质量矩阵：MA。因此，集成集中质量矩阵：。 1/20 01/2 e l = 4 1/200 03/40 001/ M = 刚度矩阵： 11 1 1 e EA K l = 。因此集成刚度矩阵： 440 462 022 KE = 因此运动方程： 3 3 ( )1/200440 03/404620 001/4022( ) Q t aEa Q t += 。 由于最左边的点位移始终为 0，因此可以省略不计，将其变为二维的方程组： 3 03/4062 ( )01/422 aEa Q t += 。 初始条件：( ,0)( ,0)0u xu x =，时间步长选择：每个部分的振动周期： 22 m Tl k =/E，因此时间步长：/ n T tE =。取 / 2 E t = - 2 - 计算力学作业 第 14 章 力 92 方楠 991788 在时间段内： 0 tt 0 4 (/)E =c， 1 1 / c E =， 2 8 / c E =， 3 / 8 E c =。 由初始条件知道：； 1 0 4 a Q = 1 000 /1 00482 t E QQE 0 at =+= tt ； ? 01 30 01 Mc McCE =+= 。 ? 201 00130 ()()2 1001 tttttt QQKc M ac McC aEaEa Q =+ t (*) 求t时刻的位移： + （*） ? 3000130 2 011001 tttttt EaQEaEa Q + =+ 为了简便起见编写程序计算：Fortran 语言 program main implicit none integer:n double precision,dimension(2):a0=0 !t-t 位置 double precision,dimension(2):a1=0 !t 位置 double precision,dimension(2):a2=0 !t+t 位置 double precision:tem1,tem2 a0(1)=0.0d0 a0(2)=0.5d0 a1(1)=0.0d0 a1(2)=0.0d0 !初始值 do n=1,15 a2(1)=(0.0+2.0*a1(2)-3.0*a0(1)/3.0d0 a2(2)=1.0+2.0*a1(1)-1.0*a0(2) write(*,*)n write(*,*)a2 a0=a1 ! a1=a2 !传递给下一个时刻 enddo end program 仿照书上得到下表：(省略了 /E Q，只有系数)。 时间 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 t 10 t 0 0.33 0.667 0.4440.2220.003-0.1480.196 0.543 0.498 - 3 - 计算力学作业 第 14 章 力 92 方楠 991788 t a 0.5 1.00 1.167 1.3330.7220.1110.352 0.593 1.043 1.494 在t时间段内：我认为这段时间内不知道其初始的速度和加速度，因此缺少必要的初始 0 t 条件，故不能用中心差分法计算。若想计算，给出以上条件后，按照同样的步骤进行计算 即可。就不再详细叙述了。 14.4 用 Newmark 方法求解题 14.3 。 解：Newmark 方法是无条件稳定的，因此t 可以任意的选择。仍取： / 2 E t =，0.5,0.25= 。 运动方程： 3 03/4062 ( )01/422 aEa Q t += 。 由初始条件知道：； 0 0 0 a = 0 0 0 a = ； 1 0 4 a Q = ； 0 16 / c E =， 1 4 / c E =， 2 8 / c E =， 3 1c =， 4 1c =，c50=， 6 / 4 E c =， 7 / 4 E c = 在时间段内： 0 tt ? 01 623/40182 16 2201/426/ KKc McCEE E =+=+= ? 023145 ()( 03/40 168 () 01/4(/)/ ttt tttttt tt t QQM c ac ac aC c ac ac a aaa QEE + =+ =+ )t (*) 最终：? ? 18203/40 168 () 2601/4(/)/ tt ttttt KaEaaaQ QEE + =+= ； ? 1 1 62 1 218104 KE = 。 计算t时间的加速度和速度： t+ - 4 - 计算力学作业 第 14 章 力 92 方楠 991788 023 67 168 ()() (/)/ / 44 tttttt tttttt tttttttttt ac aac ac aaaa EE EE aac ac aaaa + + + = =+=+ a 为了简便起见我编写了一个 Fortran 程序计算之： program main implicit none integer:n double precision,dimension(2):a1=0,a2=0,a3 !t 时刻的位置，速度，加速度 double precision,dimension(2):b1=0,b2=0,b3=0 !t+t 时刻的位置，速度，加速度 double precision:tem1,tem2 a3(1)=0.0 a3(2)=4.0 !加速度初始值 do n=1,15 tem1=0+0.75d0*(16.0*a1(1)+8.0*a2(1)+a3(1) tem2=1+0.25d0*(16.0*a1(2)+8.0*a2(2)+a3(2) b1(1)=1.0d0/104.0d0*(tem1*6.0d0+tem2*2.0d0) b1(2)=1.0d0/104.0d0*(tem1*2.0d0+tem2*18.0d0) b3=16.0*(b1-a1)-8.0*a2-a3 b2=a2+0.25*a3+0.25*b3 write(*,*)n write(*,*)b1 a1=b1 a2=b2 a3=b3 enddo end program 下表省略了 /E Q，只有系数。 时间 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 t 10 t 0.038 0.213 0.508 0.6570.4350.029 -0.14 8 0.014 0.237 0.451 t a 0.346 0.994 1.307 1.1730.8940.6080.2940.173 0.553 1.212 在t时间段内，同理可以求得，这里就不仔细写出过程了。 0 t 14.5 用振型叠加法求解题 14.3。 解：运动方程： 3 03/4062 ( )01/422 aEa Q t += 。 - 5 - 计算力学作业 第 14 章 力 92 方楠 991788 按照一般的线性代数知识： 2 11 248 322 (/)3.38(/), 33 T EE = ， 2 21 248 322 (/)12.62(/), 33 T EE + = 。 最后转化成方程： 1 1 2 3.38(/)aEaQ +=， 2 2 2 12.62(/)aEaQ += = 。 初始条件： 经过转化后：( ,0)( ,0)0u xu x = 11 ( ,0)( ,0)0x xx x =，。 22 ( ,0)( ,0)0xxxx = 利用无阻尼情况下的杜哈梅积分公式： 0 1 ( )( )sin()sincos t iiiiii i i x trtdatbt =+ ； 最后结果为： 0 tt时候： 0 22 0 1 ( )sin()cos()|1 cos t t ii iiii ii RR i i x tRtdtt = 分量： 1( ) 0.41841 cos 3.38(/)x tQE E t= ， 2( ) 0.11211 cos 12.62(/)x tQ E Et= 。 最后振型叠加： 22 0.41841 cos 3.38(/) 33 ( ) 22 0.11211 cos 12.62(/) 0.3416 1 cos 3.38(/)0.0915 1 cos 12.62(/) 0.5917 1 cos 3.38(/)0.1585 1 cos 12.62(/) QEt E a t QEt E EtE Q E t EtE = = + t 取 / 2 E t = 计算几个数值：(省略了 /E Q，只有系数)。 时间 t 2 t 3 t 4 t 5 t 6 t 7 t 8 t 9 t 10 t 0.024 0.257 0.619 0.6060.211 -0.02 6 0.0030.080 0.301 0.625 t a 0.432 1.052 1.207 1.1510.9540.3770.0080.478 1.143 1.252 0 tt时候：可以按照上面的方法求出时刻的，有因为这个时候 Q0，所以： 0 t i ( ), ( )a t a t ( )sincos iiii x tatbt=+ ，带入初始条件即可，这里就不计算了。 - 6 - 计算力学作业 第 14 章 力 92 方楠 991788 最后，我们将三种情况在t时候做一个比较，可见各种方法均能得到近似正确的结果。 0 t 近似的方法得到的解在真解附近振荡。如果想得到更精确的解可以适当的减小步长。 三种方法比较 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 12345678910 步数 最右端位移 振型叠加法 Newmark方法 中心差分法 - 7 - 计算力学作业 第 15 章 力 92 方楠 991788 计算力学（2）作业 第 15 章 材料非线性问题的有限单元法 清华大学工程力学系 力 92 方楠 991788 15.1 一维弹塑性问题如图 15.25（a）所示，作用于中面截面得轴向力，材料性质如图 12.25(b)，分别用直接叠代法，NR 法和 mNR 法求解。 （ 30P = 1 12 AA= ） 。 解：受力方向如图所示：左边为 1，右边为 2。设 12 , 均大于 0。并且有 12 =。 由于存在着位移和应变关系l=，所以这里都是取应变为变量，最后进行代换即可。 应力应变有如下的关系： 3 11 1 3 1 (1 10 10(1 10 ) E = ) P ， 。 3 22 2 3 2 (3 10 30(3 10 ) E = 根据平衡关系我们有： 12 ()0 n PKf=+=，其中 12 (),30 n Kf= = 。(*) 因此可以根据以上的式子求解： 先设 0 12 1 10 3 = 带入 （*） 并按照分段函数进行求解： 3 33 3 1 30(1 10 ) 2 1 20(1 103 10 ) (3 10 ) n E E = 111 1 23003001.5 10 2 E 3 E =。继续求解 A、 直接叠代法 当时左端进入塑性。因此：。 1 1.5 10 = 3322 2002.0 10E = 继续带入发现值不变，因此最后的解就是 3 2.0 10 =。中间点的位移0.02l= B、 NR 法 根据上面的分析有： - 1 - 计算力学作业 第 15 章 力 92 方楠 991788 3 33 3 1 (230)(1 10 ) 2 1 (20)(1 103 10 ) (3 10 ) n E E E E = （*） 简单的设。因此带入 0 0= 00 1 (230)1.5 10 2 E E 3 = ， = 100 1.5 10 3 =+=。 11 1 (20)0.5 10E 3 E = = ， 211 2 10 3 =+=。 再次带入发现 0 = ，即达到稳定，因此最终的结果：。中间点的位移 3 2.0 10 = 0.02l= C、 mNR 法 首先设。 0 0= 0 2 n TT KK= E ，因此 3 33 3 1 (230)(1 10 ) 2 1 (20)(1 103 10 ) 2 (3 10 ) n E E E E = 简单的设。因此带入 0 0= 00 1 (230)1.5 10 2 E 3 E = ， = 100 1.5 10 3 =+=。 11 1 (20)0.25 10 2 E 3 E = = ， 211 1.75 10 3 =+=。 22 1 (20)0.125 10 2 E 3 E = = ， 322 1.875 10 3 =+=。 33 1 (20)0.0625 10 2 E 3 E = ，= 433 1.9375 10 3 =+=。 44 1 (20)0.03125 10 2 E 3 E = = ，。 544 1.96875 10 =+= 3 55 1 (20)0.015625 10 2 E 3 E = ，。 = 655 1.984375 10 =+= 3 66 1 (20)0.0078125 10 2 E 3 E = ，。 = 766 1.9921875 10 =+= 3 77 1 (20)0.00390625 10 2 E 3 E = = ，。 877 1.99609375 10 =+= 3 88 1 (20)0.001953125 10 2 E 3 E = = ，。 988 1.998046875 10 =+= 3 3 988 ()/0.098%0.1%= 可以认为满足收敛条件了。 因此最终的解：。中间点的位移 3 1.998 1010 =0.01998l=。 - 2 - 计算力学作业 第 15 章 力 92 方楠 991788 15.2 用增量法求解 15.1。对于以下两种方案： （1） （2）01 0152025 3030 30 624 分别用有平衡校正和无平衡校正得 Eular 法计算。 解： （1） 0152025 【A】 、有平衡校正： 正常的增量法表示为： 11 11 ()() ( mmTmmTmmm) KfKff + = = 。 取1/2=， 11 1 (1)() 2 mmmm + =+=+ m+ ()( TmTm KK，) + = 11 11 ()()( mmTmmTmm ) m KfKff + = = + 。 下面开始计算，采用 NR 方法叠代求解： 3 33 2(1 10 ) (1 103 10 ) E K E = 01：5 13 0100010 1 0, ()()0.75 10 2 T Kfff E = = =， 3 001 1 ( )0.375 10 2 + =+=。 13 10010 1 ()()0.75 10 2 Tm Kfff E + = = = 因此。 3 1 0.75 10 = 1 111 ()( () nTn mnmmm 10) KP + = +f + 。再次叠加发现不变，即为最后解。 15： ，20 31 1211121 1 0.75 10 , ()()0.25 10 2 T Kfff 3 E = = = 3 112 1 ( )0.5 10 2 + =+=。 13 21121 1 ()()0.25 10 2 Tm Kfff E + = = = 因此。 3 2 1.00 10 = 1 111 ()( () nTn mnmmm 10) KP + = +f + 。再次叠加发现不变，即为最后解。 ：2025 31 2322232 1 1.00 10 , ()()0.25 10 2 T Kfff 3 E = = = 3 223 1 ( )0.625 10 2 + =+=。 13 32232 1 ()()0.25 10 2 Tm Kfff E + = = = 因此。 叠加公式： 3 3 1.25 10 = 1 111 ()( () nTn mnmmm 10) KPf + = + + 。 （a） 1111 333103 1 () ( ()(15)0.25 10 T nm KPfE E + = += = 3 23 3 1.5 10 = ，。 叠加： （这里要注意 f 的定义，不仅仅简单的是载荷值，必须写出平衡方程才可以） 。 （b） 2122 333103 1 () ( ()(15) T nm KPfE E + = += = 3 1.5 10= 0 ，这时叠代不变了。因此最后的结果 。 3 - 3 - 计算力学作业 第 15 章 力 92 方楠 991788 ：2530 31 3433343 1 1.5 10 , ()()0.5 10 T Kfff 3 E = = = 3 334 1 ( )1.00 10 2 + =+=。 13 43343 1 ()()0.5 10 Tm Kfff E + = = = 因此。 3 4 2.0 10 = 1 111 ()( () nTn mnmmm 10) KP + +f + = 。再次叠加发现不变，即为最后解。 最后，整理结果，中点处的位移变化为：0。 0.00750.010.0150.02 【B】 、无平衡校正： 正常的增量法表示为： 11 11 ()() ( mmTmmTmmm) KfKff + = = 。这些和上面一样， 不单独列出了 01：5 13 0100010 1 0, ()()0