필터 지우기
필터 지우기

Calculation stiffness matrix for 8-node hexahedral element.

조회 수: 34 (최근 30일)
Jan Tomaszewski
Jan Tomaszewski 2021년 4월 23일
댓글: Alehegn Tesfaye 2022년 5월 16일
Hi, I'm trying to solve stifness matrix for standard hexaheral element (100x100x100), I'm doing everything like in a literature but the final result from my solver does not match the result from the program (Midas NFX).The compared value is the nodal forces Qe . Can you help me find the bug?
clear
syms r s t % r,s,t - local coordinates
E=210000;% Modulus
v=0.28;%Poisson
G=E/(2*(1+v));%Shear modulus
coor_glob=[0 0 100
100 0 100
0 100 100
100 100 100
0 0 0
100 0 0
0 100 0
100 100 0];% global coordinates
x1= coor_glob(1,1); y1= coor_glob(1,2); z1= coor_glob(1,3);
x2= coor_glob(2,1); y2= coor_glob(2,2); z2= coor_glob(2,3);
x3= coor_glob(3,1); y3= coor_glob(3,2); z3= coor_glob(3,3);
x4= coor_glob(4,1); y4= coor_glob(4,2); z4= coor_glob(4,3);
x5= coor_glob(5,1); y5= coor_glob(5,2); z5= coor_glob(5,3);
x6= coor_glob(6,1); y6= coor_glob(6,2); z6= coor_glob(6,3);
x7= coor_glob(7,1); y7= coor_glob(7,2); z7= coor_glob(7,3);
x8= coor_glob(8,1); y8= coor_glob(8,2); z8= coor_glob(8,3);
S=[1/E -v/E -v/E 0 0 0;
-v/E 1/E -v/E 0 0 0;
-v/E -v/E 1/E 0 0 0;
0 0 0 1/(G) 0 0;
0 0 0 0 1/(G) 0 ;
0 0 0 0 0 1/(G)];
D=inv(S); %elasticity matrix
N1=1/8*((1-r)*(1-s)*(1+t));
N2=1/8*((1+r)*(1-s)*(1+t));
N3=1/8*((1-r)*(1+s)*(1+t));
N4=1/8*((1+r)*(1+s)*(1+t));
N5=1/8*((1-r)*(1-s)*(1-t));
N6=1/8*((1+r)*(1-s)*(1-t));
N7=1/8*((1-r)*(1+s)*(1-t));
N8=1/8*((1+r)*(1+s)*(1-t));% shape functions
x=N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8; %interpolate coord.
y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;
J=[diff(x,r) diff(y,r) diff(z,r);
diff(x,s) diff(y,s) diff(z,s);
diff(x,t) diff(y,t) diff(z,t)];%Jacobian matrix
d_N1=inv(J)*[diff(N1,r);diff(N1,s);diff(N1,t)];
d_N2=inv(J)*[diff(N2,r);diff(N2,s);diff(N2,t)];
d_N3=inv(J)*[diff(N3,r);diff(N3,s);diff(N3,t)];
d_N4=inv(J)*[diff(N4,r);diff(N4,s);diff(N4,t)];
d_N5=inv(J)*[diff(N5,r);diff(N5,s);diff(N5,t)];
d_N6=inv(J)*[diff(N6,r);diff(N6,s);diff(N6,t)];
d_N7=inv(J)*[diff(N7,r);diff(N7,s);diff(N7,t)];
d_N8=inv(J)*[diff(N8,r);diff(N8,s);diff(N8,t)];%differentials
B1=[d_N1(1) 0 0
0 d_N1(2) 0
0 0 d_N1(3)
d_N1(2) d_N1(1) 0
0 d_N1(3) d_N1(2)
d_N1(3) 0 d_N1(1)];%strain matrix for node 1
B2=[d_N2(1) 0 0
0 d_N2(2) 0
0 0 d_N2(3)
d_N2(2) d_N2(1) 0
0 d_N2(3) d_N2(2)
d_N2(3) 0 d_N2(1)];%strain matrix for node 2
B3=[d_N3(1) 0 0
0 d_N3(2) 0
0 0 d_N3(3)
d_N3(2) d_N3(1) 0
0 d_N3(3) d_N3(2)
d_N3(3) 0 d_N3(1)];%strain matrix for node 3
B4=[d_N4(1) 0 0
0 d_N4(2) 0
0 0 d_N4(3)
d_N4(2) d_N4(1) 0
0 d_N4(3) d_N4(2)
d_N4(3) 0 d_N4(1)];%strain matrix for node 4
B5=[d_N5(1) 0 0
0 d_N5(2) 0
0 0 d_N5(3)
d_N5(2) d_N5(1) 0
0 d_N5(3) d_N5(2)
d_N5(3) 0 d_N5(1)];%strain matrix for node 5
B6=[d_N6(1) 0 0
0 d_N6(2) 0
0 0 d_N6(3)
d_N6(2) d_N6(1) 0
0 d_N6(3) d_N6(2)
d_N6(3) 0 d_N6(1)];%strain matrix for node 6
B7=[d_N7(1) 0 0
0 d_N7(2) 0
0 0 d_N7(3)
d_N7(2) d_N7(1) 0
0 d_N7(3) d_N7(2)
d_N7(3) 0 d_N7(1)];%strain matrix for node 7
B8=[d_N8(1) 0 0
0 d_N8(2) 0
0 0 d_N8(3)
d_N8(2) d_N8(1) 0
0 d_N8(3) d_N8(2)
d_N8(3) 0 d_N8(1)];%strain matrix for node 8
B=[B1 B2 B3 B4 B5 B6 B7 B8];%main strin matrix for hexahedral
f(r,s,t)=B'*D*B*det(J);% function to integration
N=[-sqrt(1/3) sqrt(1/3)];% Gauss-Legendre quadrature points fot N=2
w=[1 1];% Gauss-Legendre quadrature weights fot N=2
wij=[w(1)*w(1)*w(2);
w(2)*w(1)*w(2);
w(1)*w(2)*w(2);
w(2)*w(2)*w(2);
w(1)*w(2)*w(2);
w(2)*w(1)*w(1);
w(1)*w(2)*w(1);
w(2)*w(2)*w(1)];%matrix of weights for future multiplikation
%% Gauss-Legendre integration:
k=f(N(1),N(1),N(2))*wij(1)+f(N(2),N(1),N(2))*wij(2)+f(N(1),N(2),N(2))*wij(3)+f(N(2),N(2),N(2))*wij(4)+f(N(1),N(1),N(1))*wij(5)+f(N(2),N(1),N(1))*wij(6)+f(N(1),N(2),N(1))*wij(7)+f(N(2),N(2),N(1))*wij(8);
K=eval(k); %elemennt stiffnes matrix
u=[0 0 0 0 0 0 0 0 0 -2.9665e-003 7.7690e-004 7.7690e-004 0 0 0 0 0 0 0 0 0 0 0 0]'; %translations from FEM program
Qe=K*u %nodal forces
%% Nodal forces from FEM program
Qe_fem=[4279.5049 2764.0085 -241.1585 -2782.678 -869.3088 -1551.515 4012.6929 1291.8174 1291.8174 -10000 0 0 3283.1487 675.8129 675.8129 -289.4951 -2069.6565 -2069.6565 4279.5049 -241.1585 2764.0085 -2782.678 -1551.515 -869.3088]';

