好久没有更新博客了~前一段又有点对人生迷茫,玩了一段三国志13,fm,很是懈怠,五一前本来想在把读书时的代码找出来,转换成python代码的任务迟迟没有完成。今天,迎来了三年股市最大跌幅,三月清仓,月定投华能国际,竟然在四月底跑赢了指数,定投的指数基金也在上月成功减仓!准备月定投中国中铁和中国石油,现在这里mark一下!
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
回归正题,今天要用Python实现了任意四边形单元的面积计算,需要理解的几个基本概念:形函数,等参单元、jacobian矩阵和高斯积分。
形函数:这个概念的提出前提是单元内所感兴趣的函数是连续的,因此可以建立一组插值函数表达式,方便获得单元内某个坐标所对应的函数值。
形函数和等餐单元如图所示:
等参单元:它是作为唯一一个标准的四边形单元存在的,其他形状的四边形都是由它乘以Jacobian矩阵变换而来的,打个不恰当的比喻,等参单元就好比是父亲,他所有的孩子都是父亲基因依照某种变换法则得到的。
Jacobian矩阵:它表征的就一个任意形状的四边形单元和等参单元之间的转换关系(也就是上面所说的变换法则),包括各个方向的拉伸和旋转,其行列式,则是该单元与等参单元面积的比值。
高斯积分:前提要求,四边形内的函数必须是一个连续函数,而采用高斯积分的目的,是将连续积分转换成离散积分,便于计算机编程实现,因此高斯积分点越多,其精度越高,当然,精度高会带来计算时间的增加。如果一个单元内要积分的函数非线性程度不是很强的话,可以考虑采用较少的高斯积分点。进一步可以得出,如果单元尺寸足够小,其内部的非线性将会进一步降低,因此可以考虑选择较少的高斯积分点进行后续的积分计算。参见下图。
# -*- coding: utf-8 -*-
"""
Created on Sat Apr 27 07:05:44 2019
@author: Yujin.Wang
"""
import numpy as np
import pandas as pd
def funcQuadShapeDiff(eps,eta):
'''形函数导数'''
NDiff = []
NDiff = [[-0.25*(1-eta),-0.25*(1+eta),0.25*(1+eta),0.25*(1-eta)],
[-0.25*(1-eps),0.25*(1-eps),0.25*(1-eps),-0.25*(1-eps)]]
return np.matrix(NDiff)
def funcQuadJacobian(QuadNode):
''''Jacobian矩阵计算和高斯积分'''
a = pow(1./3,1./2)
GaussPoint = [[-a,a],[a,a],[a,-a],[-a,-a]] #等参单元的二阶高斯积分点坐标
GaussSum = 0
for i in range(4):
J = np.matrix([[0.,0.],[0.,0.]])
NDiff = funcQuadShapeDiff(GaussPoint[i][0],GaussPoint[i][1])
J[0,0] = np.matmul(NDiff[0],QuadNode[:,1])
J[0,1] = np.matmul(NDiff[0],QuadNode[:,2])
J[1,0] = np.matmul(NDiff[1],QuadNode[:,1])
J[1,1] = np.matmul(NDiff[1],QuadNode[:,2])
GaussSum +=np.linalg.det(J)
return GaussSum
if __name__ == "__main__":
QuadNode = [[1,1.,-0.],[2,1.,3.],[3,-1.,3.],[4,-1.,0.]]#节点编号,X坐标,Y坐标
QuadNode = np.matrix(QuadNode)#转变成矩阵形式
a = funcQuadJacobian(QuadNode)#计算单元面积
print "Area:",a
对应的文献可参加I.M.Smith的有限元方法编程。
参考文献:
I.M,Smith 著,张新春等 译. 有限单元方法编程(第五版)