【正文】 )。K21(1,14,i)=AG*Bh(1)。K21(2,4,i)=AG*Bh(1)。K21(2,8,i)=AG*Bh(1)。K21(2,18,i)=AG*Bh(1)。K21(2,22,i)=AG*Bh(1)。K21(2,12,i)=AG*Bh(1)。K21(2,14,i)=AG*Bh(1)。K21(3,4,i)=*AG*Bh(1)。K21(3,22,i)=*AG*Bh(1)。K21(3,12,i)=*AG*Bh(1)。K21(3,14,i)=*AG*Bh(1)。end%前5階自振圓頻率計算w=zeros(5,1)。bl=zeros(5,1)。w0=1。ZX=zeros(5,n+1)。i1=1。i2=5。for iii=i1:i2ss=iii。 w1=w0+。dw=1。%計算動力項for i=1:r Kw(:,:,i)=w1*w1*p*Kc。endJJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。while JJss w1=w1dw。 for i=1:r Kw(:,:,i)=w1*w1*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endwhile JJss1 w1=w1+dw。 for i=1:r Kw(:,:,i)=w1*w1*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endw2=w1+dw。for i=1:r Kw(:,:,i)=w2*w2*p*Kc。endJJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。 while JJss w1=w2。 w2=w2+dw。 for i=1:r Kw(:,:,i)=w2*w2*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endJJwhile abs(w2w1) w0=(w1+w2)/。 for i=1:r Kw(:,:,i)=w0*w0*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。 if JJss w2=w0。 else w1=w0。 endendw(iii)=w0。%bl(iii)=L*sqrt(w0*sqrt(RA(1)/Dw(1)))/。%bl(iii)=L*sqrt(w0*sqrt(RA(1)/Dw(1)))/。X=zhenxing1(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。for iki=1:2*n+1 ZX(iii,iki)=X(1,iki)/X(1,1)。endend%定義節(jié)點(diǎn)坐標(biāo)（考慮虛擬段以后）xn=L:L/n:0。for i=1:n1 yn(2*i)=xn(i+1)。 yn(2*i+1)=xn(i+1)。endyn(1)=xn(1)。yn(2*n)=xn(n+1)。yn(2*n+1)=xn(n+1)。%以下為輸出部分for i=i1:i2 i0=ii1+1。i00=i2i1+1。subplot(1,i00,i0)。plot(ZX(i,:),yn)。% plot(yn,X(2,:))。% plot(yn,X(3,:))。grid on。end附錄二（算例二源程序）%箱型彎扭插值方法分析程序clear%節(jié)點(diǎn)總個數(shù)s1=3。s2=89。%廣義位移總個數(shù)m=s1+s2。%結(jié)構(gòu)高度或跨度L=18。%結(jié)構(gòu)分段情況（數(shù)組長度為分大段數(shù)，數(shù)值為每大段分段數(shù)）Ir=[18]。%斷面中線每大段分段數(shù)Cr=[90]。%斷面中線總分段數(shù)DS=0。for i=1:length(Cr) DS=DS+Cr(i)。end%斷面中線每大段中小段長度Cd=[]。%斷面每大段壁厚Bh=[]。%材料密度p=2500。%每小段節(jié)點(diǎn)信息 %一般段的節(jié)點(diǎn)信息 for i=1:83 JD(1,i)=i。 JD(2,i)=i+1。 end for i=86:89 JD(1,i)=i1。 JD(2,i)=i。 end %特殊段的節(jié)點(diǎn)信息 JD(1:2,84)=[84,1]。 JD(1:2,85)=[19,85]。 JD(1:2,90)=[89,61]。%斷面中線每小段長度、厚度k=0。for i=1:length(Cr) for j=k+1:k+Cr(i) Ad(j)=Cd(i)。 At(j)=Bh(i)。 endk=k+Cr(i)。endAE=+10。AG=+10。%分布外荷載q=zeros(2,m)。% q(1,1)=100。% q(2,1)=100。%結(jié)構(gòu)分段的大段數(shù)r=length(Ir)。%結(jié)構(gòu)分段的小段總數(shù)n=0。for i=1:r n=n+Ir(i)。end%每小段的長度dL=L/n。%賦初值 %結(jié)構(gòu)端部支承的柔度和剛度矩陣 PL=zeros(m,m)。 SR=zeros(m,m)。 %中間支承的剛度矩陣（樓板的作用可看做支承） Sf=zeros(m,m,n)。%集中荷載賦值（按實(shí)際情況輸入） %fb(2,24)=1e+6。%fb(2,60)=1e+4%支承剛度矩陣賦值% for i=1:s1% PL(i,i)=1E12。% SR(i,i)=1E+12。% end% for i=s1+1:m% PL(i,i)=1E+12。% SR(i,i)=1E12。% end%計算結(jié)構(gòu)的 K22，K21，K11,Kc（按大段分別計算）K22=zeros(m,m,r)。K21=zeros(m,m,r)。K11=zeros(m,m,r)。Kc=zeros(m,m,r)。for i=1:rK22(1,1,i)=AG*Bh(1)*Cd(1)*72。K22(2,2,i)=AG*Bh(1)*Cd(1)*18。 K22(3,3,i)=AG*。Kc(1,1,i)=Bh(1)*Cd(1)*72。Kc(2,2,i)=Bh(1)*Cd(1)*18。 Kc(3,3,i)=。for ii=1:DS j=JD(1,ii)+s1。k=JD(2,ii)+s1。 K22(j,j,i)=K22(j,j,i)+1/3*AE*At(ii)*Ad(ii)。 K22(k,k,i)=K22(k,k,i)+1/3*AE*At(ii)*Ad(ii)。 K22(j,k,i)=1/6*AE*At(ii)*Ad(ii)。 K22(k,j,i)=1/6*AE*At(ii)*Ad(ii)。 K11(j,j,i)=K11(j,j,i)+AG*At(ii)/Ad(ii)。 K11(k,k,i)=K11(k,k,i)+AG*At(ii)/Ad(ii)。 K11(j,k,i)=AG*At(ii)/Ad(ii)。 K11(k,j,i)=AG*At(ii)/Ad(ii)。 Kc(j,j,i)=Kc(j,j,i)+1/3*At(ii)*Ad(ii)。 Kc(k,k,i)=Kc(k,k,i)+1/3*At(ii)*Ad(ii)。 Kc(j,k,i)=1/6*At(ii)*Ad(ii)。 Kc(k,j,i)=1/6*At(ii)*Ad(ii)。endK21(1,4,i)=AG*Bh(1)。K21(1,40,i)=AG*Bh(1)。K21(1,46,i)=AG*Bh(1)。K21(1,82,i)=AG*Bh(1)。% K21(1,100,i)=AG*Bh(1)。% K21(1,118,i)=AG*Bh(1)。K21(2,4,i)=AG*Bh(1)。K21(2,22,i)=AG*Bh(1)。K21(2,40,i)=AG*Bh(1)。K21(2,46,i)=AG*Bh(1)。K21(2,64,i)=AG*Bh(1)。K21(2,82,i)=AG*Bh(1)。% K21(2,85,i)=AG*Bh(1)。% K21(2,100,i)=AG*Bh(1)。K21(3,4,i)=*AG*Bh(1)。K21(3,40,i)=*AG*Bh(1)。K21(3,46,i)=*AG*Bh(1)。K21(3,82,i)=*AG*Bh(1)。end%前5階自振圓頻率計算w=zeros(5,1)。bl=zeros(5,1)。w0=1。ZX=zeros(5,n+1)。i1=1。i2=5。for iii=i1:i2ss=iii。 w1=w0+。dw=1。%計算動力項for i=1:r Kw(:,:,i)=w1*w1*p*Kc。endJJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。while JJss w1=w1dw。 for i=1:r Kw(:,:,i)=w1*w1*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endwhile JJss1 w1=w1+dw。 for i=1:r Kw(:,:,i)=w1*w1*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endw2=w1+dw。for i=1:r Kw(:,:,i)=w2*w2*p*Kc。endJJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。 while JJss w1=w2。 w2=w2+dw。 for i=1:r Kw(:,:,i)=w2*w2*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。endJJwhile abs(w2w1) w0=(w1+w2)/。 for i=1:r Kw(:,:,i)=w0*w0*p*Kc。 end JJ=bzzjs(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。 if JJss w2=w0。 else w1=w0。 endendw(iii)=w0。%bl(iii)=L*sqrt(w0*sqrt(RA(1)/Dw(1)))/。%bl(iii)=L*sqrt(w0*sqrt(RA(1)/Dw(1)))/。X=zhenxing1(L,m,r,n,Ir,K22,K21,K11,Kw,Sf,SR,PL)。for iki=1:2*n+1 ZX(iii,iki)=X(1,iki)/X(1,1)。endend%定義節(jié)點(diǎn)坐標(biāo)（考慮虛擬段以后）xn=L:L/n:0。for i=1:n1 yn(2*i)=xn(i+1)。 yn(2*i+1)=xn(i+1)。endyn(1)=xn(1)。yn(2*n)=xn(n+1)。yn(2*n+1)=xn(n+1)。%以下為輸出部分for i=i1:i2 i0=ii1+1。i00=i2i1+1。subplot(1,i00,i0)。plot(ZX(i,:),yn)。% plot(yn,X(2,:))。% plot(yn,X(3,:))。grid on。end41