写了个4节点等参元法求刚度矩阵的脚本，求确认优化

whlxy2002

问题：1，请达人优化一下吧，特别是求Br，Bl矩阵的代码
          2，k矩阵如何简化
         3， 节点位移为什么这么大呢 是不是程序不对
        4，节点应力的求法对吗？
        5，这样是不是计算的时候，就不用求整体刚度矩阵了，求出每个单元的 矩阵就行了？还是要合成一个大的刚度矩阵，比如有100个节点，就要合成100维的刚度矩阵？然后在计算节点位移？

代码

1. 
   
2. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%四节点等参元计算单元刚度矩阵
   
3. syms m n  %m n 为母单元的局部坐标的m n轴，x y全局坐标系的轴
   
4. N=zeros(2,8);%
   
5. N1=(1+m)*(1+n)/4;N2=(1-m)*(1+n)/4;N3=(1-m)*(1-n)/4;N4=(1+m)*(1-n)/4;% 形函数
   
6. N=[N1 0 N2 0 N3 0 N4 0;
   
7.    0 N1 0 N2 0 N3 0 N4];%形函数矩阵
   
8. e=[1 0 0 0;
   
9.    0 0 0 1;
   
10.    0 1 1 0];%应变=e*[du/dx du/dy dv/dx dv/dy]'
    
11. Nv=N(1,:);%中间变量
    
12. Nv(find(Nv(:)==0))=[];%
    
13. XY=[0 1 1 0;0 0 1 1]';%单元节点坐标XY=[x1 x2 x3 x4; y1 y2 y3 y4]'
    
14. Jdiff=jacobian(Nv,[m n]);%雅克比矩阵
    
15. Jdiff=Jdiff';%
    
16. J=Jdiff*XY;%雅克比方阵
    
17. InvJ=inv(J);%雅克比方阵的逆矩阵
    
18. NJ=J(1,1)*J(2,2)-J(1,2)*J(2,1);%雅克比方阵的模;符号矩阵求模怎么算呢？
    
19. Br=[diff(N(1,:),m);diff(N(1,:),n);diff(N(2,:),m);diff(N(2,:),n)];%B矩阵右侧矩阵，形函数对局部坐标分别求导
    
20. Bl=e*[InvJ(1,1) InvJ(1,2) 0 0 ;
    
21.       InvJ(2,1) InvJ(2,2) 0 0 ;
    
22.       0 0 InvJ(1,1) InvJ(1,2);
    
23.       0 0 InvJ(2,1) InvJ(2,2)];%B矩阵左侧矩阵
    
24. B=Bl*Br;%矩阵B，应变=B*[u1 v1 u2 v2 u3 v3 u4 v4]'
    
25. E=1000;mu=0.3;%定义杨氏模量Pa和泊松比
    
26. D=E/(1-mu^2)*[1 mu 0;mu 1 0; 0 0 (1-mu)/2];%弹性矩阵，平面应变问题应力应变关系矩阵
    
27. k=B'*D*B*NJ;%积分号内的矩阵，得出的是符号矩阵，里面的数都是用分数表示的，如何用有限位的小数表示呢？
    
28. % 例k(1,1)=25*(-4719744281318681/85899345920-4719744281318681/85899345920*conj(n))*(-1/20-1/20*n)+25*(-6607641993846153/343597383680-6607641993846153/343597383680*conj(m))*(-1/20-1/20*m)
    
29. 
    
30. Ke=double(int(int(k,m,-1,1),n,-1,1))%单元刚度矩阵
    
31. %%%单元等效节点荷载，根据虚功原理，单元荷载做得虚功=等效节点荷载做得虚功
    
32. q=[0 1*1*1*18]';%单元荷载q=[qx qy]'，此处假设qx=0，qy=单元体积*容重
    
33. pe=N'*q*NJ %等效节点荷载局部坐标积分号内的矩阵
    
34. Pe=double(int(int(pe,m,-1,1),n,-1,1))
    
35. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%单元节点位移
    
36. Dispe=Ke\Pe %单元节点位移，报警告：Warning: Matrix is close to singular or badly scaled.
    
37.              % Results may be inaccurate. RCOND = 1.408701e-018.
    
38.              % > In test at 36
    
39. %算出来的位移出奇的大
    
40. %Dispe = 1.0e+014 *  [1.7944  -0.1056 1.7944 -0.4222  2.1111  -0.4222
    
41. %2.1111  -0.1056]=[u1 v1 u2 v2 u3 v3 u4 v4]哪儿出错了呢？
    
42. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%节点应力，没看太明白怎么求,不确定求得对不对
    
43. Se=D*B*Dispe; %Se=[Sxx Syy Sxy]'
    
44. Sen=zeros(3,4);
    
45. m=1;n=1;Sen(:,1)=subs(Se);
    
46. m=-1;n=1;Sen(:,2)=subs(Se);
    
47. m=-1;n=1;Sen(:,3)=subs(Se);
    
48. m=1;n=-1;Sen(:,4)=subs(Se);
    
49. Sen
    
50. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    

复制代码