椭圆型方程网格生成法

• 算法特征
  ①. 给定边界条件; ②. 在物理空间构建椭圆型方程; ③. 在计算空间求解椭圆型方程
• 算法推导
  以2维物理空间$(x, y)$及计算空间$(xi,eta)$为例, 不失一般性, 令存在如下变换关系,
  egin{equation}
  left{
  egin{split}
  xi = xi(x, y) \
  eta = eta(x, y)
  end{split}
  ight.
  label{eq_1}
  end{equation}
  其中, $(x, y)$为单连通任意形状区域, $(xi, eta)$为单连通矩形区域. 经过上式$eqref{eq_1}$变换, 物理空间边界对应计算空间边界, 物理空间内点对应计算空间内点.
  直接通过给定边界条件求解PDE以获取内点映射关系, 通常应当建立椭圆型方程. 进一步, 在物理空间中构建如下Poisson's equations具体描述式$eqref{eq_1}$之变换关系,
  egin{equation}
  egin{split}
  frac{partial^2xi}{partial x^2} + frac{partial^2xi}{partial y^2} &= P(x,y) \
  frac{partial^2eta}{partial x^2} + frac{partial^2eta}{partial y^2} &= Q(x,y)
  end{split}
  label{eq_2}
  end{equation}
  现考虑将物理空间之PDE转换至计算空间, 由于,
  egin{equation}
  left{
  egin{split}
  frac{partial x}{partial x} &= x_xixi_x + x_etaeta_x = 1 \
  frac{partial x}{partial y} &= x_xixi_y + x_etaeta_y = 0 \
  frac{partial y}{partial x} &= y_xixi_x + y_etaeta_x = 0 \
  frac{partial y}{partial y} &= y_xixi_y + y_etaeta_y = 1 \
  end{split}
  ight.
  label{eq_3}
  end{equation}
  求解上式$eqref{eq_3}$得,
  egin{equation}
  xi_x = frac{y_eta}{J}quad xi_y=frac{-x_eta}{J}quad eta_x=frac{-y_xi}{J}quad eta_y=frac{x_xi}{J}
  label{eq_4}
  end{equation}
  其中, $displaystyle J = |partial(x, y)/partial(xi,eta)| = x_xi y_eta - y_xi x_eta$.
  由于,
  egin{equation}
  left{
  egin{split}
  frac{partial^2x}{partial x^2}&=x_{xixi}xi_x^2 + x_{etaeta}eta_x^2 + 2x_{xieta}xi_xeta_x + x_xixi_{xx} + x_etaeta_{xx} = 0 \
  frac{partial^2x}{partial y^2}&=x_{xixi}xi_y^2 + x_{etaeta}eta_y^2 + 2x_{xieta}xi_yeta_y + x_xixi_{yy} + x_etaeta_{yy} = 0 \
  frac{partial^2y}{partial x^2}&=y_{xixi}xi_x^2 + y_{etaeta}eta_x^2 + 2y_{xieta}xi_xeta_x + y_xixi_{xx} + y_etaeta_{xx} = 0 \
  frac{partial^2y}{partial y^2}&=y_{xixi}xi_y^2 + y_{etaeta}eta_y^2 + 2y_{xieta}xi_yeta_y + y_xixi_{yy} + y_etaeta_{yy} = 0 \
  end{split}
  ight.
  label{eq_5}
  end{equation}
  结合式$eqref{eq_2}$、式$eqref{eq_4}$, 上式$eqref{eq_5}$转换如下,
  egin{equation}
  egin{split}
  alpha x_{xixi} - 2eta x_{xieta} + gamma x_{etaeta} + J^2(Px_xi+Qx_eta) = 0 \
  alpha y_{xixi} - 2eta y_{xieta} + gamma y_{etaeta} + J^2(Py_xi+Qy_eta) = 0 \
  end{split}
  label{eq_6}
  end{equation}
  其中, $alpha = x_eta^2 + y_eta^2$, $eta=x_xi x_eta+y_xi y_eta$, $gamma=x_xi^2 + y_xi^2$. 上式$eqref{eq_6}$即计算空间中待求解之椭圆型方程.
  由于大型方程组求解的本质困难, 本文拟采用含时演化技术将椭圆型方程转换为抛物型方程进行求解. 如果含时演化的抛物型方程最终趋于稳定, 则收敛于相应椭圆型方程的解. 式$eqref{eq_6}$椭圆型方程转换后的抛物型方程如下,
  egin{equation}
  egin{split}
  x_t &= alpha x_{xixi} - 2eta x_{xieta} + gamma x_{etaeta} + J^2(Px_xi+Qx_eta) = 0 \
  y_t &= alpha y_{xixi} - 2eta y_{xieta} + gamma y_{etaeta} + J^2(Py_xi+Qy_eta) = 0 \
  end{split}
  label{eq_7}
  end{equation}
  实际操作中, 通常将计算空间$(xi,eta)$取为物理空间$(x,y)$离散化后的下标值, 因此各方向相邻网格节点之间距$Deltaxi = Deltaeta = 1$. 本文拟采用如下差分格式离散化式$eqref{eq_7}$各微分项(以$x$为例),
  egin{equation}
  left{
  egin{split}
  x_xi &= frac{x_{i+1,j} - x_{i-1,j}}{2} \
  x_eta &= frac{x_{i,j+1} - x_{i, j-1}}{2} \
  x_{xixi} &= x_{i+1,j} + x_{i-1,j} - 2x_{i,j} \
  x_{etaeta} &= x_{i,j+1} + x_{i,j-1} - 2x_{i,j} \
  x_{xieta} &= frac{x_{i+1,j+1} + x_{i-1,j-1} - x_{i+1,j-1} - x_{i-1,j+1}}{4} \
  x_t &= frac{x^{k+1}-x^k}{h_t} \
  end{split}
  ight.
  label{eq_8}
  end{equation}
  注意, 式$eqref{eq_7}$之求解需要配合初始条件及边界条件, 本文重在获取其稳态解, 因此对初始条件不作具体要求, 边界条件则根据物理空间边界与计算空间边界之映射关系灵活给出. 若,
  egin{equation}
  |X^{k+1} - X^k| < epsilonquad &quad |Y^{k+1} - Y^k|<epsilon
  label{eq_9}
  end{equation}
  其中, $epsilon$取值足够小正数, 则$X^k$、$Y^k$均趋于稳定并收敛于式$eqref{eq_6}$的解, 此即物理空间中根据椭圆型方程生成的网格节点.
  
• 代码实现
  现以如下翼型为例进行算法实施,
  egin{equation*}
  y=0.1781sqrt{x}-0.0756x-0.2122x^2+0.1705x^3-0.0609x^4, quad ext{where }xin[0,1]
  end{equation*}
  令$P=Q=0$, 通过求解无源Laplace's eqautions进行网格生成, 具体实现如下,
  

  # 椭圆型方程网格生成法
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  # 对称翼型的上半部
  def half_wing(x):
      funcVal = 0.1781 * numpy.sqrt(x) - 0.0756 * x - 0.2122 * x ** 2 + 0.1705 * x ** 3 - 0.0609 * x ** 4
      return funcVal
      
  
  # 构造物理空间与计算空间边界映射关系
  def build_bdy_maps():
      p2c_xi = [
          ((1.0, 0.0), (0, 15)),
          ((0.9, half_wing(0.9)), (1, 15)),
          ((0.8, half_wing(0.8)), (2, 15)),
          ((0.7, half_wing(0.7)), (3, 15)),
          ((0.6, half_wing(0.6)), (4, 15)),
          ((0.5, half_wing(0.5)), (5, 15)),
          ((0.4, half_wing(0.4)), (6, 15)),
          ((0.3, half_wing(0.3)), (7, 15)),
          ((0.2, half_wing(0.2)), (8, 15)),
          ((0.1, half_wing(0.1)), (9, 15)),
          ((0.0, half_wing(0.0)), (10, 15)),
          ((0.1, -half_wing(0.1)), (11, 15)),
          ((0.2, -half_wing(0.2)), (12, 15)),
          ((0.3, -half_wing(0.3)), (13, 15)),
          ((0.4, -half_wing(0.4)), (14, 15)),
          ((0.5, -half_wing(0.5)), (15, 15)),
          ((0.6, -half_wing(0.6)), (16, 15)),
          ((0.7, -half_wing(0.7)), (17, 15)),
          ((0.8, -half_wing(0.8)), (18, 15)),
          ((0.9, -half_wing(0.9)), (19, 15)),
          ((1.0, 0.0), (20, 15)),
          ((4.0, 0.0), (0, 0)),
          ((4.0, 1.0), (1, 0)),
          ((4.0, 2.0), (2, 0)),
          ((3.0, 2.0), (3, 0)),
          ((2.0, 2.0), (4, 0)),
          ((1.0, 2.0), (5, 0)),
          ((0.0, 2.0), (6, 0)),
          ((-1.0, 2.0), (7, 0)),
          ((-2.0, 2.0), (8, 0)),
          ((-2.0, 1.0), (9, 0)),
          ((-2.0, 0.0), (10, 0)),
          ((-2.0, -1.0), (11, 0)),
          ((-2.0, -2.0), (12, 0)),
          ((-1.0, -2.0), (13, 0)),
          ((0.0, -2.0), (14, 0)),
          ((1.0, -2.0), (15, 0)),
          ((2.0, -2.0), (16, 0)),
          ((3.0, -2.0), (17, 0)),
          ((4.0, -2.0), (18, 0)),
          ((4.0, -1.0), (19, 0)),
          ((4.0, 0.0), (20, 0))
      ]
      
      p2c_eta = [
          ((4.0, 0.0), (0, 0)),
          ((3.8, 0.0), (0, 1)),
          ((3.6, 0.0), (0, 2)),
          ((3.4, 0.0), (0, 3)),
          ((3.2, 0.0), (0, 4)),
          ((3.0, 0.0), (0, 5)),
          ((2.8, 0.0), (0, 6)),
          ((2.6, 0.0), (0, 7)),
          ((2.4, 0.0), (0, 8)),
          ((2.2, 0.0), (0, 9)),
          ((2.0, 0.0), (0, 10)),
          ((1.8, 0.0), (0, 11)),
          ((1.6, 0.0), (0, 12)),
          ((1.4, 0.0), (0, 13)),
          ((1.2, 0.0), (0, 14)),
          ((1.0, 0.0), (0, 15)),
          ((4.0, 0.0), (20, 0)),
          ((3.8, 0.0), (20, 1)),
          ((3.6, 0.0), (20, 2)),
          ((3.4, 0.0), (20, 3)),
          ((3.2, 0.0), (20, 4)),
          ((3.0, 0.0), (20, 5)),
          ((2.8, 0.0), (20, 6)),
          ((2.6, 0.0), (20, 7)),
          ((2.4, 0.0), (20, 8)),
          ((2.2, 0.0), (20, 9)),
          ((2.0, 0.0), (20, 10)),
          ((1.8, 0.0), (20, 11)),
          ((1.6, 0.0), (20, 12)),
          ((1.4, 0.0), (20, 13)),
          ((1.2, 0.0), (20, 14)),
          ((1.0, 0.0), (20, 15))
      ]
      
      return p2c_xi, p2c_eta
      
      
  class EllipticGrid(object):
      
      def __init__(self, p2c_xi, p2c_eta):
          self.__p2c_xi = p2c_xi                           # 物理空间与计算空间xi轴边界映射关系
          self.__p2c_eta = p2c_eta                         # 物理空间与计算空间eta轴边界映射关系
      
          self.__n_xi, self.__n_eta = self.__get_n()       # 计算空间子区间划分数
          self.__ht = 0.1
          
          self.__bdyX, self.__bdyY = self.__get_bdy()
      
      
      def get_solu(self, maxIter=1000000, epsilon=1.e-9):
          '''
          数值求解
          maxIter: 最大迭代次数
          epsilon: 收敛判据
          '''
          X0 = self.__get_initX()
          Y0 = self.__get_initY()
          
          for i in range(maxIter):
              Xx = self.__calc_Ux(X0)
              Xe = self.__calc_Ue(X0)
              Xxx = self.__calc_Uxx(X0)
              Xee = self.__calc_Uee(X0)
              Xxe = self.__calc_Uxe(X0)
              Yx = self.__calc_Ux(Y0)
              Ye = self.__calc_Ue(Y0)
              Yxx = self.__calc_Uxx(Y0)
              Yee = self.__calc_Uee(Y0)
              Yxe = self.__calc_Uxe(Y0)
              
              alpha = self.__calc_alpha(Xe, Ye)
              beta = self.__calc_beta(Xx, Xe, Yx, Ye)
              gamma = self.__calc_gamma(Xx, Yx)
              
              Kx = alpha * Xxx - 2 * beta * Xxe + gamma * Xee
              Ky = alpha * Yxx - 2 * beta * Yxe + gamma * Yee
              
              Xk = self.__calc_Uk(X0, Kx, self.__ht)
              Yk = self.__calc_Uk(Y0, Ky, self.__ht)
              self.__fill_bdyX(Xk)
              self.__fill_bdyY(Yk)
              
              # print(i, numpy.linalg.norm(Xk - X0, numpy.inf), numpy.linalg.norm(Yk - Y0, numpy.inf))
              if self.__converged(Xk - X0, Yk - Y0, epsilon):
                  break
                  
              X0 = Xk
              Y0 = Yk
          else:
              raise Exception(">>> Not converged after {} iterations! <<<".format(maxIter))
          
          return Xk, Yk        
                  
                  
      def __converged(self, deltaX, deltaY, epsilon):
          normVal1 = numpy.linalg.norm(deltaX, numpy.inf)
          normVal2 = numpy.linalg.norm(deltaY, numpy.inf)
          if normVal1 < epsilon and normVal2 < epsilon:
              return True
          return False
          
          
      def __calc_Ux(self, U):
          Ux = numpy.zeros(U.shape)
          Ux[:, 1:-1] = (U[:, 2:] - U[:, :-2]) / 2
          return Ux
          
          
      def __calc_Ue(self, U):
          Ue = numpy.zeros(U.shape)
          Ue[1:-1, :] = (U[2:, :] - U[:-2, :]) / 2
          return Ue
          
          
      def __calc_Uxx(self, U):
          Uxx = numpy.zeros(U.shape)
          Uxx[:, 1:-1] = U[:, 2:] + U[:, :-2] - 2 * U[:, 1:-1]
          return Uxx
          
          
      def __calc_Uee(self, U):
          Uee = numpy.zeros(U.shape)
          Uee[1:-1, :] = U[2:, :] + U[:-2, :] - 2 * U[1:-1, :]
          return Uee
          
          
      def __calc_Uxe(self, U):
          Uxe = numpy.zeros(U.shape)
          Uxe[1:-1, 1:-1] = (U[2:, 2:] + U[:-2, :-2] - U[:-2, 2:] - U[2:, :-2]) / 4
          return Uxe
          
          
      def __calc_alpha(self, Xe, Ye):
          alpha = Xe ** 2 + Ye ** 2
          return alpha
          
          
      def __calc_beta(self, Xx, Xe, Yx, Ye):
          beta = Xx * Xe + Yx * Ye
          return beta
          
          
      def __calc_gamma(self, Xx, Yx):
          gamma = Xx ** 2 + Yx ** 2
          return gamma
          
          
      def __calc_Uk(self, U, K, ht):
          Uk = U + K * ht
          return Uk
          
          
      def __get_bdy(self):
          '''
          获取边界条件
          '''
          bdyX = numpy.zeros((self.__n_eta + 1, self.__n_xi + 1))
          bdyY = numpy.zeros((self.__n_eta + 1, self.__n_xi + 1))
          
          for XY, XE in self.__p2c_xi:
              bdyX[XE[1], XE[0]] = XY[0]
              bdyY[XE[1], XE[0]] = XY[1]
              
          for XY, XE in self.__p2c_eta:
              bdyX[XE[1], XE[0]] = XY[0]
              bdyY[XE[1], XE[0]] = XY[1]
              
          return bdyX, bdyY
          
          
      def __get_initX(self):
          '''
          获取X之初始条件
          '''
          X0 = numpy.zeros(self.__bdyX.shape)
          self.__fill_bdyX(X0)
          return X0
          
          
      def __get_initY(self):
          '''
          获取Y之初始条件
          '''
          Y0 = numpy.zeros(self.__bdyY.shape)
          self.__fill_bdyY(Y0)
          return Y0
          
          
      def __fill_bdyX(self, U):
          '''
          填充X之边界条件
          '''
          U[:, 0] = self.__bdyX[:, 0]
          U[:, -1] = self.__bdyX[:, -1]
          U[0, :] = self.__bdyX[0, :]
          U[-1, :] = self.__bdyX[-1, :]
          
          
      def __fill_bdyY(self, U):
          '''
          填充Y之边界条件
          '''
          U[:, 0] = self.__bdyY[:, 0]
          U[:, -1] = self.__bdyY[:, -1]
          U[0, :] = self.__bdyY[0, :]
          U[-1, :] = self.__bdyY[-1, :]
      
          
      def __get_n(self):
          arr_xi = numpy.array(list(item[1] for item in self.__p2c_xi))
          n_xi = numpy.max(arr_xi[:, 0])
          n_eta = numpy.max(arr_xi[:, 1])
          return n_xi, n_eta
          
          
  class EGPlot(object):
      
      @staticmethod
      def plot_fig(egObj):
          maxIter = 1000000
          epsilon = 1.e-9
          
          X, Y = egObj.get_solu(maxIter, epsilon)
          
          fig = plt.figure(figsize=(6, 4))
          ax1 = plt.subplot()
          
          ax1.plot(X[:, 0], Y[:, 0], c="red", lw=1, label="mapping to $\eta$")
          ax1.plot(X[:, -1], Y[:, -1], c="red", lw=1)
          n_eta, n_xi = X.shape
          for colIdx in range(1, n_xi-1):
              tmpX = X[:, colIdx]
              tmpY = Y[:, colIdx]
              ax1.plot(tmpX, tmpY, "k-", lw=1)
              
          ax1.plot(X[0, :], Y[0, :], c="green", lw=1, label="mapping to $\xi$")
          ax1.plot(X[-1, :], Y[-1, :], c="green", lw=1)
          for rowIdx in range(1, n_eta-1):
              tmpX = X[rowIdx, :]
              tmpY = Y[rowIdx, :]
              ax1.plot(tmpX, tmpY, "k-", lw=1)
          ax1.legend()
          
          fig.tight_layout()
          fig.savefig("plot_fig.png", dpi=100)
  
  
  
  if __name__ == "__main__":
      p2c_xi, p2c_eta = build_bdy_maps()
      egObj = EllipticGrid(p2c_xi, p2c_eta)
      
      EGPlot.plot_fig(egObj)

• 结果展示
   
  左图为翼型所处物理空间, 右图为利用椭圆型方程在此物理空间生成之网格. 图中绿色与红色曲线分别映射至计算空间相应边界.
• 使用建议
  ①. 时间步长$h_t$选择过小收敛速度慢, 选择过大容易发散;
  ②. 物理空间与计算空间边界映射关系可灵活选择;
  ③. 多连通物理空间需适当剖分转单连通空间.
• 参考文档
  Thompson, Joe F., Zahir UA Warsi, and C. Wayne Mastin. Numerical grid generation: foundations and applications. Elsevier North-Holland, Inc., 1985.