浅谈有限单元法中的形函数

易公子 2023年7月9日浏览：1916 评论：3 收藏：4

相信大家正在使用商软的同时，对于有限元的一些基础概念有些许淡忘，都在专注如何适用于复杂、高大上的场景，对于基础的理论多多少少不太关注，比如今天所要分享的有关 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