梁邊界元程序資料

1、6.12.1 梁靜態(tài)拉壓邊界元程序BEM1 % * 數(shù)據(jù)準(zhǔn)備 * % = 梁的幾何參數(shù) = L=1;D = 0.05;AREA = pi*D*D/4;% = 材料特性參數(shù) = E=200E+9;% = 梁的載荷條件 = NForce=1;Force(1:NForce)=100000;%力的大小，參數(shù)：力的序號F_Position(1:NForce)=0.5;%力的位置。參數(shù)：力的序號% = 梁的邊界條件 = IB(1:4)=0;U(1:4)=0;% 邊界編號% _% | POSITION | 受力（N） | 變形(m) |% -% | X=0(首端） | 1 | 3 |% | X=L(末端)

2、| 2 | 4 |% -U(2)=0;U(3)=0;IB(2)=1;IB(3)=1;% = 欲求解的內(nèi)點坐標(biāo) = N=1000; %將梁等分為N等份I=1:N+1;X(I)=L/N*(I-1);% 內(nèi)點中包含兩端邊界點% * 數(shù)據(jù)準(zhǔn)備結(jié)束 *%*生成力系統(tǒng)矩陣 G 陣*G(1:2,1:2) = 0;G=0 -L/(2*AREA*E);-L/(2*AREA*E),0;%* 生成位移系統(tǒng)矩陣 H 陣 *H(1:2,1:2) = 0;H=-0.5 0.5; 0.5 -0.5;%* 生成力矩陣 P 陣 *P(1:2,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*F_P

3、osition(I)/(2*AREA*E);P(2)=P(2)+Force(I)*(L-F_Position(I)/(2*AREA*E);end% *轉(zhuǎn)換成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% * 將已知與末知邊界所對應(yīng)的G，H陣分解成兩部分，其中AB陣為未知邊界所對應(yīng)的系數(shù)陣;*A=eye(2) -K;AB=;for j=1:4 if IB(j)=0 AB=AB A(:,j); else for i=1:2 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:4 if

4、 IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 計算內(nèi)點值 *UI=; %保存內(nèi)點值for I=1:N+1% X(I)=L/N*I; AA=-X(I)/(2*AREA*E) (X(I)-L)/(2*AREA*E) 0.5 0.5;-0.5/AREA 0.5/AREA 0 0; % PI(1)=0; PI(2)=0; for J=1:NForce PI(1) = PI(1) + (Force(J)*(abs(X(I)-F_Position(J)/(-2*AREA*E); PI(2) = PI(2) + (Force(J)*FUHAO(X(I),F_Pos

5、ition(J)/2); end UI = UI;X(I) (AA*U'+PI')'end% * 結(jié)果輸出 *message='邊界值： 力變形'reshape(U,2,2)'message='內(nèi)點值：位置坐標(biāo) 變形應(yīng)力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),'.-')title('圖1 撓度') xlabel('長度L / m') ylabel('變形U/ m') grid onsubplot(2,1,2),plot(UI(:,

6、1),UI(:,3),'*-'); title('圖2 撓角') xlabel('長度L / m') ylabel('撓角 / rad') grid onreturn;6.12.2 梁縱向振動特性邊界元程序BEM2% 梁的縱向振動特性邊界元法計算% 輸入梁的材料參數(shù)E=input('輸入梁結(jié)構(gòu)材料彈性模量(GPa)：E=');p=input('輸入梁結(jié)構(gòu)材料密度(g/cm3)：p=');E=1.E9*E;p=p*1000.; % 輸入梁的幾何結(jié)構(gòu)參數(shù)type=input(' 輸入梁的截面形狀

7、(圓形=1，矩形=2),type=');if type=1 D=input('輸入梁截面直徑(mm)：D='); L=input('輸入梁長度(mm)：L='); A=0.25*pi*D*D*1.E-6; L=L*1.E-3;elseif type=2 width=input('輸入梁截面的寬度(mm)：width='); high=input('輸入梁截面的高度（mm）:high='); L=input('輸入梁長度（mm）:L='); A=width*high*1.E-6; L=L*1.E-3;end%

8、 邊界條件輸入%MM=input('輸入約束的數(shù)目(< 2)：MM=');for k=1:MM M(k)=input('輸入約束的編號(1,2)：=');end% 輸入頻率掃描計算參數(shù)%N=input('輸入頻率掃描計算點數(shù)目，N=');W0=input('輸入起始頻率(Hz)，W0=');WW=input('輸入頻率間隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1); Lmd=(2*pi*freq(k)*(p/E)0.5; SL=sin(Lmd*L); CL=cos(Lm

9、d*L); AES=Lmd*A*E*SL; r(k,1,1)=-CL/AES; r(k,1,2)=-1/AES; r(k,2,1)=-1/AES; r(k,2,2)=-CL/AES; % % % % % % % % % % if MM > 0 if M(1)=1 RR(k)=-r(k,2,1)*r(k,1,2)/r(k,1,1)+r(k,2,2); elseif M(1)=2 RR(k)=r(k,1,1)+r(k,1,2)*(-r(k,2,1)/r(k,2,2); end endend%for kk=1:N if MM=1 xx=abs(RR(kk); elseif MM=0 xx=ab

10、s(r(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('頻率（Hz）');ylabel('柔度（dB）');title('梁縱向振動頻率特性曲線');for kk=1:N-2 if MM=0 if (r(kk,1,1)<r(kk+1,1,1)&(r(kk+1,1,1)>r(kk+2,1,1) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end elseif MM=

11、1 if (RR(kk)<RR(kk+1)&(RR(kk+1)>RR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend6.12.3 梁的扭轉(zhuǎn)邊界元程序BEM3% * 數(shù)據(jù)準(zhǔn)備 * % = 梁的幾何參數(shù) = L=1;D = 0.05;%直徑IP=pi*D4/32.;AREA = pi*D*D/4;% = 材料特性參數(shù) = GGP=80e+9;% = 梁的載荷條件 = NForce=1;Force(1:NForce)=100000;%力的大小，參數(shù)：力的序號F

12、_Position(1:NForce)=0.5;%力的位置。參數(shù)：力的序號% = 梁的邊界條件 = IB(1:4)=0;U(1:4)=0;% 邊界編號% _% | POSITION | 受力（N） | 變形(m) |% -% | X=0(首端） | 1 | 3 |% | X=L(末端) | 2 | 4 |% -U(2)=0;U(3)=0;IB(2)=1;IB(3)=1;% = 欲求解的內(nèi)點坐標(biāo) = N=1000; %將梁等分為N等份I=1:N+1;X(I)=L/N*(I-1);% 內(nèi)點中包含兩端邊界點% * 數(shù)據(jù)準(zhǔn)備結(jié)束 *%*生成力系統(tǒng)矩陣 G 陣*G(1:2,1:2) = 0;G=0 -L

13、/(2*GGP*IP);-L/(2*GGP*IP),0;%*生成位移系統(tǒng)矩陣 H 陣 *H(1:2,1:2) = 0;H=-0.5 0.5; 0.5 -0.5;%* 生成力矩陣 P 陣 *P(1:2,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*F_Position(I)/(2*GGP*IP);P(2)=P(2)+Force(I)*(L-F_Position(I)/(2*GGP*IP);end% *轉(zhuǎn)換成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% *將已知與末知邊界所對應(yīng)的G，H陣分解成兩部分，其中AB陣為未知邊界

14、所對應(yīng)的系數(shù)陣;A=eye(2) -K;AB=;for j=1:4 if IB(j)=0 AB=AB A(:,j); else for i=1:2 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:4 if IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 計算內(nèi)點值 *UI=; %保存內(nèi)點值for I=1:N+1 AA=-X(I)/(2*GGP*IP) (X(I)-L)/(2*GGP*IP) 0.5 0.5;-0.5 0.5 0 0; % PI(1)=0; PI(2)=0; fo

15、r J=1:NForce PI(1) = PI(1) + (Force(J)*(abs(X(I)-F_Position(J)/(-2*GGP*IP); PI(2) = PI(2) + (Force(J)*FUHAO(X(I),F_Position(J)/2); end UI = UI;X(I) (AA*U'+PI')'end% * 結(jié)果輸出 *message='邊界值： 力變形'reshape(U,2,2)'message='內(nèi)點值：位置坐標(biāo) 變形內(nèi)力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),&

16、#39;.-')title('圖1 撓度') xlabel('長度L / m') ylabel('變形U/ m') grid onsubplot(2,1,2),plot(UI(:,1),UI(:,3),'*-'); title('圖2 撓角') xlabel('長度L / m') ylabel('撓角 / rad') grid onreturn;6.12.4 梁扭轉(zhuǎn)振動邊界元程序BEM4% 梁的扭轉(zhuǎn)振動特性邊界元法計算% 輸入梁的材料參數(shù)%G=input('輸入梁結(jié)構(gòu)

17、材料剪切模量(GPa)：G=');p=input('輸入梁結(jié)構(gòu)材料密度(g/cm3)：p=');G=1.E9*G;p=p*1000.; % 輸入梁的幾何結(jié)構(gòu)參數(shù)%type=input(' 輸入梁的截面形狀(圓形=1，其它形狀=2),type=');if type=1 D=input('輸入梁截面直徑(mm)：D='); L=input('輸入梁長度(mm)：L='); Ip=pi*D4/32*1.E-12; L=L*1.E-3;elseif type=2 Ip=input('輸入梁截面的極慣性矩(mm4)：Ip=&

18、#39;); L=input('輸入梁長度（mm）:L='); L=L*1.E-3;end% 邊界條件輸入%MM=input('輸入約束的數(shù)目(< 2)：MM=');for k=1:MM M(k)=input('輸入約束的編號(1,2)：=');end% 輸入頻率掃描計算參數(shù)N=input('輸入頻率掃描計算點數(shù)目，N=');W0=input('輸入起始頻率(Hz)，W0=');WW=input('輸入頻率間隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1);

19、 Lmd=(2*pi*freq(k)*(p/G)0.5; SL=sin(Lmd*L); CL=cos(Lmd*L); AES=Lmd*Ip*G*SL; r(k,1,1)=-CL/AES; r(k,1,2)=-1/AES; r(k,2,1)=-1/AES; r(k,2,2)=-CL/AES; % % % % % % % % % % if MM > 0 if M(1)=1 RR(k)=-r(k,2,1)*r(k,1,2)/r(k,1,1)+r(k,2,2); elseif M(1)=2 RR(k)=r(k,1,1)+r(k,1,2)*(-r(k,2,1)/r(k,2,2); end ende

20、nd%for kk=1:N if MM=1 xx=abs(RR(kk); elseif MM=0 xx=abs(r(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('頻率（Hz）');ylabel('柔度（dB）');title('梁扭轉(zhuǎn)振動頻率特性曲線');for kk=1:N-2 if MM=0 if (r(kk,1,1)<r(kk+1,1,1)&(r(kk+1,1,1)>r(kk+2,1,1) str=num2str(freq(kk+1),'

21、 Hz' text(freq(kk+1),aa(kk+1),str); end elseif MM=1 if (RR(kk)<RR(kk+1)&(RR(kk+1)>RR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend6.12.5 梁的彎曲邊界元程序BEM5% * 數(shù)據(jù)準(zhǔn)備 * % = 梁的幾何參數(shù) = L = 0.3;I=pi*0.054/64.% = 材料特性參數(shù) = E=200E+9;EI=E*I;% = 梁的載荷條件 = NForce=1;

22、Force(1:NForce)=1000;%力的大小，參數(shù)：力的序號F_Position(1:NForce)=0.2;%力的位置。參數(shù)：力的序號% = 梁的邊界條件 = IB(1:8)=0;U(1:8)=0;% 邊界編號% _% | POSITION | 力（N）|彎矩(Nm)| 撓度(m)|撓角(rad)|% -% | X=0(首端） | 1 | 2 | 5 | 6 |% | X=L(末端) | 3 | 4 | 7 | 8 |% -U(7)=0;U(4)=0;U(5)=0;U(6)=0;IB(5)=1;IB(6)=1;IB(7)=1;IB(4)=1;% = 欲求解的內(nèi)點坐標(biāo) = N=10;

23、%將梁等分為N等份I=1:N+1;X(I)=L/N*(I-1);% 內(nèi)點中包含兩端邊界點% * 數(shù)據(jù)準(zhǔn)備結(jié)束 * % *生成力系統(tǒng)矩陣 G 陣*G(1:4,1:4) = 0;G=-L3/(6*EI), 0, 0, -L2/(4*EI); 0, L2/(4*EI), -L3/(6*EI), 0; 0, L/(2*EI), -L2/(4*EI), 0; L2/(4*EI), 0, 0, L/(2*EI);%* 生成位移系統(tǒng)矩陣 H 陣 *H(1:4,1:4) = 0;H=0.5, -L/2 -0.5, 0; -0.5, 0, 0.5 L/2; 0 ,-0.5, 0 ,0.5;0, 0.5, 0,

24、-0.5;%* 生成力矩陣 P 陣 *P(1:4,1)=0;for I = 1:NForceP(1)=P(1)+Force(I)*(2*L3+F_Position(I)3-3*L*F_Position(I)2)/(12*EI);P(2)=P(2)+Force(I)*(2*L3+(L-F_Position(I)3 -3*L*(L-F_Position(I)2) /(12*EI);P(3)=P(3)+Force(I)*( F_Position(I)*(-F_Position(I)+2*L)/(4*EI);P(4)=P(4)+Force(I)*(L-F_Position(I)*(-F_Positio

25、n(I)-L)/(4*EI);end% *轉(zhuǎn)換成FK*X+P的形式*H = -H;K = inv(G)*H;P = inv(G)*P;% *將已知與末知邊界所對應(yīng)的G，H陣分解成兩部分，其中AB陣為未知邊界所對應(yīng)的系數(shù)陣A=eye(4) -K;AB=;for j=1:8 if IB(j)=0 AB=AB A(:,j); else for i=1:4 P(i)=P(i)-A(i,j)*U(j); end endendResult=inv(AB)*P;ID=0;for j=1:8 if IB(j)=1 ID=ID+1; U(j)=Result(ID); endend% * 計算內(nèi)點值 *UI=;f

26、or I=1:N+1AA=-(2*L3+X(I)3-3*L*X(I)2)/(12*EI),-(X(I)*(X(I)-2*L)/(4*EI),-(2*L3+(L-X(I)3 -3*L*(L-X(I)2)/(12*EI), (L-X(I)*(L-X(I)-2*L)/(4*EI), -0.5, (X(I)-L)/2., -0.5, X(I)/2;(X(I)*(X(I)-2*L)/(4*EI), (X(I)-L)/(2*EI), -(L-X(I)*(L-X(I)-2*L) /(4*EI), -X(I)/(2*EI), 0, -0.5, 0, -0.5 ; PI(1)=0; PI(2)=0; for J

27、=1:NForce PI(1)=PI(1)+(Force(J)*(2*L3+(abs(F_Position(J)-X(I)3-3*L*(F_Position (J)-X(I)2)/(12*EI);PI(2)=PI(2) + (-Force(J)*(abs(F_Position(J)-X(I) *(abs(F_Position(J)-X(I)-2*L)/(-4*EI) *FUHAO(F_Position(J),X(I); end UI = UI;X(I) (PI'-AA*U')'end% * 結(jié)果輸出 *message='邊界值： 力變形'%Ureshap

28、e(U,4,2)'message='內(nèi)點值：位置坐標(biāo) 變形應(yīng)力'UIsubplot(2,1,1),plot(UI(:,1),UI(:,2),'.-')title('圖1 撓度') xlabel('長度L / m') ylabel('變形U/ m') grid onsubplot(2,1,2),plot(UI(:,1),UI(:,3),'*-'); title('圖2 撓角') xlabel('長度L / m') ylabel('撓角 / rad'

29、;) grid onreturn;6.12.6 梁彎曲振動邊界元程序BEM6% 梁的彎曲振動特性邊界元法計算clear all% % 輸入計算初始條件% 輸入梁的材料參數(shù)%E=input('輸入梁結(jié)構(gòu)材料彈性模量(GPa)：E=');p=input('輸入梁結(jié)構(gòu)材料密度(g/cm3)：p=');E=1.E9*E;p=p*1000.; % 輸入梁的幾何結(jié)構(gòu)參數(shù)%type=input(' 輸入梁的截面形狀(圓形=1，矩形=2),type=');if type=1 D=input('輸入梁截面直徑(mm)：D='); L=input(&

30、#39;輸入梁長度(mm)：L='); A=0.25*pi*D*D*1.E-6; II=pi*D4/64.*1.E-12; L=L*1.E-3;elseif type=2 width=input('輸入梁截面的寬度(mm)：width='); high=input('輸入梁截面的高度（mm）:high='); L=input('輸入梁長度（mm）:L='); A=width*high*1.E-6; II=width*high3/12*1.E-12; L=L*1.E-3;end% 邊界條件輸入%MM=input('輸入約束的數(shù)目(&

31、lt; 4)：MM=');for k=1:MM M(k)=input('輸入約束的編號(1,2,3,4)：=');endkkk=0; % 自由度已知元素的數(shù)目：kkk if MM=1 for k1=1:4 if M(1)=k1 kkk=kkk+1; MK(kkk)=k1; end endelseif MM=2 for k1=1:4 if M(1)=k1 & M(2)=k1 kkk=kkk+1; MK(kkk)=k1; end endelseif MM=3 for k1=1:4 if M(1)=k1 & M(2)=k1 & M(3)=k1 kkk=

32、kkk+1; MK(kkk)=k1; end endend % 輸入頻率掃描計算參數(shù)N=input('輸入頻率掃描計算點數(shù)目，N=');W0=input('輸入起始頻率(Hz)，W0=');WW=input('輸入頻率間隔（Hz），WW=');%for k=1:N freq(k)=W0+WW*(k-1); Lmd=(2*pi*freq(k)0.5*(A*p/E/II)0.25; SL=sin(Lmd*L); CHL=cosh(Lmd*L); CL=cos(Lmd*L); SHL=sinh(Lmd*L); COCH=E*II*(CL*CHL-1);

33、 r(1,1)=(SL*CHL-CL*SHL)/(Lmd3*COCH); r(1,2)=SL*SHL/(Lmd2*COCH); r(1,3)=(SL-SHL)/(Lmd3*COCH); r(1,4)=(CHL-CL)/(Lmd2*COCH); r(2,1)=r(1,2); r(2,2)=(SL*CHL+CL*SHL)/(Lmd*COCH); r(2,3)=(CL-CHL)/(Lmd2*COCH); r(2,4)=(SL+SHL)/(Lmd*COCH); r(3,1)=r(1,3); r(3,2)=r(2,3); r(3,3)=(SL*CHL-CL*SHL)/(Lmd3*COCH); r(3,4

34、)=-SL*SHL/(Lmd2*COCH); r(4,1)=r(1,4); r(4,2)=r(2,4); r(4,3)=r(3,4); r(4,4)=(SL*CHL+CL*SHL)/(Lmd*COCH); % % % % % % % % % % for k1=1:MM for k2=1:MM if MM>1 QR(k1,k2)=r(M(k1),M(k2); else QR(1)=r(M(k1),M(k2); end end % for k2=1:kkk QP(k1,k2)=r(M(k1),MK(k2); end end if MM>0 QRQP=inv(QR)*QP; end %

35、% % % % % % % % % % % for k1=1:kkk for k2=1:kkk if kkk>1 R(k1,k2)=r(MK(k1),MK(k2); else R(1)=r(MK(k1),MK(k2); end end % for k2=1:MM RRP(k1,k2)=r(MK(k1),M(k2); end end % % % % % % % % % % if kkk>1 RRR(k,:,:)=R-RRP*QRQP; elseif kkk=1 RRR(k)=R-RRP*QRQP; else RRR(k,:,:)=r; endend%for kk=1:N if kkk

36、>1 xx=abs(RRR(kk,1,1); elseif kkk=1 xx=abs(RRR(kk); else xx=abs(RRR(kk,1,1); end aa(kk)=log(xx);end figure;plot(freq,aa);xlabel('頻率（Hz）');ylabel('柔度（dB）');title(' 梁的彎曲振動頻率特性曲線');for kk=1:N-2 if kkk>1 | kkk=0 if (RRR(kk,1,1)<RRR(kk+1,1,1)&(RRR(kk+1,1,1)>RRR(kk+

37、2,1,1) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end elseif kkk=1 if (RRR(kk)<RRR(kk+1)&(RRR(kk+1)>RRR(kk+2) str=num2str(freq(kk+1),' Hz' text(freq(kk+1),aa(kk+1),str); end endend%NN=input('輸入梁結(jié)構(gòu)內(nèi)部等份數(shù),NN=');dlt=L/NN;kjj=input('輸入要繪制的振型曲線的數(shù)目: kj

38、j=');for kj=1:kjj fq=input('輸入固有頻率（Hz），fq='); Lmd4=(2*pi*fq)0.5*(A*p/E/II)0.25; SL4=sin(Lmd4*L); CHL4=cosh(Lmd4*L); CL4=cos(Lmd4*L); SHL4=sinh(Lmd4*L); COCH4=E*II*(CL4*CHL4-1); % r(1,1)=(SL4*CHL4-CL4*SHL4)/(Lmd43*COCH4); r(1,2)=SL4*SHL4/(Lmd42*COCH4); r(1,3)=(SL4-SHL4)/(Lmd43*COCH4); r(1,4)=(CHL4-CL4)/(Lmd42*COCH4); r(2,1)=r(1,2); r(2,2)=(SL4*CHL4+CL4*SHL4)/(Lmd4*COCH4); r(2,3)=(CL4-CHL4)/(Lmd42*COCH4); r(2,4)=(SL4+SHL4)/(Lmd4*COCH4); r(3,1)=r(1,3); r(3,2)=r(2,3); r(3,3)=(SL4*CHL4-CL4*SHL4)/(Lmd43*COCH4); r(3,4)=-SL4*SHL