matlab - 编程有限元方法

我正在尝试自学有限元方法。

我的所有代码均改编自以下链接页面 16-20
http://homepages.cae.wisc.edu/~suresh/ME964Website/M964Notes/Notes/introfem.pdf

我正在 Matlab 中编程以对单个 8 节点立方体元素执行有限元分析。我已经定义了 xi,eta,zeta 局部轴(我们现在可以将其视为 x、y、z)，因此我得到以下形状函数:

%%shape functions

zeta = 0:.01:1;
eta = 0:.01:1;
xi = 0:.01:1;

N1 = 1/8*(1-xi).*(1-eta).*(1-zeta);
N2 = 1/8*(1+xi).*(1-eta).*(1-zeta);
N3 = 1/8*(1+xi).*(1+eta).*(1-zeta);
N4 = 1/8*(1-xi).*(1+eta).*(1-zeta);
N5 = 1/8*(1-xi).*(1-eta).*(1+zeta);
N6 = 1/8*(1+xi).*(1-eta).*(1+zeta);
N7 = 1/8*(1+xi).*(1+eta).*(1+zeta);
N8 = 1/8*(1-xi).*(1+eta).*(1+zeta);

根据我正在阅读的文本，[N] Matrix 将像这样排列:

%N Matrix
N= [N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8 0 0;
    0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8 0;
    0 0 N1 0 0 N2 0 0 N3 0 0 N4 0 0 N5 0 0 N6 0 0 N7 0 0 N8];

要找到 [B] 矩阵，我必须使用以下 [D] 矩阵:

%%Del Matrix for node i
%[  d/dx   0     0
%    0    d/dy   0
%    0     0    d/dz        . . .
%   d/dy  d/dx   0
%    0    d/dz  d/dy
%   d/dz   0    d/dx   ]

这是继续 [N] 的运算符。 ( B=DN )

稍后，如文本所示，我将进行涉及该 [B] 矩阵对该元素体积的积分的计算。

所以，我的问题是，如何将这些多项式形状函数存储在矩阵中，对它们进行微分运算，然后对它们进行数值积分。我可以通过我现在设置的方式判断它不会工作，因为我已经将函数定义为一个区间 [0,1] 上的向量，然后将这些向量存储在 [N] 矩阵中。然后使用 diff() 函数进行适当的微分以找到 [B] 矩阵。
但是由于 [B] 的矩阵元素现在是区间 [0,1] 上的向量，我认为这会导致问题。你们会如何进行我上面发布的教科书中描述的这些计算？

最佳答案

使用匿名函数并将多项式存储在符号矩阵中解决了我的问题。例子:

syms xi eta zeta
N1= ... %type out in terms of xi eta and zeta
.
.
.
dN1dXi = diff(N1,xi) %symbolic differentiation with respect to xi

还可以在需要时执行符号积分:

intN1 = int(N1,xi,lowerLimit,upperLimit) %symbolic integration with respect to xi

当准备好替换实际值以评估符号函数时:

subs(N1,{xi,eta,zeta},{value1,value2,value3})