一种针对椭圆区域的高质量网格划分方法

在处理一些积分域为椭圆区域的数值积分问题时，通常需要将积分域进行离散处理。本文介绍了一种针对椭圆区域的网格离散方法，该方法基于椭圆坐标系变换实现，无需借助网格划分软件，即可生成高质量网格。

这里首先介绍一种比较常规的处理方法，该方法基于椭圆的仿射变换实现。仿射变换，又称仿射映射，是指在几何中，一个向量空间进行一次线性变换并接上一个平移，变换为另一个向量空间。

对于椭圆

其中a和b分别为椭圆的长轴和短轴。

如果我们定义变换

则可以将椭圆变换为平面内的单位圆：

以上变换被称为圆锥曲线标准变换。在经过变换后，点仍然是点，线仍然是线，在直线上的点仍然在直线上。

这提示我们可以首先构造如下所示的圆形区域网格，对于圆形区域可以较容易获得非常规则的计算网格，如下图所示。

用于变换的网格（半径r=1）

基于上述变换，对网格节点坐标进行坐标变换，即可以得到如下所示的椭圆形网格。

变换之后的椭圆区域网格（a=2，b=1.2）

上述坐标映射过程实际上就是对圆形进行了简单的伸缩变换，可以看到变换之后网格的正交性将不再得到满足。因此如果需要生成具有高长宽比的椭圆，则很容易导致变换之后出现网格畸形的情况。下面，本文介绍一种新的变换方式，使得变换之后的网格仍然能满足正交性。

首先，我们引入一个椭圆坐标系(ξ, η)，从椭圆坐标系到笛卡尔坐标系的变换可表示为：

其中c为椭圆的焦点，可表示为：

下图给出了椭圆坐标系的示意图，其中红色实线代表坐标η，蓝色虚线表示ξ。

椭圆坐标系统的坐标线（a=2，b=1.2）

因此，我们可以首先构造一个矩形区域并划分网格，用该矩形区域的节点坐标来表示ξ和η，然后应用上述坐标变换得到椭圆区域的网格。当然，这里需要确定矩形的尺寸，也就是ξ和η的坐标范围。显然η的取值范围为[0, 2π]，因此下面需要确定ξ的取值范围。为此我们可以首先研究从笛卡尔坐标系到椭圆坐标系的变换。

将上述变换写为如下形式：

根据

可得

如果将η视为常数，则上式代表焦距为c，离心率e=sec(η)的双曲线。

如果令

代入上式可得：

化简后有：

同理，如果我们将前面给出的椭圆坐标变换表示为如下形式：

因此有：

如果将ξ视为常数，则上述方程代表了焦距为c，离心率e=cosh^-1(ξ)的椭圆。

如果令

代入可得：

化简后有：

可以看到上式与前面表示变量p的一元二次方程在形式上是完全一致的。

注意到变量p和q的取值范围分别为0≤p≤1和q≤0，因此有p≥q，由此可以解得p和q分别为：

而根据前面变量p的定义，可得：

同理，根据前面变量q的定义，可得：

化简后可得：

解得

注意到q≤0，因此上式的两个根均为实根，并且由于ξ是非负的，因此可解得：

上式建立了笛卡尔坐标与ξ之间的联系。如果已知待构造椭圆的长轴和短轴分别为a和b，则我们可以考虑该椭圆在笛卡尔坐标系中的两个特殊点：（a, 0）和（0, b），显然这在这两个点上可以得到同一个q值：

将其代入上式即可求得对应椭圆网格区域边界的ξ值。

仍以前面的网格为例，取待生成的椭圆网格区域的长轴为a=2，短轴为b=1.2，则由上式可得计算得到对应的ξ取值约为0.6931，用于生成椭圆区域网格的矩形网格域的取值范围为ξ∈[0, 0.6931]，η∈[0, 2π]。因此，下面首先划分一块用于生成椭圆网格的规则网格区域，如下图所示。

用于坐标变换的正方形网格（边长为1）

这里为了保持程序的通用性，同时为了减小数值误差，可以将用于坐标变换的矩形网格的长度和宽度均设为1，这样在导入网格后首先进行一次比例变换即可使其满足相应的取值范围，然后再通过上述变换方法即可获得如下图所示的椭圆区域网格。

变换之后的椭圆区域网格（a=2，b=1.2）

