差分程序.txt

设 计未知表格行数 -未知表格列数 -系数矩阵 -常数数组 -表格间距 -边界点 值 -边界点 值 -边界点 值 -未知点处 值数组 所有点处 X方向正应一个解决一般弹性力学应力问题的程序

差分程序设计
理论应用力学03-1 

0302100108 

刘琨 

2005-12-7 

程序设计目标： 

设计一个解决一般弹性力学应力问题的程序。要求针对一道题能够在任意的精度下进行计算，即网格数目可以任意确定。给予一定的边界条件能计算出其上各点处的应力函数以及应力分量值。并且具有一定的移植性。

程序符号约定： -未知表格行数 -未知表格列数 -系数矩阵 -常数数组 -表格间距 -边界点 值 -边界点 值 -边界点 值 -未知点处 值数组 所有点处 X方向正应力 所有点处 Y方向正应力

程序设计思路：

 （一） 差分计算应力问题包括以下几个步骤： 

    1. 实体-几何-进行划分表格，确定未知点的数目。（手算）    

    2. 分析边界条件，确定边界基点，计算边界点处 、 、 值 

    3. 列出各未知点的平衡方程。（电算）

    4. 解线性方程组。（电算） 

    5. 对号入座，利用差分关系求出各点的应力值。（电算）

    6. 分析应力值，对整个实体做出评价。  

（二）系数矩阵A 经分析平衡方程系数矩阵A随未知数数目的增加有规律地变化。利用此规律设计出在不同表格数目下系数矩阵的生成程序。在不同的表格数下的系数矩阵有如下变化： m=2,n=2 m=2,n=3 m=2,n=3 m=3,n=2 由以上的矩阵变化可以看出系数矩阵的变化规律：对于m行n列的未知矩阵由m*m个子块构成，每个是n*n矩阵。子块分为三类： 单位矩阵 以上为系数矩阵的形成规律。

（三）手算得常数数组A1 

（四） 解矩阵方程 ，原本用C语言写解方程程序，由于时间原因，再此用MATLAB软件计算。

（五） 用差分原理由X得到应力矩阵 yx yy 以上就是整个程序的思路。

 弹性力学实例计算例：

     设有矩形截面简支梁，深度为 ，长度为 ，体力可以不计，在上面受有均布荷载 ，由两端的反力 维持平衡如下图所示，为了研究方便，取单位宽度梁来考虑。

    此题用半逆解法计算得： 以下用差分法进行计算（取5行6列为例）分格及标号如图（按要求布点，见程序） 计算边界点 （见程序）。运行程序以及结果见下 程序：

 %以下程序均在matlab实现 %要求左上角点为基点，排点依次从左往右从上往下进行 ，点号与之相对应。 

%边界参数

 m=5; %未知行数

 n=6; %未知列数

 h=1; %步长 

q=10; l=3.5; 

H=6; 

%上边界 

for i=1:n+1 fai(i)=-(i*h)^2*q/2; %定义边界应力函数数组

 daox(i)=-q*i*h; %定义边界应力函数对X求导数组 

daoy(i)=0; %定义边界应力函数对Y求导数组

 end

 %右边界

 for i=n+2:m+n+2

 fai(i)=-q*2*l^2;

 daox(i)=(i-(n+1))*h*q*l/H-q*2*l;

 daoy(i)=0;

 end 

%下边界

 for i=m+n+3:2*n+m+3

 fai(i)=-q*l*(i-(m+n+2))*h+q*((i-(m+n+2))*h)^2/2-q*(2*l-(i-(m+n+2))*h)^2/2; 

daox(i)=-q*l;

 daoy(i)=0;

 end 

%左边界

 for i=2*n+m+4:2*n+2*m+4 fai(i)=0;

 daox(i)=(i-2*n-m-3)*h*q*l/H-q*l;

 daoy(i)=0; 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为已知量（手算）

 if m>4&n>4 %未知行列均须大于3

 %定义单元矩阵E1 

