Calculating Tetrahedral Stiffness Matrix for FEM

26 次查看（过去 30 天）
Nathan Welch
Nathan Welch 2020-12-4
评论： Biswabhanu Puhan 2023-3-30
I'm writing some FEM simulations and have the mesh sorted.
I'm now putting the stiffness matrix together and feel that this is such a generic task that surely it's already been done - but after a quick google and look through the file exchange it appears not! There are a few things that mention trusses and things but nothing that takes an obvious input and gives a straightforward stiffness matrix.
Does anyone know of any good code that takes a mesh tetrahedralisation matrix T(Ne,4) and set of vertices V(Nv,3) and produces the linear stiffness matrix K(Nv,Nv)? Obviously higher orders would also be great but linear is enough to get started.
  2 个评论
Jason Nicholson
Jason Nicholson 2020-12-29
I am not aware of any automatic global stiffness matrix generator. I think you have to role your own. That is all I have. Maybe, someone who studies this will comment. However, most FEA code writer don't work in MATLAB except for high level interface. FEA is written in C and C++.
Nathan Welch
Nathan Welch 2020-12-29
Thanks - yeah this seems to be the standard method. Also quite a nice intro to FEM to work out how to do this I suppose!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Nathan Welch
Nathan Welch 2020-12-29
As there doesn't seem to be any obvious sources, I've created my own. It calculates the stiffness matrix for linear elements of a tetrahedral mesh. As with all online code - use at your own risk. If anyone can describe how to better calculate this, esp how to expand to 2nd order or (the dream) arbitrary order, I'd still be very much interested to hear it.
function K = LinearStiffnessMatrix3(nodes, elements)
%get sizes
Nn = size(nodes,1);
Ne = size(elements,1);
%create sparse matrix by storing (i,j) indices and the values
inds = 0;
Si = zeros(16*Ne,1); %ith indices of stiffness matrix
Sj = zeros(16*Ne,1); %jth indices of stiffness matrix
Sk = zeros(16*Ne,1); %value of (i,j) element of stiffness matrix
%4 indices vs 4 indices
ii = 1:4; ji = 1:4; [ii, ji] = ndgrid(ii,ji);
ii = ii(:); ji = ji(:);
% Loop through every element in the mesh - would be nice to do this without
% looping - or this maybe why it's best to use a compiled language
for i = 1:Ne
eleNodes = elements(i,1:4);
%uses simple integration of linear terms method
B = [nodes(eleNodes,:)'; [1, 1, 1, 1]];
C = inv(B);
ve = abs(det(B));
Kl2 = C(:,1:3)*C(:,1:3)'.*abs(ve)/6;
%store all 16 terms
inds = inds(end) + (1:16);
Si(inds) = eleNodes(ii);
Sj(inds) = eleNodes(ji);
Sk(inds) = Kl2(:);
end
%Build stiffness matrix - MATLABs sparse operator naturally adds together
%all terms with the same i and j - i.e. all contributions to the ith node
%from nodes j
K = sparse(Si,Sj,Sk,Nn,Nn);
end
  1 个评论
Biswabhanu Puhan
Biswabhanu Puhan 2023-3-30
Thank you Mr. Welch, did you find the formulation for 2nd order element (C3D10 in abaqus)? Kindly please let me know.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Modeling 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by