清空筛选条件
清空筛选条件

Calculation stiffness matrix for 8-node hexahedral element.

45 次查看（过去 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 Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by