写了个4节点等参元法求刚度矩阵的脚本，求确认优化
问题：1，请达人优化一下吧，特别是求Br，Bl矩阵的代码2，k矩阵如何简化
3， 节点位移为什么这么大呢 是不是程序不对
4，节点应力的求法对吗？
5，这样是不是计算的时候，就不用求整体刚度矩阵了，求出每个单元的 矩阵就行了？还是要合成一个大的刚度矩阵，比如有100个节点，就要合成100维的刚度矩阵？然后在计算节点位移？
代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%四节点等参元计算单元刚度矩阵
syms m n%m n 为母单元的局部坐标的m n轴，x y全局坐标系的轴
N=zeros(2,8);%
N1=(1+m)*(1+n)/4;N2=(1-m)*(1+n)/4;N3=(1-m)*(1-n)/4;N4=(1+m)*(1-n)/4;% 形函数
N=[N1 0 N2 0 N3 0 N4 0;
0 N1 0 N2 0 N3 0 N4];%形函数矩阵
e=[1 0 0 0;
0 0 0 1;
0 1 1 0];%应变=e*'
Nv=N(1,:);%中间变量
Nv(find(Nv(:)==0))=[];%
XY=';%单元节点坐标XY='
Jdiff=jacobian(Nv,);%雅克比矩阵
Jdiff=Jdiff';%
J=Jdiff*XY;%雅克比方阵
InvJ=inv(J);%雅克比方阵的逆矩阵
NJ=J(1,1)*J(2,2)-J(1,2)*J(2,1);%雅克比方阵的模;符号矩阵求模怎么算呢？
Br=;%B矩阵右侧矩阵，形函数对局部坐标分别求导
Bl=e*[InvJ(1,1) InvJ(1,2) 0 0 ;
InvJ(2,1) InvJ(2,2) 0 0 ;
0 0 InvJ(1,1) InvJ(1,2);
0 0 InvJ(2,1) InvJ(2,2)];%B矩阵左侧矩阵
B=Bl*Br;%矩阵B，应变=B*'
E=1000;mu=0.3;%定义杨氏模量Pa和泊松比
D=E/(1-mu^2)*;%弹性矩阵，平面应变问题应力应变关系矩阵
k=B'*D*B*NJ;%积分号内的矩阵，得出的是符号矩阵，里面的数都是用分数表示的，如何用有限位的小数表示呢？
% 例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)
Ke=double(int(int(k,m,-1,1),n,-1,1))%单元刚度矩阵
%%%单元等效节点荷载，根据虚功原理，单元荷载做得虚功=等效节点荷载做得虚功
q=';%单元荷载q='，此处假设qx=0，qy=单元体积*容重
pe=N'*q*NJ %等效节点荷载局部坐标积分号内的矩阵
Pe=double(int(int(pe,m,-1,1),n,-1,1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%单元节点位移
Dispe=Ke\Pe %单元节点位移，报警告：Warning: Matrix is close to singular or badly scaled.
% Results may be inaccurate. RCOND = 1.408701e-018.
% > In test at 36
%算出来的位移出奇的大
%Dispe = 1.0e+014 *[1.7944-0.1056 1.7944 -0.42222.1111-0.4222
%2.1111-0.1056]=哪儿出错了呢？
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%节点应力，没看太明白怎么求,不确定求得对不对
Se=D*B*Dispe; %Se='
Sen=zeros(3,4);
m=1;n=1;Sen(:,1)=subs(Se);
m=-1;n=1;Sen(:,2)=subs(Se);
m=-1;n=1;Sen(:,3)=subs(Se);
m=1;n=-1;Sen(:,4)=subs(Se);
Sen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
页:
[1]