梁边界元程序解读.doc

梁静态拉压边界元程序BEM1 % ************************* 数据准备 ********************************* % ======================== 梁的几何参数 ============================ L=1; D = 0.05; AREA = pi*D*D/4; % ======================== 材料特性参数 ============================ E=200E+9; % ======================== 梁的载荷条件 ============================ NForce=1; Force(1:NForce)=100000;%力的大小，参数：力的序号 F_Position(1:NForce)=0.5;%力的位置。参数：力的序号 % ======================== 梁的边界条件 ============================ 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=1000; %将梁等分为N等份 I=1:N+1; X(I)=L/N*(I-1);% 内点中包含两端边界点 % ************************ 数据准备结束 *********************************** %　*********************　生成力系统矩阵 G 阵　*************************** G(1:2,1:2) = 0; G=[0 -L/(2*AREA*E);-L/(2*AREA*E),0]; %　************************ 生成位移系统矩阵 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:NForce P(1)=P(1)+Force(I)*F_Position(I)/(2*AREA*E); P(2)=P(2)+Force(I)*(L-F_Position(I))/(2*AREA*E); end % ***********************　转换成F＝K*X+P的形式　************************ H = -H; K = inv(G)*H; P = inv(G)*P; % **** 将已知与末知边界所对应的G，H阵分解成两部分，其中AB阵为未知边界所对应的系数阵;**** 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 end end Result=inv(AB)*P; ID=0; for j=1:4 if IB(j)~=1 ID=ID+1; U(j)=Result(ID); end end % **************************** 计算内点值 ***************************** UI=[]; %保存内点值 for I=1:N+1 % X(I)=L/N*I; AA=[-X(I)/(2*AREA*E) (X(I)-L)