• 周六. 7月 2nd, 2022

5G编程聚合网

5G时代下一个聚合的编程学习网

热门标签

椭圆型方程网格生成法

admin

11月 28, 2021
  • 算法特征
    ①. 给定边界条件; ②. 在物理空间构建椭圆型方程; ③. 在计算空间求解椭圆型方程
  • 算法推导
    以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进行网格生成, 具体实现如下,

      1 # 椭圆型方程网格生成法
      2 
      3 import numpy
      4 from matplotlib import pyplot as plt
      5 
      6 
      7 # 对称翼型的上半部
      8 def half_wing(x):
      9     funcVal = 0.1781 * numpy.sqrt(x) - 0.0756 * x - 0.2122 * x ** 2 + 0.1705 * x ** 3 - 0.0609 * x ** 4
     10     return funcVal
     11     
     12 
     13 # 构造物理空间与计算空间边界映射关系
     14 def build_bdy_maps():
     15     p2c_xi = [
     16         ((1.0, 0.0), (0, 15)),
     17         ((0.9, half_wing(0.9)), (1, 15)),
     18         ((0.8, half_wing(0.8)), (2, 15)),
     19         ((0.7, half_wing(0.7)), (3, 15)),
     20         ((0.6, half_wing(0.6)), (4, 15)),
     21         ((0.5, half_wing(0.5)), (5, 15)),
     22         ((0.4, half_wing(0.4)), (6, 15)),
     23         ((0.3, half_wing(0.3)), (7, 15)),
     24         ((0.2, half_wing(0.2)), (8, 15)),
     25         ((0.1, half_wing(0.1)), (9, 15)),
     26         ((0.0, half_wing(0.0)), (10, 15)),
     27         ((0.1, -half_wing(0.1)), (11, 15)),
     28         ((0.2, -half_wing(0.2)), (12, 15)),
     29         ((0.3, -half_wing(0.3)), (13, 15)),
     30         ((0.4, -half_wing(0.4)), (14, 15)),
     31         ((0.5, -half_wing(0.5)), (15, 15)),
     32         ((0.6, -half_wing(0.6)), (16, 15)),
     33         ((0.7, -half_wing(0.7)), (17, 15)),
     34         ((0.8, -half_wing(0.8)), (18, 15)),
     35         ((0.9, -half_wing(0.9)), (19, 15)),
     36         ((1.0, 0.0), (20, 15)),
     37         ((4.0, 0.0), (0, 0)),
     38         ((4.0, 1.0), (1, 0)),
     39         ((4.0, 2.0), (2, 0)),
     40         ((3.0, 2.0), (3, 0)),
     41         ((2.0, 2.0), (4, 0)),
     42         ((1.0, 2.0), (5, 0)),
     43         ((0.0, 2.0), (6, 0)),
     44         ((-1.0, 2.0), (7, 0)),
     45         ((-2.0, 2.0), (8, 0)),
     46         ((-2.0, 1.0), (9, 0)),
     47         ((-2.0, 0.0), (10, 0)),
     48         ((-2.0, -1.0), (11, 0)),
     49         ((-2.0, -2.0), (12, 0)),
     50         ((-1.0, -2.0), (13, 0)),
     51         ((0.0, -2.0), (14, 0)),
     52         ((1.0, -2.0), (15, 0)),
     53         ((2.0, -2.0), (16, 0)),
     54         ((3.0, -2.0), (17, 0)),
     55         ((4.0, -2.0), (18, 0)),
     56         ((4.0, -1.0), (19, 0)),
     57         ((4.0, 0.0), (20, 0))
     58     ]
     59     
     60     p2c_eta = [
     61         ((4.0, 0.0), (0, 0)),
     62         ((3.8, 0.0), (0, 1)),
     63         ((3.6, 0.0), (0, 2)),
     64         ((3.4, 0.0), (0, 3)),
     65         ((3.2, 0.0), (0, 4)),
     66         ((3.0, 0.0), (0, 5)),
     67         ((2.8, 0.0), (0, 6)),
     68         ((2.6, 0.0), (0, 7)),
     69         ((2.4, 0.0), (0, 8)),
     70         ((2.2, 0.0), (0, 9)),
     71         ((2.0, 0.0), (0, 10)),
     72         ((1.8, 0.0), (0, 11)),
     73         ((1.6, 0.0), (0, 12)),
     74         ((1.4, 0.0), (0, 13)),
     75         ((1.2, 0.0), (0, 14)),
     76         ((1.0, 0.0), (0, 15)),
     77         ((4.0, 0.0), (20, 0)),
     78         ((3.8, 0.0), (20, 1)),
     79         ((3.6, 0.0), (20, 2)),
     80         ((3.4, 0.0), (20, 3)),
     81         ((3.2, 0.0), (20, 4)),
     82         ((3.0, 0.0), (20, 5)),
     83         ((2.8, 0.0), (20, 6)),
     84         ((2.6, 0.0), (20, 7)),
     85         ((2.4, 0.0), (20, 8)),
     86         ((2.2, 0.0), (20, 9)),
     87         ((2.0, 0.0), (20, 10)),
     88         ((1.8, 0.0), (20, 11)),
     89         ((1.6, 0.0), (20, 12)),
     90         ((1.4, 0.0), (20, 13)),
     91         ((1.2, 0.0), (20, 14)),
     92         ((1.0, 0.0), (20, 15))
     93     ]
     94     
     95     return p2c_xi, p2c_eta
     96     
     97     
     98 class EllipticGrid(object):
     99     
    100     def __init__(self, p2c_xi, p2c_eta):
    101         self.__p2c_xi = p2c_xi                           # 物理空间与计算空间xi轴边界映射关系
    102         self.__p2c_eta = p2c_eta                         # 物理空间与计算空间eta轴边界映射关系
    103     
    104         self.__n_xi, self.__n_eta = self.__get_n()       # 计算空间子区间划分数
    105         self.__ht = 0.1
    106         
    107         self.__bdyX, self.__bdyY = self.__get_bdy()
    108     
    109     
    110     def get_solu(self, maxIter=1000000, epsilon=1.e-9):
    111         '''
    112         数值求解
    113         maxIter: 最大迭代次数
    114         epsilon: 收敛判据
    115         '''
    116         X0 = self.__get_initX()
    117         Y0 = self.__get_initY()
    118         
    119         for i in range(maxIter):
    120             Xx = self.__calc_Ux(X0)
    121             Xe = self.__calc_Ue(X0)
    122             Xxx = self.__calc_Uxx(X0)
    123             Xee = self.__calc_Uee(X0)
    124             Xxe = self.__calc_Uxe(X0)
    125             Yx = self.__calc_Ux(Y0)
    126             Ye = self.__calc_Ue(Y0)
    127             Yxx = self.__calc_Uxx(Y0)
    128             Yee = self.__calc_Uee(Y0)
    129             Yxe = self.__calc_Uxe(Y0)
    130             
    131             alpha = self.__calc_alpha(Xe, Ye)
    132             beta = self.__calc_beta(Xx, Xe, Yx, Ye)
    133             gamma = self.__calc_gamma(Xx, Yx)
    134             
    135             Kx = alpha * Xxx - 2 * beta * Xxe + gamma * Xee
    136             Ky = alpha * Yxx - 2 * beta * Yxe + gamma * Yee
    137             
    138             Xk = self.__calc_Uk(X0, Kx, self.__ht)
    139             Yk = self.__calc_Uk(Y0, Ky, self.__ht)
    140             self.__fill_bdyX(Xk)
    141             self.__fill_bdyY(Yk)
    142             
    143             # print(i, numpy.linalg.norm(Xk - X0, numpy.inf), numpy.linalg.norm(Yk - Y0, numpy.inf))
    144             if self.__converged(Xk - X0, Yk - Y0, epsilon):
    145                 break
    146                 
    147             X0 = Xk
    148             Y0 = Yk
    149         else:
    150             raise Exception(">>> Not converged after {} iterations! <<<".format(maxIter))
    151         
    152         return Xk, Yk        
    153                 
    154                 
    155     def __converged(self, deltaX, deltaY, epsilon):
    156         normVal1 = numpy.linalg.norm(deltaX, numpy.inf)
    157         normVal2 = numpy.linalg.norm(deltaY, numpy.inf)
    158         if normVal1 < epsilon and normVal2 < epsilon:
    159             return True
    160         return False
    161         
    162         
    163     def __calc_Ux(self, U):
    164         Ux = numpy.zeros(U.shape)
    165         Ux[:, 1:-1] = (U[:, 2:] - U[:, :-2]) / 2
    166         return Ux
    167         
    168         
    169     def __calc_Ue(self, U):
    170         Ue = numpy.zeros(U.shape)
    171         Ue[1:-1, :] = (U[2:, :] - U[:-2, :]) / 2
    172         return Ue
    173         
    174         
    175     def __calc_Uxx(self, U):
    176         Uxx = numpy.zeros(U.shape)
    177         Uxx[:, 1:-1] = U[:, 2:] + U[:, :-2] - 2 * U[:, 1:-1]
    178         return Uxx
    179         
    180         
    181     def __calc_Uee(self, U):
    182         Uee = numpy.zeros(U.shape)
    183         Uee[1:-1, :] = U[2:, :] + U[:-2, :] - 2 * U[1:-1, :]
    184         return Uee
    185         
    186         
    187     def __calc_Uxe(self, U):
    188         Uxe = numpy.zeros(U.shape)
    189         Uxe[1:-1, 1:-1] = (U[2:, 2:] + U[:-2, :-2] - U[:-2, 2:] - U[2:, :-2]) / 4
    190         return Uxe
    191         
    192         
    193     def __calc_alpha(self, Xe, Ye):
    194         alpha = Xe ** 2 + Ye ** 2
    195         return alpha
    196         
    197         
    198     def __calc_beta(self, Xx, Xe, Yx, Ye):
    199         beta = Xx * Xe + Yx * Ye
    200         return beta
    201         
    202         
    203     def __calc_gamma(self, Xx, Yx):
    204         gamma = Xx ** 2 + Yx ** 2
    205         return gamma
    206         
    207         
    208     def __calc_Uk(self, U, K, ht):
    209         Uk = U + K * ht
    210         return Uk
    211         
    212         
    213     def __get_bdy(self):
    214         '''
    215         获取边界条件
    216         '''
    217         bdyX = numpy.zeros((self.__n_eta + 1, self.__n_xi + 1))
    218         bdyY = numpy.zeros((self.__n_eta + 1, self.__n_xi + 1))
    219         
    220         for XY, XE in self.__p2c_xi:
    221             bdyX[XE[1], XE[0]] = XY[0]
    222             bdyY[XE[1], XE[0]] = XY[1]
    223             
    224         for XY, XE in self.__p2c_eta:
    225             bdyX[XE[1], XE[0]] = XY[0]
    226             bdyY[XE[1], XE[0]] = XY[1]
    227             
    228         return bdyX, bdyY
    229         
    230         
    231     def __get_initX(self):
    232         '''
    233         获取X之初始条件
    234         '''
    235         X0 = numpy.zeros(self.__bdyX.shape)
    236         self.__fill_bdyX(X0)
    237         return X0
    238         
    239         
    240     def __get_initY(self):
    241         '''
    242         获取Y之初始条件
    243         '''
    244         Y0 = numpy.zeros(self.__bdyY.shape)
    245         self.__fill_bdyY(Y0)
    246         return Y0
    247         
    248         
    249     def __fill_bdyX(self, U):
    250         '''
    251         填充X之边界条件
    252         '''
    253         U[:, 0] = self.__bdyX[:, 0]
    254         U[:, -1] = self.__bdyX[:, -1]
    255         U[0, :] = self.__bdyX[0, :]
    256         U[-1, :] = self.__bdyX[-1, :]
    257         
    258         
    259     def __fill_bdyY(self, U):
    260         '''
    261         填充Y之边界条件
    262         '''
    263         U[:, 0] = self.__bdyY[:, 0]
    264         U[:, -1] = self.__bdyY[:, -1]
    265         U[0, :] = self.__bdyY[0, :]
    266         U[-1, :] = self.__bdyY[-1, :]
    267     
    268         
    269     def __get_n(self):
    270         arr_xi = numpy.array(list(item[1] for item in self.__p2c_xi))
    271         n_xi = numpy.max(arr_xi[:, 0])
    272         n_eta = numpy.max(arr_xi[:, 1])
    273         return n_xi, n_eta
    274         
    275         
    276 class EGPlot(object):
    277     
    278     @staticmethod
    279     def plot_fig(egObj):
    280         maxIter = 1000000
    281         epsilon = 1.e-9
    282         
    283         X, Y = egObj.get_solu(maxIter, epsilon)
    284         
    285         fig = plt.figure(figsize=(6, 4))
    286         ax1 = plt.subplot()
    287         
    288         ax1.plot(X[:, 0], Y[:, 0], c="red", lw=1, label="mapping to $\eta$")
    289         ax1.plot(X[:, -1], Y[:, -1], c="red", lw=1)
    290         n_eta, n_xi = X.shape
    291         for colIdx in range(1, n_xi-1):
    292             tmpX = X[:, colIdx]
    293             tmpY = Y[:, colIdx]
    294             ax1.plot(tmpX, tmpY, "k-", lw=1)
    295             
    296         ax1.plot(X[0, :], Y[0, :], c="green", lw=1, label="mapping to $\xi$")
    297         ax1.plot(X[-1, :], Y[-1, :], c="green", lw=1)
    298         for rowIdx in range(1, n_eta-1):
    299             tmpX = X[rowIdx, :]
    300             tmpY = Y[rowIdx, :]
    301             ax1.plot(tmpX, tmpY, "k-", lw=1)
    302         ax1.legend()
    303         
    304         fig.tight_layout()
    305         fig.savefig("plot_fig.png", dpi=100)
    306 
    307 
    308 
    309 if __name__ == "__main__":
    310     p2c_xi, p2c_eta = build_bdy_maps()
    311     egObj = EllipticGrid(p2c_xi, p2c_eta)
    312     
    313     EGPlot.plot_fig(egObj)

    View Code

  • 结果展示

     
    左图为翼型所处物理空间, 右图为利用椭圆型方程在此物理空间生成之网格. 图中绿色与红色曲线分别映射至计算空间相应边界.


  • 使用建议

    ①. 时间步长$h_t$选择过小收敛速度慢, 选择过大容易发散;
    ②. 物理空间与计算空间边界映射关系可灵活选择;
    ③. 多连通物理空间需适当剖分转单连通空间.

  • 参考文档
    Thompson, Joe F., Zahir UA Warsi, and C. Wayne Mastin. Numerical grid generation: foundations and applications. Elsevier North-Holland, Inc., 1985.

发表评论

您的电子邮箱地址不会被公开。