这篇文章会将FEM全流程走一遍，包括网格、矩阵组装、求解、后处理。内容是大三时的大作业，今天拿出来回顾下。

1. 问题简介

涡轮机叶片需要冷却以提高涡轮的性能和涡轮叶片的寿命。我们现在考虑一个如上图所示的叶片，叶片处在一个高温环境中，中间通有四个冷却孔。

(资料图)

假设为稳态，那么叶片内导热微分方程为：

内部区域： （扩散方程）

边界：

（外表面）

（内部冷却孔）

2.模型2.1几何模型

我们简化为二维模型，如下图所示：

点坐标：

1:0.0,0.0 6:597.6,45.9 11:344.7,50.0

2:20.9,28.8 7:870.0,0.0 12:435.8,44.5

3:117.4,62.9 8:85.0,40.0 13:521.2,37.0

4:240.4,69.6 9:174.5,49.4 14:605.0,30.0

5:417.5,62.4 10:260.0,50.0 15:694.7,22.2

2.2 单位系统和物性

长度单位：mm

温度单位： K

功率单位：W

k=14.7*10^-3 W/mm. K

h_ext=205.8*10^-6 W/m².K

h_int=65.8*10^-6 W/m².K

注意：在后面的矩阵组装中的h_wall = h / k

2.3网格

用开源软件Gmsh生成网格。

首先写geo文件

注意要把外表面和中间空洞用Physical Line定义

lc = 10;  Point(1) = {0, 0, 0, lc};Point(2) = {20.9,28.8,0, lc};Point(3) = {117.4,62.9,0, lc};Point(4) = {240.4,69.6,0, lc};Point(5) =  {417.5,62.4,0, lc};Point(6) = {597.6,45.9,0, lc};Point(7) = {870.0,0.0, 0,lc};Point(8) = {85.0,40.0, 0,lc};Point(9) = {174.5,49.4,0, lc};Point(10) = {260.0,50.0,0, lc};Point(11) = {344.7,50.0,0, lc};Point(12) = {435.8,44.5,0, lc};Point(13) = {521.2,37.0,0, lc};Point(14) = {605.0,30.0,0, lc};Point(15) = {694.7,22.2,0, lc};//+Spline(1) = {1, 2, 3, 4, 5, 6, 7};//+Symmetry {0, 1, 0, 0} {  Duplicata { Point{1}; Point{2}; Point{3}; Point{4}; Point{5}; Line{1}; Point{6}; Point{7}; Point{8}; Point{9}; Point{10}; Point{11}; Point{12}; Point{13}; Point{14}; Point{15}; }}//+Line Loop(1) = {1, -2};//+Line(3) = {8, 9};//+Line(4) = {9, 33};//+Line(5) = {33, 32};//+Line(6) = {32, 8};//+Line(7) = {10, 11};//+Line(8) = {11, 35};//+Line(9) = {35, 34};//+Line(10) = {34, 10};//+Line(11) = {12, 13};//+Line(12) = {13, 37};//+Line(13) = {37, 36};//+Line(14) = {36, 12};//+Line(15) = {14, 15};//+Line(16) = {15, 39};//+Line(17) = {39, 38};//+Line(18) = {38, 14};//+Line Loop(2) = {3, 4, 5, 6};//+Line Loop(3) = {7, 8, 9, 10};//+Line Loop(4) = {11, 12, 13, 14};//+Line Loop(5) = {15, 16, 17, 18};//+Physical Line(0)={1, -2};Physical Line(1)={3, 4, 5, 6};Physical Line(2)={7, 8, 9, 10};Physical Line(3)={11, 12, 13, 14};Physical Line(4)={15, 16, 17, 18};Plane Surface(1) = {1, 2, 3, 4, 5};Physical Surface(1) = {1};

边界信息如何保存？

边界edges需要用标记，

在matlab程序中用bedge(3,Nbc)存储边界信息，前两个数字代表边的两端节点编号，第三个数字代表属于哪一个边。

生成网格后导出为“blademesh.m”用以后续使用，注意不要勾选Save all elements，否则会没有边界信息。

