4

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

我的所有代码均改编自以下链接页面 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]我正在阅读的文字，矩阵将像这样排列：

%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]我认为这会导致问题。你们将如何进行我上面发布的教科书中描述的这些计算？

4

4 回答 4

3

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

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})
于 2016-01-09T03:52:18.503 回答
3

您应该查看第 24 页，了解如何从参数域 ([0,1]^) 映射到物理域。

于 2016-12-28T02:56:59.070 回答
2

虽然我认为你可以按照你说的做，使用符号。我认为 Matlab 中的符号计算非常耗时。

我会手动导出 N 并将其存储为 dN，并在需要时使用它。

问候，

德语

于 2017-06-16T18:25:10.513 回答
0

获得形状函数后，您需要将其替换为刚度矩阵，刚度矩阵应为 24x24，因为您有 24 个自由度。要求解，您需要构建一个线性系统 (Ax=b)，右侧基于您正在求解的 PDE，并且您必须在右侧包括纽曼边界条件加上源项。在 python 中，用于 2d 元素（4 DOF）将如下所示：

    def shapefxncoef (Valxy):
    #creating a temporary metrix to store zeros and get the size of the shape
    #function matrix.
    n_temp = np.zeros((4,4))
    #filling the values of the matrix with a loop.
    for i in range(4):
        #the values used in the matrix are from the Valxy x and y components.
        xi = Valxy [0, i];
        yi = Valxy [1, i];
        n_temp[i, 0] = 1;
        n_temp[i, 1] = xi;
        n_temp[i, 2] = yi;
        n_temp[i, 3] = xi*yi;
        #this gives an identity matrix and the stiffness matric can be derived 
        #if we take the inverse.
    n = np.linalg.inv(n_temp);
    return n;

def N (Valxy, x, y):
    n = shapefxncoef (Valxy);
    res = n[0, :] + n[1, :]*x + n[2, :]*y + n[3, :]*x*y;
    return res;

def Be (Valxy, x, y):
    res = np.zeros ((2,4));
    res_temp = shapefxncoef (Valxy);
    for i in range (4):
        res_tempi = res_temp[:, i];
        dNix = res_tempi[1] + res_tempi[3]*y;
        dNiy = res_tempi[2] + res_tempi[3]*x;
        res[0, i] = dNix;
        res[1, i] = dNiy;
    return res;

def Ke (Valxy, conduct):
    a = lambda x, y: conduct * np.dot ((Be(Valxy, x, y)).T, Be(Valxy, x, y));
    k = intr.integrateOnQuadrangle (Valxy.T, a, np.zeros((4,4)));
    return k;
于 2021-03-03T17:40:57.763 回答