답변 (1개)

Boris Kotsyubuk
Boris Kotsyubuk 2022년 3월 15일
It seems like your 4th row in your strain matrices (B1, B2, etc.) is incorrect. I believe the 4th row should be your last row in the matrix, while the 5th and 6th rows move up one row. I think the strain matrix for B1 should look like this:
B1 = [ d_N1(1), 0, 0;
0, d_N1(2), 0;
0, 0, d_N1(3);
0, d_N1(3), d_N1(2);
d_N1(3), 0, d_N1(1);
d_N1(2), d_N1(1), 0 ];
Hope this helps!
  댓글 수: 1
Alehegn Tesfaye
Alehegn Tesfaye 2022년 5월 16일
I did B matrix on octave exactly as you did but it still won't work, any help.
% Coordinates parameters
N=32;
x1=0; y1=2; z1=0;
x2=0; y2=0; z2=0;
x3=2; y3=0; z3=0;
x4=2; y4=2; z4=0;
x5=0; y5=2; z5=2;
x6=0; y6=0; z6=2;
x7=2; y7=0; z7=2;
x8=2; y8=2; z8=2;
% Element parameters
E=(170+N);
u=0.3;
% Plane Stress or Strain parameters
D=(E/((1+u)+(1-(2*u))))*[1-u u u 0 0 0;u 1-u u 0 0 0;u u 1-u 0 0 0;0 0 0 (1-(2*u))/2 0 0;0 0 0 0 (1-(2*u))/2 0; 0 0 0 0 0 (1-(2*u))/2];
% Load parameters
% Shape Function
syms s t w real;
N1=(1-s)*(1-t)*(1-w)/8;
N2=(1+s)*(1-t)*(1-w)/8;
N3=(1+s)*(1+t)*(1-w)/8;
N4=(1-s)*(1+t)*(1-w)/8;
N5=(1-s)*(1-t)*(1+w)/8;
N6=(1+s)*(1-t)*(1+w)/8;
N7=(1+s)*(1+t)*(1+w)/8;
N8=(1-s)*(1+t)*(1+w)/8;
% Coordinate Mapping
x=simplify(N1*x1+N2*x2+N3*x3+N4*x4+N5*x5+N6*x6+N7*x7+N8*x8);
y=simplify(N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8);
z=simplify(N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8);
% N Matrix
% Jacobian Matrix
syms J J0 real;
dxds=diff(x,s);
dxdt=diff(x,t);
dxdw=diff(x,w);
dyds=diff(y,s);
dydt=diff(y,t);
dydw=diff(y,w);
dzds=diff(z,s);
dzdt=diff(z,t);
dzdw=diff(z,w);
J=[dxds dyds dzds;dxdt dydt dzdt;dxdw dydw dzdw];
Jdet=simplify(det(J));
Jinv=simplify(inv(J));
% B Matrix
syms B real;
d_N1=Jinv*[diff(N1,s);diff(N1,t);diff(N1,w)];
d_N2=Jinv*[diff(N2,s);diff(N2,t);diff(N2,w)];
d_N3=Jinv*[diff(N3,s);diff(N3,t);diff(N3,w)];
d_N4=Jinv*[diff(N4,s);diff(N4,t);diff(N4,w)];
d_N5=Jinv*[diff(N5,s);diff(N5,t);diff(N5,w)];
d_N6=Jinv*[diff(N6,s);diff(N6,t);diff(N6,w)];
d_N7=Jinv*[diff(N7,s);diff(N7,t);diff(N7,w)];
d_N8=Jinv*[diff(N8,s);diff(N8,t);diff(N8,w)];
B1=[d_N1(1) 0 0
0 d_N1(2) 0
0 0 d_N1(3)
0 d_N1(3) d_N1(2)
d_N1(3) 0 d_N1(1)
d_N1(2) d_N1(1) 0];
B2=[d_N2(1) 0 0
0 d_N2(2) 0
0 0 d_N2(3)
0 d_N2(3) d_N2(2)
d_N2(3) 0 d_N2(1)
d_N2(2) d_N2(1) 0];
B3=[d_N3(1) 0 0
0 d_N3(2) 0
0 0 d_N3(3)
0 d_N3(3) d_N3(2)
d_N3(3) 0 d_N3(1)
d_N3(2) d_N3(1) 0];
B4=[d_N4(1) 0 0
0 d_N4(2) 0
0 0 d_N4(3)
0 d_N4(3) d_N4(2)
d_N4(3) 0 d_N4(1)
d_N4(2) d_N4(1) 0];
B5=[d_N5(1) 0 0
0 d_N5(2) 0
0 0 d_N5(3)
0 d_N5(3) d_N5(2)
d_N5(3) 0 d_N5(1)
d_N5(2) d_N5(1) 0];
B6=[d_N6(1) 0 0
0 d_N6(2) 0
0 0 d_N6(3)
0 d_N6(3) d_N6(2)
d_N6(3) 0 d_N6(1)
d_N6(2) d_N6(1) 0];
B7=[d_N7(1) 0 0
0 d_N7(2) 0
0 0 d_N7(3)
0 d_N7(3) d_N7(2)
d_N7(3) 0 d_N7(1)
d_N7(2) d_N7(1) 0];
B8=[d_N8(1) 0 0
0 d_N8(2) 0
0 0 d_N8(3)
0 d_N8(3) d_N8(2)
d_N8(3) 0 d_N8(1)
d_N8(2) d_N8(1) 0];
B=simplify([B1 B2 B3 B4 B5 B6 B7 B8]);
BT=B';
% Integration function
syms F real;
F=BT*D*B*Jdet;
% Stiffness Matrix
p=(1/sqrt(3));
o=function_handle(F,'vars',[s,t,w]);
K=(o(-p,-p,-p)+o(p,-p,-p)+o(p,p,-p)+o(-p,p,-p)+o(-p,-p,p)+o(p,-p,p)+o(p,p,p)+o(-p,p,p))

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Calculus에 대해 자세히 알아보기

제품


릴리스

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by