for i=1:n

   for j=1:n 

     if i==j E1(i,j)=-8; 

       elseif i==j+1 E1(i,j)=2;

        elseif j==i+1 E1(i,j)=2; 

        else E1(i,j)=0; 

     end 

   end 

end 

% 定义单元矩阵E2(外）

  for i=1:n

    for j=1:n

      if i==j

           if i==n|i==1 

            E2(i,j)=22; 

          else E2(i,j)=21; 

          end

              elseif i==j+1

                        E2(i,j)=-8;

               elseif j==i+1 

                        E2(i,j)=-8;

               elseif i==j+2  

                        E2(i,j)=1; 

               elseif j==i+2 

                       E2(i,j)=1;

              else E2(i,j)=0; 

                 end 

           end

     end 

% 定义单元矩阵E3(内）

  for i=1:n 

      for j=1:n

            if i==j 

               if i==n|i==1 

                     E3(i,j)=21;

                          else E3(i,j)=20;

                 end 

              elseif i==j+1 

                     E3(i,j)=-8; 

              elseif j==i+1

                     E3(i,j)=-8;

               elseif i==j+2 

                      E3(i,j)=1; 

               elseif j==i+2 

                      E3(i,j)=1;

                else

                       E3(i,j)=0;  

               end 

           end 

      end 

%定义单元矩阵E4

      E4=eye(n); 

%主系数 

 for L1=1:n

      for L2=1:n 

           A(L1,L2)=E2(L1,L2);

       end 

   end 

  for L1=1:n 

        for L2=n+1:2*n

             L2y=L2-n;

             A(L1,L2)=E1(L1,L2y); 

      end

  end

  for L1=1:n

          for L2=2*n+1:3*n 

               L2y=L2-2*n; 

                A(L1,L2)=E4(L1,L2y);

      end

 end 

  for L1=(m-1)*n+1:(m-1)*n+n

           L1y=L1-(m-1)*n; 

         for L2=(m-2)*n+1:(m-2)*n+n 

            L2y=L2-(m-2)*n; 

           A(L1,L2)=E1(L1y,L2y); 

      end 

end 

    for L1=(m-1)*n+1:(m-1)*n+n 

           L1y=L1-(m-1)*n; 

    for L2=(m-1)*n+1:(m-1)*n+n 

           L2y=L2-(m-1)*n; 

             A(L1,L2)=E2(L1y,L2y);

         end

     end 

   for L1=(m-1)*n+1:(m-1)*n+n 

        L1y=L1-(m-1)*n;

       for L2=(m-3)*n+1:(m-3)*n+ n

          L2y= L2-(m-3)*n; 

          A(L1,L2)=E4(L1y,L2y); 

      end 

end 

   for k=2:m-1

       for L1=(k-1)*n+1:(k-1)*n+n 

           L1y=L1-(k-1)*n; 

       for L2=(k-2)*n+1:(k-2)*n+n

            L2y=L2-(k-2)*n; 

                A(L1,L2)=E1(L1y,L2y); 

           end

       end 

   for L1=(k-1)*n+1:(k-1)*n+n

           L1y=L1-(k-1)*n;

          for L2=(k-1)*n+1:(k-1)*n+n    

             L2y=L2-(k-1)*n;

            A(L1,L2)=E3(L1y,L2y);

     end

 end

   for L1=(k-1)*n+1:(k-1)*n+n 

          L1y=L1-(k-1)*n;

         for L2=k*n+1:k*n+n

           L2y=L2-k*n; 

                   A(L1,L2)=E1(L1y,L2y);

       end

 end 

  if k>2

         for L1=(k-1)*n+1:(k-1)*n+n 

            L1y=L1-(k-1)*n;

          for L2=(k-3)*n+1:(k-3)*n+n 

             L2y=L2-(k-3)*n; 

             A(L1,L2)=E4(L1y,L2y);

         end 

    end 

end