我用gmsh-4.4.1-Windows64版本，可以导出边界信息。但是新版的gmsh导出为.m文件时，边界信息无法保存。

3. 矩阵组装3.1 控制方程3.2 系统矩阵

其中的Ω代表全域，我们将全域分解为一个个单元，这就是有限元的思想。

计算每个单元（Ωe）的刚度矩阵，然后每项加到整体刚度矩阵：

4. 代码实现（matlab）

步骤	工具或函数
定义求解域并生成网格	gmsh导出网格为blademesh.m
读入网格信息，数据转换	bladeread
矩阵和矢量组装，线性方程组求解	bleadheat
查看结果	bladeplot

主程序：bladeheat.m

% Clear variablesclear all;% Set gas temperature and wall heat transfer coefficients at% boundaries of the blade.  Note: Tcool(i) and hwall(i) are the% values of Tcool and hwall for the ith boundary which are numbered% as follows:  %%   1 = external boundary (airfoil surface)%   2 = 1st internal cooling passage (from leading edge)%   3 = 2nd internal cooling passage (from leading edge)%   3 = 3rd internal cooling passage (from leading edge)%   3 = 4th internal cooling passage (from leading edge)% Tcool = [1300, 200, 200, 200, 200];% hwall = [14, 4.7, 4.7, 4.7, 4.7];Tcool = [1573, 473, 473, 473, 473];h = [205.8*10^-6, 65.8*10^-6, 65.8*10^-6, 65.8*10^-6, 65.8*10^-6];k = 14.7*10^-3;hwall = h / k;% Load in the grid file% NOTE:  after loading a gridfile using the load(fname) command,%        three important grid variables and data arrays exist.  These are:%% Nt: Number of triangles (i.e. elements) in mesh%% Nv: Number of nodes (i.e. vertices) in mesh%% Nbc: Number of edges which lie on a boundary of the computational%      domain.% % tri2nod(3,Nt):  list of the 3 node numbers which form the current%                 triangle.  Thus, tri2nod(1,i) is the 1st node of%                 the i"th triangle, tri2nod(2,i) is the 2nd node%                 of the i"th triangle, etc.%% xy(2,Nv): list of the x and y locations of each node.  Thus,%           xy(1,i) is the x-location of the i"th node, xy(2,i)%           is the y-location of the i"th node, etc.%% bedge(3,Nbc): For each boundary edge, bedge(1,i) and bedge(2,i) %               are the node numbers for the nodes at the end%               points of the i"th boundary edge.  bedge(3,i) is an%               integer which identifies which boundary the edge is%               on. In this solver, the third value has the%               following meaning:%%               bedge(3,i) = 0: edge is on the airfoil%               bedge(3,i) = 1: edge is on the first cooling passage%               bedge(3,i) = 2: edge is on the second cooling passage%               bedge(3,i) = 3: edge is on the third cooling passage%               bedge(3,i) = 4: edge is on the fourth cooling passage% bladeread;% Start timerTime0 = cputime;% Zero stiffness matrixK = zeros(Nv, Nv);b = zeros(Nv, 1);% Zero maximum element sizehmax = 0;% Loop over elements and calculate residual and stiffness matrixfor ii = 1:Nt,    kn(1) = tri2nod(1,ii);  kn(2) = tri2nod(2,ii);  kn(3) = tri2nod(3,ii);      xe(1) = xy(1,kn(1));  xe(2) = xy(1,kn(2));  xe(3) = xy(1,kn(3));  ye(1) = xy(2,kn(1));  ye(2) = xy(2,kn(2));  ye(3) = xy(2,kn(3));  % Calculate circumcircle radius for the element  % First, find the center of the circle by intersecting the median  % segments from two of the triangle edges.    dx21 = xe(2) - xe(1);  dy21 = ye(2) - ye(1);  dx31 = xe(3) - xe(1);  dy31 = ye(3) - ye(1);  x21  = 0.5*(xe(2) + xe(1));  y21  = 0.5*(ye(2) + ye(1));  x31  = 0.5*(xe(3) + xe(1));  y31  = 0.5*(ye(3) + ye(1));  b21 = x21*dx21 + y21*dy21;  b31 = x31*dx31 + y31*dy31;  xydet = dx21*dy31 - dy21*dx31;    x0 = (dy31*b21 - dy21*b31)/xydet;  y0 = (dx21*b31 - dx31*b21)/xydet;    Rlocal = sqrt((xe(1)-x0)^2 + (ye(1)-y0)^2);  if (hmax < Rlocal),    hmax = Rlocal;  end    % Calculate all of the necessary shape function derivatives, the  % Jacobian of the element, etc.    % Derivatives of node 1"s interpolant   dNdxi(1,1) = -1.0; % with respect to xi1  dNdxi(1,2) = -1.0; % with respect to xi2    % Derivatives of node 2"s interpolant  dNdxi(2,1) =  1.0; % with respect to xi1  dNdxi(2,2) =  0.0; % with respect to xi2  % Derivatives of node 3"s interpolant  dNdxi(3,1) =  0.0; % with respect to xi1  dNdxi(3,2) =  1.0; % with respect to xi2    % Sum these to find dxdxi (note: these are constant within an element)  dxdxi = zeros(2,2);  for nn = 1:3,    dxdxi(1,:) = dxdxi(1,:) + xe(nn)*dNdxi(nn,:);    dxdxi(2,:) = dxdxi(2,:) + ye(nn)*dNdxi(nn,:);  end    % Calculate determinant for area weighting  J = dxdxi(1,1)*dxdxi(2,2) - dxdxi(1,2)*dxdxi(2,1);  A = 0.5*abs(J); % Area is half of the Jacobian    % Invert dxdxi to find dxidx using inversion rule for a 2x2 matrix  dxidx = [ dxdxi(2,2)/J, -dxdxi(1,2)/J; ...       -dxdxi(2,1)/J,  dxdxi(1,1)/J];    % Calculate dNdx   dNdx = dNdxi*dxidx;  % Add contributions to stiffness matrix for node 1 weighted residual  K(kn(1), kn(1)) = K(kn(1), kn(1)) + (dNdx(1,1)*dNdx(1,1) + dNdx(1,2)*dNdx(1,2))*A;  K(kn(1), kn(2)) = K(kn(1), kn(2)) + (dNdx(1,1)*dNdx(2,1) + dNdx(1,2)*dNdx(2,2))*A;  K(kn(1), kn(3)) = K(kn(1), kn(3)) + (dNdx(1,1)*dNdx(3,1) + dNdx(1,2)*dNdx(3,2))*A;    % Add contributions to stiffness matrix for node 2 weighted residual  K(kn(2), kn(1)) = K(kn(2), kn(1)) + (dNdx(2,1)*dNdx(1,1) + dNdx(2,2)*dNdx(1,2))*A;  K(kn(2), kn(2)) = K(kn(2), kn(2)) + (dNdx(2,1)*dNdx(2,1) + dNdx(2,2)*dNdx(2,2))*A;  K(kn(2), kn(3)) = K(kn(2), kn(3)) + (dNdx(2,1)*dNdx(3,1) + dNdx(2,2)*dNdx(3,2))*A;    % Add contributions to stiffness matrix for node 3 weighted residual  K(kn(3), kn(1)) = K(kn(3), kn(1)) + (dNdx(3,1)*dNdx(1,1) + dNdx(3,2)*dNdx(1,2))*A;  K(kn(3), kn(2)) = K(kn(3), kn(2)) + (dNdx(3,1)*dNdx(2,1) + dNdx(3,2)*dNdx(2,2))*A;  K(kn(3), kn(3)) = K(kn(3), kn(3)) + (dNdx(3,1)*dNdx(3,1) + dNdx(3,2)*dNdx(3,2))*A;  end% Loop over boundary edges and account for bc"s% Note: the bc"s are all convective heat transfer coefficient bc"s% so the are of "Robin" form.  This requires modification of the% stiffness matrix as well as impacting the right-hand side, b.%for ii = 1:Nbc,  % Get node numbers on edge  kn(1) = bedge(1,ii);  kn(2) = bedge(2,ii);    % Get node coordinates  xe(1) = xy(1,kn(1));  xe(2) = xy(1,kn(2));    ye(1) = xy(2,kn(1));  ye(2) = xy(2,kn(2));  % Calculate edge length  ds = sqrt((xe(1)-xe(2))^2 + (ye(1)-ye(2))^2);    % Determine the boundary number  bnum = bedge(3,ii) + 1;  % Based on boundary number, set heat transfer bc  K(kn(1), kn(1)) = K(kn(1), kn(1)) + hwall(bnum)*ds*(1/3);  K(kn(1), kn(2)) = K(kn(1), kn(2)) + hwall(bnum)*ds*(1/6);  b(kn(1))        = b(kn(1))        + hwall(bnum)*ds*0.5*Tcool(bnum);    K(kn(2), kn(1)) = K(kn(2), kn(1)) + hwall(bnum)*ds*(1/6);  K(kn(2), kn(2)) = K(kn(2), kn(2)) + hwall(bnum)*ds*(1/3);  b(kn(2))        = b(kn(2))        + hwall(bnum)*ds*0.5*Tcool(bnum);    end% Solve for temperatureTsol = K\b;% Finish timerTime1 = cputime;% Plot solutionbladeplot;% Report outputsTmax = max(Tsol);Tmin = min(Tsol);fprintf("Number of nodes      = %i\n",Nv);fprintf("Number of elements   = %i\n",Nt);fprintf("Maximum element size = %5.3f\n",hmax);fprintf("Minimum temperature  = %6.1f\n",Tmin);fprintf("Maximum temperature  = %6.1f\n",Tmax);fprintf("CPU Time (secs)      = %f\n",Time1 - Time0); 

