我对六面体单元的刚度矩阵计算有什么问题？我正在使用有限元方法。 Matlab 中的 MWE

Question

为了我自己的兴趣，我正在开发一个有限元法 (FEM) 代码来解决小问题。我已经有了 2D 元素、trias 和 quadriliterals 的工作代码。

但是我无法理解六面体单元的刚度矩阵。我确信我的形函数是正确的，而且导数也是正确的。

下面是我迄今为止在 Matlab 中取得的进展的一个最小工作示例。我用高斯积分。

它只有 1 个元素，一个 20×20×20 的立方体，在其中一个顶角施加 100e3 N 的点载荷，并固定在底部。这对应于 DOF 的 1-12 为零。 material 是线弹性的，E=210e3，泊松比 nu=0.3。

我认为不正确的地方：

我是否以错误的方式应用了边界条件？
应变位移矩阵B是不是不对？
我认为我不需要围绕刚度矩阵和 assemble 全局刚度矩阵进行调整，因为我只有 1 个元素。这不正确吗？

我的单元刚度计算逻辑是，在伪代码中：

对于每个高斯点：

   Find values for gauss points xhi, eta, my

   Find values for shape functions at current Gauss point

   Find values for shape function derivates wrt xhi, eta and my at current Gauss point

   Calculate Jacobian

   Calculate shape functions derivates wrt x,y and z

   Calculate strain displacement matrix B

   Find element stiffness matrix Ke

clc
clear all
% 8 gauss points in total
xhi = [ -0.57735,-0.57735,-0.57735,-0.57735, 0.57735, 0.57735,0.57735,0.57735];
eta = [ -0.57735,-0.57735, 0.57735,0.57735,-0.57735,-0.57735,0.57735,0.57735];
my = [ -0.57735, 0.57735,-0.57735,0.57735,-0.57735, 0.57735,-0.57735,0.57735];
ngp=8; % 2 gauss points in each direction
coord = [0. 0.  0.;
   20. 0.  0.;
   20. 20. 0.;
   0.  20. 0.;
   0.  0.  20.;
   20. 0.  20.;
   20. 20. 20.;
   0.  20. 20];

E=210e3;
poisson = 0.3;
Dcoeff = E/((1 + poisson) * (1 - 2 * poisson));
D = Dcoeff * [ (1 - poisson), poisson, poisson, 0, 0, 0;
   poisson, (1 - poisson), poisson, 0, 0, 0;
   poisson, poisson, (1 - poisson), 0, 0, 0;
   0, 0, 0, ((1 - 2 * poisson)/2), 0, 0;
   0, 0, 0, 0, ((1 - 2 * poisson)/2), 0;
   0, 0, 0, 0, 0, ((1 - 2 * poisson)/2)];
