FEM-Course-Matlab/9.四面体单元悬臂梁Matlab有限元编程/Ke_C3D4.m

20 lines
869 B
Mathematica
Raw Permalink Normal View History

2024-01-28 16:46:36 +00:00
function [Ke]=Ke(D, ElementNodeCoordinate)
%<EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>ն<EFBFBD><EFBFBD><EFBFBD>
Ke=zeros(12,12);
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>NDerivative<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>JacobiDET<EFBFBD><EFBFBD>
[NDxyz, JacobiDET] = ShapeFunction( ElementNodeCoordinate);%[DN1Dx DN2Dx DN3Dx;DN1Dy DN2Dy DN3Dy;<EFBFBD><EFBFBD><EFBFBD><EFBFBD>]
Ve = JacobiDET/6;%
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B<EFBFBD><EFBFBD><EFBFBD><EFBFBD> <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>NDxyz<EFBFBD><EFBFBD><EFBFBD><EFBFBD>B<EFBFBD><EFBFBD><EFBFBD>м<EFBFBD><EFBFBD><EFBFBD>
B=zeros(6,12);
for i=1:4
sub=(i-1)*3+1:(i-1)*3+3;%<EFBFBD>Ӿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Χ
B(:,sub)=[NDxyz(1,i) 0 0;%NDx
0 NDxyz(2,i) 0;%NDy
0 0 NDxyz(3,i);%NDz
NDxyz(2,i) NDxyz(1,i) 0;
0 NDxyz(3,i) NDxyz(2,i);
NDxyz(3,i) 0 NDxyz(1,i)];
end
Ke=Ke+Ve*B'*D*B; %<EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD>֣<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
end