4

私は有限要素法について自分自身に教えようとしています。

私のコードはすべて、次のリンクページ 16-20 から改作されています http://homepages.cae.wisc.edu/~suresh/ME964Website/M964Notes/Notes/introfem.pdf

単一の 8 ノード キューブ要素に対して有限要素解析を実行するために、Matlab でプログラミングしています。ローカル軸 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

パラメトリック ドメイン ([0,1]^) から物理ドメインへのマッピング方法については、24 ページを確認してください。

于 2016-12-28T02:56:59.070 に答える
2

あなたが言ったように、シンボリックを使ってできると思いますが。Matlab での記号計算は非常に時間がかかると思います。

派生 N を手動で取得し、dN として保存し、必要なときに使用します。

よろしく、

ドイツ人

于 2017-06-16T18:25:10.513 に答える
0

形状関数を剛性マトリックスに代入する必要がある場合、24 の自由度があるため、剛性マトリックスは 24x24 である必要があります。解くには、線形システム (Ax=b) を構築する必要があります。右側は、解く PDE に基づいており、右側とソース項にノイマン境界条件を含める必要があります。2D 要素 (4 DOF) の Python では、次のようになります。

    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 に答える