精仪学院陶会荣1011202017临界力求解的直接迭代算法算法原理.DOC

精仪学院 陶会荣 1011202017 临界力求解的直接迭代算法 一、 算法原理 针对变截面压杆临界载荷的问题，基于一阶临界载荷作用下的结构稳定平衡条件，可以进行杆件的离散化，以杆件各离散点的未知挠度和压杆临界力为未知量，以每个离散点的挠度与临界力所满足的差分方程，构建求解临界力非线性方程组。 应用材料力学理论中，针对常截面，可以通过欧拉公式Fcr=π2EI/(μl)2 和经验公式Fcr=a-bλ来确定临界力。但是受压杆截面为变截面时，没有相应临界力表达式。而迭代算法可以应用于该领域的计算。 算例：如图所示两端铰支受压杆，压杆长l，截面为圆台形，材料弹性模量E，确定临界力。 图示1 问题分析：当压杆处于稳定状态的时候，应满足: 在该临界力的作用下，杆件处于弯曲状态下的平衡状态。若计算工程实际中的临界力，也就是杆件的一阶临界力，杆件正好处于图示2所示状态。设未知的挠曲线方程为y=f(x)。 图示2 若将杆件沿横坐标x轴离散，取Δx=l/n，离散横坐标为xi=iΔx，i=0,1,2,…n。对应挠度线的纵坐标为yi，i=0,1,2,…n,(xi，yi)，i = 0,1,2,…n为离散点。对于简支结构y0=yn=0。 针对压杆挠曲轴线微分方程: EIy″=M(x)。应用差分原理，y″i=yi+1-2yi + yi-1/Δx2，代入变截面条件I(xi)=Ii转化为: (1) 方程(1)为含yi，i=1,2,…,n-1和Fcr未知量为n个的非线性方程组。由于未知量数量大于方程的数量，因此不能直接求解临界力。在临界力作用下，挠度线存在一种形态，挠度线y 值是不定的。但是若给出其中的一点xm处的挠度ym，则变形曲线就是唯一确定的。此时共有yi，i= 1,2,…,n-1,i≠m和临界力Fcr的n-1个未知量，可以由方程(1)唯一确定。 迭代算法推导：设，,则公式（1）转化为： (2) 设其中一点为常数ym=c,获得i=0,1,2,…,n-1,n的迭代方程： 针对yn+1让公式(2)求和得到： (3) 获得收敛迭代式： 约束条件： 解本题： 由于 可知求出即可求出 则 设m=n/2,初始值。分段数n=42，迭代精度为e-6。 二、程序及结果分析 其MATLAB程序如下： function [k,r]=ddfch(n,x0,r0,eps) %k为迭代次数，r为迭代结果，n为分段数，n为偶数，x0为迭代初始值，r0为的值，exp %为迭代精度 n=42; x0=0.2; r0=1; eps=0.000001; r=0; for i=1:1:n+2 y(i)=x0; end y(1)=0; y(n+1)=0; y((n+2)/2)=r0; k=1; tol=1; while toleps for i=2:1:n a(i)=(y(i+1)+y(i-1))*((1+(i-1)/n)^4); end b=0; c=0; for i=2:1:n b=b+((y(i+1)-2*y(i)+y(i-1))*((1+(i-1)/n))^4); c=c+y(i); end d=b/c; for i=2:1:n z(i)=a(i)/(d+2*(1+(i-1)/n)^4); end z(1)=0; z((n+2)/2)=1; z(n+1)=0; z(n+2)=-d; for i=1:1:n+2 e=norm(y(i)-z(i))； end tol=max(e)； k=k+1; if(n100000) disp(迭代次数过多，方程可能不收敛); return; end for i=1:1:n+2 y(i)=z(i); end end r=y(n+2); 迭代结果：k=668，r=0.0225 通过以上迭代程序结果，可见迭代算法应用于工程结构的临界载荷分析计算的正确性和实用性。