zero = zeros(1,8);
K_cpp = zeros(24,24);
% if this was a real problem i would have to loop over each element
% but i only have 1 element, so only need to loop over the Gauss points
for i=1:8
   %     // shape functions, 1/8=0.125
   N(1) = (1-xhi(i))*(1-eta(i))*(1-my(i))*0.125;
   N(2) = (1+xhi(i))*(1-eta(i))*(1-my(i))*0.125;
   N(3) = (1+xhi(i))*(1+eta(i))*(1-my(i))*0.125;
   N(4) = (1-xhi(i))*(1+eta(i))*(1-my(i))*0.125;
   N(5) = (1-xhi(i))*(1-eta(i))*(1+my(i))*0.125;
   N(6) = (1+xhi(i))*(1-eta(i))*(1+my(i))*0.125;
   N(7) = (1+xhi(i))*(1+eta(i))*(1+my(i))*0.125;
   N(8) = (1-xhi(i))*(1+eta(i))*(1+my(i))*0.125;
   %     // derive shape functions wrt xhi
   dNdXhi(1) = -(1-eta(i))*(1-my(i))*0.125;
   dNdXhi(2) =  (1-eta(i))*(1-my(i))*0.125;
   dNdXhi(3) =  (1+eta(i))*(1-my(i))*0.125;
   dNdXhi(4) = -(1+eta(i))*(1-my(i))*0.125;
   dNdXhi(5) = -(1-eta(i))*(1+my(i))*0.125;
   dNdXhi(6) =  (1-eta(i))*(1+my(i))*0.125;
   dNdXhi(7) =  (1+eta(i))*(1+my(i))*0.125;
   dNdXhi(8) = -(1+eta(i))*(1+my(i))*0.125;
   %     // derive shape functions wrt eta
   dNdEta(1) = -(1-xhi(i))*(1-my(i))*0.125;
   dNdEta(2) = -(1+xhi(i))*(1-my(i))*0.125;
   dNdEta(3) =  (1+xhi(i))*(1-my(i))*0.125;
   dNdEta(4) =  (1-xhi(i))*(1-my(i))*0.125;
   dNdEta(5) = -(1-xhi(i))*(1+my(i))*0.125;
   dNdEta(6) = -(1+xhi(i))*(1+my(i))*0.125;
   dNdEta(7) =  (1+xhi(i))*(1+my(i))*0.125;
   dNdEta(8) =  (1-xhi(i))*(1+my(i))*0.125;
   %     // derive shape functions wrt my
   dNdMy(1)  = -(1-xhi(i))*(1-eta(i))*0.125;
   dNdMy(2)  = -(1+xhi(i))*(1-eta(i))*0.125;
   dNdMy(3)  = -(1+xhi(i))*(1+eta(i))*0.125;
   dNdMy(4)  = -(1-xhi(i))*(1+eta(i))*0.125;
   dNdMy(5)  =  (1-xhi(i))*(1-eta(i))*0.125;
   dNdMy(6)  =  (1+xhi(i))*(1-eta(i))*0.125;
   dNdMy(7)  =  (1+xhi(i))*(1+eta(i))*0.125;
   dNdMy(8)  =  (1-xhi(i))*(1+eta(i))*0.125;
   %     // find Jacobian for each Gauss point
   dNdXhidEtadMy(1,:) = dNdXhi;
   dNdXhidEtadMy(2,:) = dNdEta;
   dNdXhidEtadMy(3,:) = dNdMy;
   %     J = [dNdXhi*coord(:,1),dNdXhi*coord(:,1),dNdXhi*coord(:,1);
   J = dNdXhidEtadMy*coord;
   
   %     // find shape functions derivate matrix wrt x, y & z
   dNdxdydz = inv(J)*dNdXhidEtadMy;
   %     // construct B matrix
   q_x = dNdxdydz(1,:);
   q_y = dNdxdydz(2,:);
   q_z = dNdxdydz(3,:);
   B_cpp = [  q_x, zero, zero;
       zero,  q_y, zero;
       zero, zero,  q_z;
       q_y,  q_x, zero;
       zero,  q_z,  q_y;
       q_z, zero, q_x];
   %     // compute Gauss point contribution to element Ke matrix
   K_cpp = K_cpp + (B_cpp'*D*B_cpp*det(J));
end
% nod 1-4 dof 1-3 = 0:
% apply boundary conditions on nodes 0-4 --> dofs 1-12
for i=1:12
   K_cpp(:,i)=zeros(1,24);
   K_cpp(i,:)=zeros(24,1);
   K_cpp(i,i)=1;
end
f=[zeros(1,23) 100e3]';
u=inv(K_cpp)*f;

这产生了解决方案 u:

u =
1.0e+14*[0  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0, -1.0293, -1.0531, -1.0511, -1.0503, -1.3835, -1.3835, -0.3304, -0.3304, -1.3791, -1.3895, -0.3304, -0.3304,

所以前12个DOF上规定的零位移边界条件似乎可行，但显然位移完全不正常。根据我仔细检查过的商业 FE 代码，它应该是 0.25mm 左右。

Answer 1

在 reddit 的好人的帮助下，我设法解决了这个问题。com/r/fea。

问题是应变-位移矩阵 B。

这是不正确的 B 矩阵，我有：

      [N1x N2x ... zeros(1,2*8)]
      [zeros(1,2*8)   N1y, N2y,... zeros(1,1*8)]
[B] = [zeros(1,2*8) N1z N2z ... ]
      [N1y N1x ... N1x N2x ... zeros(1,1*8)]
      [0 0 ... N1z N2z ... N1y N2y ...]
      [N1z N2z ... zeros(1,1*8) N1x N2x ...]

这是正确的定义。

      [N1x 0   0   N2x 0   0   ... ]
      [0   N1y 0   0   N2y 0   ... ]
[B] = [0   0   N1z 0   0   N2z ... ]
      [N1y N1x 0   N2y N2x 0   ... ]
      [0   N1z N1y 0   N2z N2y ... ]
      [N1z 0   N1x N2z 0   N2x ... ]

我对六面体单元的刚度矩阵计算有什么问题？我正在使用有限元方法。 Matlab 中的 MWE

Whas is wrong with my stiffness matrix calculation for a hexahedron element? I'm using finite element method. MWE in Matlab

matlab

finite-element-analysis