浅谈有限单元法中的形函数
最近在翻阅有限元相关的书籍，看到一本国内还不错的编程教材：《有限元法与 MATLAB——理论、体验与实践》。书中内容以理论结合编程实现，深入浅出，推荐大家学习，由于是 2022 年出版的新书，没有电子版，可在下方链接进行购买纸质版。
进入我们今天的主题：浅谈有限单元法中的形函数。相信大家正在使用商软的同时，对于有限元的一些基础概念有些许淡忘，都在专注如何适用于复杂、高大上的场景，对于基础的理论多多少少不太关注，比如今天所要分享的有关 FEM 中形函数的概念以及如何构造。
在基本的结构有限元编程中，大多是直接移值已有的形函数的形式，如四节点等参单元的形函数公式，从兴趣学习的角度来讲，搞明白形函数构造的方法或许比“直接拿来用”更有意义，
什么是形函数？
目前主流有限元分析力学问题时，位移 作为基本未知量，即有限元求解（位移）-几何方程（应变）--本构方程（应力）。以平面 3 节点三角形单元为例，分析单元的位移模式，引出形函数的概念。
单元位移模式可表示为：
单元内任一点的位移 可由节点的位移 通过形函数 进行内积求和得出。从数学的角度分析，形函数就是插值函数 。
形函数的性质
上一节引出了形函数的定义，那它有哪些性质（适用于任何单元类型）呢？（工科思维：讲完概念谈性质）
在 节点处值为1，在 、 节点处值为0，即满足 （Kronecker符号）：
在单元中任一点形函数之和等于1 为 、 的坐标函数，与位移函数有相同的阶次
构造单元形函数
本小节将简要介绍一下如何推导平面3节点三角形单元的形函数（理论＋代码实操），在后续的推文中，木木将分享从一维到三维各个常见单元类型的形函数推导方法。
广义坐标法
形函数 可表示为:
为三角形的面积，为保证 取正值，三个节点须按逆时针进行编号：
、 、 的计算式为：
以上通过节点坐标来推导单元形函数的方法被称为广义坐标法 。接下来，以书中的一个实例，基于Matlab符号语言使用广义坐标法推导三角形形函数。例1：一平面3节点三角形单元的节点坐标分别为： 、 、 ，求节点1、2、3的形函数 、 、 。
clear;
clc;
syms x y
xy = [0,0;1,0;0,1]
[A,abc] = SHP3(xy)
N1 = abc(1,1)+abc(1,2)*x+abc(1,3)*y;
N1 = N1/(2*A)
N2 = abc(2,1)+abc(2,2)*x+abc(2,3)*y;
N2 = N2/(2*A)
N3 = abc(3,1)+abc(3,2)*x+abc(3,3)*y;
N3 = N3/(2*A)
function [A,abc] = SHP3(xy)
% shape function for triangle element with 3 nodes
% xy(3,2)-- node coordinates of an element
A = [1, xy(1,1), xy(1,2);
1, xy(2,1), xy(2,2);
1, xy(3,1), xy(3,2)]
A = 0.5*det(A);
%-----------------------
a1 = [xy(2,1),xy(2,2);
xy(3,1),xy(3,2)];
a1 = det(a1);
a2 = [xy(3,1),xy(3,2);
xy(1,1),xy(1,2)];
a2 = det(a2);
a3 = [xy(1,1),xy(1,2);
xy(2,1),xy(2,2)];
a3 = det(a3);
%-------------------
b1 = [1,xy(2,2);
1,xy(3,2)];
b1 = -det(b1);
b2 = [1,xy(3,2);
1,xy(1,2)];
b2 = -det(b2);
b3 = [1,xy(1,2);
1,xy(2,2)];
b3 = -det(b3);
%--------------------
c1 = [1,xy(2,1);
1,xy(3,1)];
c1 = det(c1);
c2 = [1,xy(3,1);
1,xy(1,1)];
c2 = det(c2);
c3 = [1,xy(1,1);
1,xy(2,1)];
c3 = det(c3);
%---------------------
abc = [a1,b1,c1;
a2,b2,c2;
a3,b3,c3];
end
可以得出：
面积坐标法
根据形函数的基本性质，可由面积坐标 进行表示：
其中， 、 、 为下图所示三角形内任意一点 的面积坐标，且
其中，
继续上一小节的例题，使用面积坐标可进行如下推导：
clear;
clc;
syms x y
assume (x, 'real')
assume (y, 'real')
x1=0;
y1=0;
x2=1;
y2=0;
x3=0;
y3=1;
A = [1,x1,y1;
1,x2,y2;
1,x3,y3];
A = 1/2*det(A)
A1 = [1,x,y;
1,x2,y2;
1,x3,y3];
A1 = 1/2*det(A1)
A2 = [1,x1,y1;
1,x,y;
1,x3,y3];
A2 = 1/2*det(A2)
A3 = [1,x1,y1;
1,x2,y2;
1,x,y];
A3 = 1/2*det(A3)
N1 = A1/A
N2 = A2/A
N3 = A3/A
可以得出：
形函数云图绘制
绘制以上形函数 在单元内的云图，代码由《有限元法与 MATLAB——理论、体验与实践》提供：
clear;clc;
p1 = [0,0];
p2 = [1,0];
p3 = [0,1];
xy31 = Pdivide(p3,p1,5)
xy32 = Pdivide(p3,p2,5)
xy1 = xy31(1,:)
xy2 = Pdivide(xy31(2,:),xy32(2,:),2)
xy3 = Pdivide(xy31(3,:),xy32(3,:),3)
xy4 = Pdivide(xy31(4,:),xy32(4,:),4)
xy5 = Pdivide(xy31(5,:),xy32(5,:),5)
xy = [xy1;xy2;xy3;xy4;xy5]
N1 = 1 - xy(:,1) - xy(:,2)
Enod = delaunay(xy(:,1),xy(:,2))
Nxy = [(1:15)',xy]
Enod = [(1:16)',Enod]
N1 = [(1:15)',N1]
pltvc2d3n_x(Enod,Nxy,N1)
function xy = Pdivide(P1,P2,n)
% Point partition function
% P1(1,2) - [x1,y1]
% P2(1,2) - [x2,y2]
% xy(n,2)
x = linspace(P1(1),P2(1),n);
y = linspace(P1(2),P2(2),n);
xy = [x',y']
end
function pltvc2d3n_x(Enod,Nxy,Nvar)
figure
hold on
axis equal
axis off
m = size(Enod,1)
for i=1:m
k = Enod(i,2:4)
x = Nxy(k,2)
y = Nxy(k,3)
c = Nvar(k,2) %
h = fill(x,y,c)
set(h, 'linestyle','none')
end
colorbar('location','eastoutside')
end
由云图可知， 在节点1处的值为1，在别的节点处值为0。
以上内容就是木木本次分享有关有限元中有关形函数的一些浅显理解，由于本人知识面较窄，对理论的理解未必都是正确的，欢迎大神、同行批评指正！