bladeread.m

% Read three important grid variables and data arrays% Nt: Number of triangles (i.e. elements) in mesh% Nv: Number of nodes (i.e. vertices) in mesh% Nbc: Number of edges which lie on a boundary of the computational%      domain.% tri2nod(3,Nt):  list of the 3 node numbers which form the current%                 triangle.  Thus, tri2nod(1,i) is the 1st node of%                 the i"th triangle, tri2nod(2,i) is the 2nd node%                 of the i"th triangle, etc.% xy(2,Nv): list of the x and y locations of each node.  Thus,%           xy(1,i) is the x-location of the i"th node, xy(2,i)%           is the y-location of the i"th node, etc.% bedge(3,Nbc): For each boundary edge, bedge(1,i) and bedge(2,i) %               are the node numbers for the nodes at the end%               points of the i"th boundary edge.  bedge(3,i) is an%               integer which identifies which boundary the edge is%               on. In this solver, the third value has the%               following meaning:%%               bedge(3,i) = 0: edge is on the airfoil%               bedge(3,i) = 1: edge is on the first cooling passage%               bedge(3,i) = 2: edge is on the second cooling passage%               bedge(3,i) = 3: edge is on the third cooling passage%               bedge(3,i) = 4: edge is on the fourth cooling passage% clcrun("blademesh.m");Nv=msh.nbNod;Nt=size(msh.TRIANGLES,1);Nbc=size(msh.LINES,1);for i=1:Nt    tri2nod(1,i)=msh.TRIANGLES(i,1);    tri2nod(2,i)=msh.TRIANGLES(i,2);    tri2nod(3,i)=msh.TRIANGLES(i,3);endfor i=1:Nv    xy(1,i)=msh.POS(i,1);    xy(2,i)=msh.POS(i,2);endfor i=1:Nbc    bedge(1,i)=msh.LINES(i,1);    bedge(2,i)=msh.LINES(i,2);    bedge(3,i)=msh.LINES(i,3);end

bladeplot.m

% Plot T in trianglesfigure;for ii = 1:Nt,  for nn = 1:3,    xtri(nn,ii) = xy(1,tri2nod(nn,ii));    ytri(nn,ii) = xy(2,tri2nod(nn,ii));    Ttri(nn,ii) = Tsol(tri2nod(nn,ii));  endendHT = patch(xtri,ytri,Ttri);axis("equal");set(HT,"LineStyle","none");title("Temperature (K)");% caxis([298,1573]);colormap(jet);HC = colorbar;hold on; bladeplotgrid; hold off;

5. 计算结果