我正在努力自学有限元方法。
我的所有代码都改编自以下链接第16-20页http://homepages.cae.wisc.edu/~suresh/ME964网站/M64Notes/Notes/introfme.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 ]
这是在CCD_ 4上进行的运算符。(B=DN
)
稍后,正如本文所示,我将进行计算,涉及[B]
矩阵在该元素体积上的积分。
所以,我的问题是,如何将这些多项式形状函数存储在矩阵中,用微分运算,然后用数值积分。根据我现在设置的方式,我可以看出它不起作用,因为我已经将函数定义为区间[0,1]
上的向量,然后将这些向量存储在[N]
矩阵中。然后利用CCD_ 9函数进行适当的微分,得到CCD_。但由于[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})
您应该查看第24页,了解如何从参数域([0,1]^)映射到物理域。
尽管我认为您可以按照您所说的那样使用符号。我认为Matlab中的符号计算非常耗时。
我会手动获取导数N并存储为dN,并在需要时使用它。
问候,
德国
在获得需要在刚度矩阵中替换的形状函数后,刚度矩阵应为24x24,因为您有24个自由度。为了求解,你需要建立一个线性系统(Ax=b),右手边是基于你正在求解的PDE,你必须在右手边加上源项中包括neuman边界条件。在python中,二维元素(4自由度)将类似于:
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;