Метод конечных элементов программирования

Я пытаюсь научить себя, как насчет методов конечных элементов.

Весь мой код адаптирован со следующих страниц ссылки 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 ответа

Решение

Решил мою проблему с помощью анонимных функций и хранения полиномов в символьной матрице. пример:

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})
Источник

Вы должны проверить страницу 24 о том, как отобразить из параметрической области ([0,1]^) в физическую область.

Источник

Хотя я думаю, что вы можете сделать, как вы сказали, используя символическое. Я думаю, что символьные вычисления в Matlab очень трудоемки.

Я бы пошел для получения N вручную и сохранить как dN, и использовать его, когда это необходимо.

С Уважением,

Немецкий

Источник

после того, как у вас есть функции формы, вам нужно заменить их в матрице жесткости, матрица жесткости должна быть 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;
Источник
Другие вопросы по тегам