Matlab求解平面应力问题

danielree

最近在学习有限单元法。但是感觉单纯看书太枯燥了，因此尝试通过编程来加深理解。在这里贴出学习求解平面应力问题时所编的程序供大家批评指正。希望能够有更多人来讨论。

形函数：输入面积坐标，输出形函数和形函数对自然坐标的偏导数。

function [shape,pdnatural] = shapetri3(l1, l2, l3)

%------------------------------------------------------------------------------

if nargin == 3

   if l1+l2+l3 ~= 1

      error('The sum of area coordinates mustequal 1!');

   else

      % No action

   end

elseif nargin == 2

   l3 = 1 - l1 - l2;

else

   error('Input error!');

end

%------------------------------------------------------------------------------

syms L1 L2 L3;

N1 = L1;

N2 = L2;

N3 = L3;

shapef = [N1; N2; N3];

pdnaturalf =[diff(N1,L1)-diff(N1,L3), diff(N2,L1)-diff(N2,L3), ...

              diff(N3,L1)-diff(N3,L3);...

              diff(N1,L2)-diff(N1,L3),diff(N2,L2)-diff(N2,L3), ...

              diff(N3,L2)-diff(N3,L3)];

% ------------------------------------------------------------------------------

shape =subs(shapef,{L1,L2,L3},{l1,l2,l3});

pdnatural =subs(pdnaturalf,{L1,L2,L3},{l1,l2,l3});            

end

雅可比矩阵：

function [jm, invjm,pdc] = jacobimatrix(pdn, node)

jm = pdn * node;

invjm = inv(jm);

pdc = invjm * pdn;

end

单元刚度矩阵，type选项：1为平面应力为题，其他为平面应变问题。

function ke =elemstiffpl(E, Nu, t, pdc, jm, type)

B1 = [pdc(1,1) 0; 0pdc(2,1);...

      pdc(2,1) pdc(1,1)];

B2 = [pdc(1,2) 0; 0pdc(2,2);...

      pdc(2,2) pdc(1,2)];      

B3 = [pdc(1,3) 0; 0pdc(2,3);...

      pdc(2,3) pdc(1,3)];      

B = [B1 B2 B3];

switch type

case 1                                            

   D = E/(1-Nu^2)*[1 Nu 0; Nu 1 0; 0 0(1-Nu)/2];

otherwise

   D = E/(1+Nu)/(1-2*Nu)*[1-Nu Nu 0; Nu 1-Nu 0;0 0 (1-2*Nu)/2];

end

ke=B'*D*B*t*1/2*det(jm);

end

集成总体刚度矩阵。

function Ke =assemblestiffpl(Ke, ke, i, j, k)

Ke([2*i-1 2*i 2*j-1 2*j2*k-1 2*k],[2*i-1 2*i 2*j-1 2*j 2*k-1 2*k]) = ...

Ke([2*i-1 2*i 2*j-1 2*j2*k-1 2*k],[2*i-1 2*i 2*j-1 2*j 2*k-1 2*k]) + ke

end

1.1   验证

n=[-1.000000,   -1.000000;

    1.000000,  -1.000000;

   -0.333333,  -1.000000;

    0.333333,  -1.000000;

    1.000000,   1.000000;

    1.000000,  -0.333333;

    1.000000,   0.333333;

   -1.000000,   1.000000;

    0.333333,   1.000000;

   -0.333333,   1.000000;

   -1.000000,   0.333333;

   -1.000000,  -0.333333;

   -0.422650,   0.000000;

    0.000000,   0.555556;

    0.415812,   0.018987;

   -0.001292,  -0.537154;

   -0.551425,  -0.577145;

    0.549601,  -0.573348;

    0.549829,   0.581575;

   -0.551197,   0.577778;];

e=[    18.00,      16.00,        4.00;

       18.00,       15.00,       16.00;

       19.00,       15.00,        7.00;

       19.00,       14.00,      15.00;

       20.00,       14.00,       10.00;

       20.00,       13.00,       14.00;

       17.00,       16.00,       13.00;

       17.00,        3.00,       16.00;

       13.00,       20.00,       11.00;

       20.00,        8.00,       11.00;

       10.00,        8.00,       20.00;

       14.00,       19.00,        9.00;

       19.00,        5.00,        9.00;

        7.00,        5.00,       19.00;

       15.00,       18.00,        6.00;

       18.00,        2.00,        6.00;

        4.00,       2.00,       18.00;

       12.00,       17.00,       13.00;

       17.00,        1.00,        3.00;

       12.00,        1.00,       17.00;

       16.00,       15.00,       13.00;

       14.00,       13.00,       15.00;

        3.00,        4.00,       16.00;

        6.00,        7.00,       15.00;

        9.00,       10.00,       14.00;

       11.00,       12.00,       13.00;];

nelem=size(e,1);

nnode=size(n,1);

K=zeros(nnode*2);

E = 2.1e8;

Nu = 1/3;

t = 1;

for ei=1:nelem        

node =[n(e(ei,1),:);n(e(ei,2),:);n(e(ei,3),:)]

[shape, pdn] =shapetri3(1/3, 1/3, 1/3);

[jm, invjm, pdc] =jacobimatrix(pdn, node);

ke = elemstiffpl(E, Nu,t, pdc, jm, 1);

K=assemblestiffpl(K,ke,e(ei,1),e(ei,2),e(ei,3));

end

fixedNodes=find(abs(n(:,1)+1)<1e-5);

Dof=[2*fixedNodes-1;2*fixedNodes];

force=zeros(nnode*2,1);

rightedge=find(abs(n(:,1)-1)<1e-5);

force(2*rightedge-1)=1000000;

activeDof=setdiff([1:nnode*2]',Dof);

U=K(activeDof,activeDof)\force(activeDof);

disp=zeros(nnode*2,1);

disp(activeDof)=U;

参考文献：

[1]    P.I. Kattan. Matlab有限元分析与应用.

[2]    王勖成. 有限单元法.

[3]    A.J.MFerreira. Matlab Codes for Finite Element Analysis.

liuweiji2311

不错:D:D

wanmiao

厉害，厉害

王鹏张猛

新人   学习下

快乐宙斯

受教了:)