可以看到，网格在变换之后仍然保持了非常好的正交性。利用有限元软件的网格质量工具对其检查可以发现其网格质量不存在任何问题。当然，这样的变换方法 会在椭圆焦点处生成少量的三角形单元，这些三角形单元是由四边形单元的某条单元边被折叠形成，因此本质上仍为四边形单元，其可能会对该区域的计算精度产生一定影响。

[1] C. Sun, Explicit equations to transform from Cartesian to elliptic coordinates, Mathematical Modelling and Applications. 2(4) (2017) 43-46.

[2] K. Yuan, Y. Jiang, J. Liu, et al, 2D weight functions of stress intensity factors for high aspect ratio semi-elliptical surface cracks in finite thickness plate, Theoretical and Applied Fracture Mechanics. 110 (2020) 102808.

%程序用于演示生成椭圆区域的网格
clear,clc
%*****************************参数输入*************************************
a=2;                         %椭圆长轴
b=1.2;                       %椭圆短轴
Length_slice_num=40;         %用于变换的矩形网格的长度分段数
Width_slice_num=40;          %用于变换的矩形网格的宽度分段数
%**************************************************************************

%***************************构造矩形网格域**********************************
Length_coord=linspace(0,1,Length_slice_num+1);
Width_coord=linspace(0,1,Width_slice_num+1);
[Length_COORD,Width_COORD]=meshgrid(Length_coord,Width_coord);
%生成网格节点构成信息
Elem_connect=zeros(Length_slice_num*Width_slice_num,4);
k=1;
for i=1:Width_slice_num
  for j=1:Length_slice_num
    Elem_connect(k,1)=((Length_slice_num+1)*(i-1))+j;
    Elem_connect(k,2)=((Length_slice_num+1)*(i-1))+j+1;
    Elem_connect(k,3)=((Length_slice_num+1)*i)+j+1;
    Elem_connect(k,4)=((Length_slice_num+1)*i)+j;
    k=k+1;
  end
end
%生成节点坐标信息
Xi=reshape(Length_COORD',[(Length_slice_num+1)*(Width_slice_num+1),1]);
Eta=reshape(Width_COORD',[(Length_slice_num+1)*(Width_slice_num+1),1]);
Node_Data=[Xi Eta];
%绘制矩形网格区域
Patch_Xi=zeros(4,size(Elem_connect,1));
Patch_Eta=zeros(4,size(Elem_connect,1));
for i=1:size(Elem_connect,1)   %针对单元的循环
  for j=1:4                      %针对单元构成节点的循环
    Patch_Xi(j,i)=Node_Data(Elem_connect(i,j),1);
    Patch_Eta(j,i)=Node_Data(Elem_connect(i,j),2);
  end
end
patch(Patch_Xi,Patch_Eta,1,...
  "FaceColor",[234/256 234/256 234/256]);axis equal;
title("变换之前的网格");
%**************************************************************************

%*********************基于坐标变换构造椭圆网格域****************************
c=sqrt(a^2-b^2);          %椭圆焦距
q=-b^2/c^2;
Xi_Bound=log(1-2*q+2*sqrt(q^2-q))/2;  %椭圆区域边界对应的椭圆坐标
Node_Data_Scale(:,1)=Node_Data(:,1)*Xi_Bound;
Node_Data_Scale(:,2)=Node_Data(:,2)*2*pi;
Node_Data_trans(:,1)=c*...
  cosh(Node_Data_Scale(:,1)).*cos(Node_Data_Scale(:,2));
Node_Data_trans(:,2)=c...
  *sinh(Node_Data_Scale(:,1)).*sin(Node_Data_Scale(:,2));
%绘制椭圆形网格域
Patch_X=zeros(4,size(Elem_connect,1));
Patch_Y=zeros(4,size(Elem_connect,1));
for i=1:size(Elem_connect,1)   %针对单元的循环
  for j=1:4                      %针对单元构成节点的循环
    Patch_X(j,i)=Node_Data_trans(Elem_connect(i,j),1);
    Patch_Y(j,i)=Node_Data_trans(Elem_connect(i,j),2);
  end
end
figure;patch(Patch_X,Patch_Y,1,...
  "FaceColor",[234/256 234/256 234/256]);axis equal;
title("变换之后的网格");
%**************************************************************************