Search

Article

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

Effect of equilibrium distribution positivity on numerical performance of lattice Boltzmann method

Ye Huan-Feng Jin Di Kuang Bo Yang Yan-Hua

Citation:

Effect of equilibrium distribution positivity on numerical performance of lattice Boltzmann method

Ye Huan-Feng, Jin Di, Kuang Bo, Yang Yan-Hua
PDF
HTML
Get Citation
  • Traditionally, the numerical performance of the lattice Boltzmann method is mainly determined by the moment degree of a discrete equilibrium distribution. The equilibrium distribution positivity is merely considered as an ancillary property which is used to constrict the numerical configuration. With the newly-developed partial Gaussian-Hermite quadrature scheme, the positivity of equilibrium distribution is validated as an independent property like moment degree which can be adjusted by discrete velocities. Researchers speculated that the positivity should also be significant for the numerical performance by the lattice Boltzmann method and can be used to improve the performance. Comparing with the classical improvement through moment degree, the positivity approach will not bring additional computation. However, due to the lack of boundary treatment, the speculation has not been validated by detailed numerical simulations. In this paper, through employing a periodic case, the Taylor-Green vortex, to avoid the boundary issue, we in depth analyze the numerical effect of the model positivity, including the numerical accuracy in the model positive range, the influence of positivity on the numerical performance, and the significance comparison between positivity and moment degree. The results show that for a given model, the numerical accuracy is not consistent in the whole positive range. As the configuration is close to the border of positive range, the accuracy will degrade though it is still acceptable. The numerical performance of a model depends on both moment degree and positivity. The role that the moment degree plays lies mainly in the qualification of a model on Galilean invariance. Once a model fulfills the Galilean invariance, its numerical performance is solely dependent on the positivity. Hence, the improvement approach through modifying the model positivity is a viable solution, and a Galilean invariant model with wider positive range does possess a better numerical performance regardless of its moment degree. Furthermore, based on the numerical results in this paper, all DnHm models with high moment degree are better than the classical D2Q9 model. Of the above models, the D2H3-2 model has the best performance and deserves to be further studied
      Corresponding author: Ye Huan-Feng, huanfye@163.com
    • Funds: Project supported by the National Science and Technology Major Project of the Ministry of Science and Technology of China (Grant No. 2017ZX06002002)
    [1]

    Ammar S, Pernaudat G, Trepanier J 2017 J. Comput. Phys. 343 73Google Scholar

    [2]

    Li Q, Luo K H, Kang Q J, He Y L, Chen Q, Liu Q 2016 Prog. Energ. Combust. 52 62Google Scholar

    [3]

    Shan X W 2010 Phys. Rev. E 81 045701Google Scholar

    [4]

    Shan X W, Chen H D 1993 Phys. Rev. E 47 1815Google Scholar

    [5]

    Sheikholeslami M, Gorji-Bandpy M, Ganji D D 2014 Powder Technol. 254 82Google Scholar

    [6]

    Zhang X L, Xiao L Z, Shan X W, Guo L 2014 Sci. Rep. 4 4843

    [7]

    Shan X W 2011 IMA J. Appl. Math. 76 650Google Scholar

    [8]

    Kang Q J, Zhang D X, Chen S Y, He X Y 2002 Phys. Rev. E 65 036318Google Scholar

    [9]

    Mountrakis L, Lorenz E, Malaspinas O, Alowayyed S, Chopard B, Hoekstra A G 2015 J. Comput. Sci. 9 45Google Scholar

    [10]

    Ladd A J C, Verberg R 2001 J. Stat. Phys. 104 1191Google Scholar

    [11]

    Chopard B, Ouared R, Deutsch A, Hatzikirou H, Wolf-Gladrow D 2010 Acta Biotheor. 58 329Google Scholar

    [12]

    Hirabayashi M, Ohta M, Rjfenacht D A, Chopard B 2003 Phys. Rev. E 68 021918Google Scholar

    [13]

    Latt J, Chopard B 2007 Int. J. Mod. Phys. C 18 619Google Scholar

    [14]

    Yu H D, Luo L S, Girimaji S S 2006 Comput. Fluids 35 957Google Scholar

    [15]

    Yu H D, Girimaji S S, Luo L S 2005 J. Comput. Phys. 209 599Google Scholar

    [16]

    He X Y, Chen S Y, Doolen G D 1998 J. Comput. Phys. 146 282Google Scholar

    [17]

    He X Y, Luo L S 1997 Phys. Rev. E 56 6811Google Scholar

    [18]

    Ye H F, Kuang B, Yang Y H 2019 Chin. Phys. B 28 014701Google Scholar

    [19]

    Dellar P J 2014 J. Comput. Phys. 259 270Google Scholar

    [20]

    Nie X B, Shan X W, Chen H D 2008 EPL-Europhys. Lett. 81 34005Google Scholar

    [21]

    Hazi G, Kavran P 2006 J. Phys. A 39 3127Google Scholar

    [22]

    Shan X W 2016 J. Comput. Sci. 17 475Google Scholar

    [23]

    Shan X W 2010 Phys. Rev. E 81 036702Google Scholar

    [24]

    Shim J W 2013 Phys. Rev. E 87 013312Google Scholar

    [25]

    Ye H F, Gan Z C, Kuang B, Yang Y H 2019 Chin. Phys. B 28 054702Google Scholar

    [26]

    Shan X W, He X Y 1998 Phys. Rev. Lett. 80 65Google Scholar

  • 图 1  格子Boltzmann离散模型正值性分析 虚线框表示模型正值区域, 其涵盖区域为$ [-U_{\rm pos},U_{\rm pos}]\times[-U_{\rm pos},U_{\rm pos}] $, 模型$ U_{\rm pos} $标注在图中模型名称之后

    Figure 1.  Positivity analysis of lattice Boltzmann discrete models. The dash squares denote the positive areas of each lattice Boltzmann discrete model, which are constructed as $ [-U_{\rm pos},U_{\rm pos}]\times[-U_{\rm pos},U_{\rm pos}] $. The $ U_{\rm pos} $ values of each model are annotated after the model names.

    图 2  各格子Boltzmann模型在Taylor-Green涡设计计算工况下初始时刻计算流域$ [0,2{\text{π}}]\times[0,2{\text{π}}] $$ \kappa $值色阶图 图中列对应格子Boltzmann模型, 行则对应表3中计算工况; 每个色阶图上方的数值$ \Delta u $值, 即正值区域值$ u_{\rm pos} $设计计算工况$ u_{L{\rm B} ,0} $差值;为方便观察正值性随计算工况遭到破坏的变化过程, 图中对$ \kappa $值做了截断处理, 仅区别$ [1.0,2.0] $范围变化

    Figure 2.  The initial $ \kappa $ filled contours of designed Taylor-Green vortex simulations for each lattice Boltzmann model in the fluid domain $ [0,2{\text{π}}]\times[0,2{\text{π}}] $. Each column is a simulation set of a lattice Boltzmann model with different designed cases. Each row is a set of a simulation case with different lattice Boltzmann models. The values in the contour panels are the values of $ \Delta u $, i.e. the difference between the $ u_{\rm pos} $ and the designed $ u_{L{\rm B} ,0} $. For the sake of identifying the positivity violation process of a lattice Boltzmann model, the filled contours truncate the range of $ \kappa $, limiting into $ [1.0,2.0] $.

    图 3  Taylor-Green涡半衰时刻$ t = 3.4657359\; \rm{s} $, $ [0.3775,0.6225]{\text{π}} \times[0.3775,0.6225]{\text{π}} $流域内, 不同计算工况下理论解(Theoretical) 与各格子Boltzmann模型的速度平方和$ U^2 $等高线图; 图中列对应格子Boltzmann模型, 行则对应表3中计算工况

    Figure 3.  The $ U^2 $ contour lines of Taylor-Green vortex in fluid domain $ [0.3775,0.6225]{\text{π}} \times[0.3775,0.6225]{\text{π}} $ at half-value decay time $ t = 3.4657359 \; \rm{s} $, including the theoretical solution (denoted as “Theoretical”) and the numerical results of lattice Boltzmann models. Each column is a set of simulations under the designed cases. Each row is a set of simulations under selected lattice Boltzmann models including the theoretical solution.

    图 4  各Taylor-Green涡设计计算工况下, 各格子Boltzmann模型计算误差随时间演变, 图中标号a, b, c, d分别对应表3各设计计算工况; 为方便表示, 图中横坐标单位为$ 0.173286795\; \rm{s} $

    Figure 4.  The error evolution of each simulation. The panel a, b, c, d renders each designed case in Table 3 respectively. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

    图 5  Taylor-Green涡设计计算工况下, $ \Delta u $值为正值的各格子Boltzmann模型数值模拟误差对比. 图中实线, 虚线和虚点线分别为工况a, b, c (见表3)结果, 离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s

    Figure 5.  The error comparison of Taylor-Green vortex simulations with postive value of $ \Delta u $. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed and dash-dotted line indicates case a, b and c in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

    图 6  Taylor-Green涡数值模拟误差与$ \Delta u $值相关性分析 (a)对比了$ \Delta u $接近的D2H4模型a, b, c工况与D2H5模型b, c, d工况计算误差; (b)对比了$ \Delta u $$ -0.30 $附近所有数值模拟计算误差. 所有图中均保留了D2H3-2模型在a, b, c, d工况下的计算误差作为参照. 曲线上标注值为该数值模拟的$ \Delta u $值. 图中a, b, c, d工况分别用实线, 虚线, 虚点线及点线标注, 而离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s

    Figure 6.  The numerical performance of Taylor-Green vortex simulations vs the value of $ \Delta u $. Panel a plots the numerical errors of model D2H4 under case a, b, c and model D2H5 under case b, c, d, which possess close value of $ \Delta u $ respectively. Panel b plots the numerical errors of simulations with a value of $ \Delta u $ around $ -0.30 $. The numbers labeled on the curves are their values of $ \Delta u $. The simulation configurations are denoted with line style, in which solid, dashed, dash-dotted and dotted line indicates case a, b, c and d respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. In all panels, the results of model D2H3-2 with four designed configurations are plotted for a reference. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed, dash-dotted and dotted line indicates case a, b, c and d in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

    图 7  Taylor-Green涡数值模拟误差与初始$ \kappa $值相关性分析. 图中对比了在$ \Delta u $值为负情况下, 初始$ \kappa $值偏离$ 1 $幅值可忽略的算例计算误差, 包括模型D2H3-1工况c, 模型D2H3-2工况d, 模型D2H4工况b, 模型D2H5工况c, d和模型D2H6工况b, c. 图中均保留了D2H3-2模型在a, b, c, d工况下的计算误差作为参照. 曲线上标注值为该数值模拟的$ \Delta u $值. 图中a, b, c, d工况分别用实线, 虚线, 虚点线及点线标注, 而离散模型则用颜色和符号标注. 为方便表示, 图中横坐标单位为0.173286795 s

    Figure 7.  The numerical performance of Taylor-Green vortex simulations vs initial $ \kappa $. The error evolutions of numerical simulations with negative values of $ \Delta u $ but negligible departure of initial $ \kappa $ from $ 1 $ are plotted, including model D2H3-1 under case c, model D2H3-2 under case d, model D2H4 under case b, model D2H5 under case c,d and model D2H6 under case b, c. The results of model D2H3-2 with four designed configurations are also plotted for a reference. The designed configurations of Taylor-Green vortex are denoted with line styles, in which the solid, dashed and dash-dotted line indicates case a, b and c in Table 3 respectively, meanwhile the lattice Boltzmann models are labeled with line colors and markers. For the sake of rendering, the unit of the time axis is scaled as 0.173286795 s.

    表 1  前五阶物理Hermite多项式

    Table 1.  The first five physicists′ Hermite polynomials.

    Hermite 多项式阶数
    The degree of Hermite polynomial
    多项式表达式
    Hermite polynomial
    0 ${ { H}_0}\left( x \right) = 1$
    1 ${ { H}_1}\left( x \right) = 2x$
    2 ${ { H}_2}\left( x \right) = 4x^2-2$
    3 ${ { H}_3}\left( x \right) = 8x^3-12x$
    4 ${ {H}_4}\left( x \right) = 16x^4-48x^2+12$
    DownLoad: CSV

    表 2  格子Boltzmann离散模型. 为方便展示, 表格罗列的是用于对应构造高阶模型的一维模型参数. 实际计算中需根据(9)式和(10)式对表格中的参数进行张量构造. 另外需注意的是表格中罗列的是网格常数c, 其与格子声速$ c_{\rm s} $存在如下换算关系$ c = 1/\sqrt{2}c_{\rm s} $

    Table 2.  Discrete model description of lattice Boltzmann method. For the seeking of space saving, the parameters illustrated in this table are of the corresponding unidimensional models. In numerical implementations, they require tensor product illuminated in Eq. (9) and Eq. (10). It should be noted that the table lists the lattice constant c instead of the lattice sonic speed $ c_{\rm s} $, which can be expressed as $ c = 1/\sqrt{2}c_{\rm s} $.

    Model name Discrete velocityset {$v_\alpha$}  Lattice constant c Weights {${w}_\alpha$}
    D2H2 {$0, \pm 1$} $1.2247\times10^0$ {$6.6667\times10^{-1}$, 1.6667$\times10^{-1}$}
    D2H3-1 {$0, \pm 1, \pm 3$} $5.5343\times10^{-1}$ {$7.4464\times10^{-2}$, $4.1859\times10^{-1}$, $4.4182\times10^{-2}$}
    D2H3-2 {$0, \pm 2, \pm 5$} $3.4420\times10^{-1}$ {$3.1044\times10^{-1}$, $3.0997\times10^{-1}$, $3.4812\times10^{-2}$}
    D2H4 {$0, \pm 1, \pm 2, \pm 3$} $8.4639\times10^{-1}$ {$4.7667\times10^{-1}$, $2.3391\times10^{-1}$, $2.6938\times10^{-2}$, $8.1213\times10^{-4}$}
    D2H5 {$0, \pm 1, \pm 2, \pm 3, \pm 5$} $4.7940\times10^{-1}$ {$1.6724\times10^{-1}$, $3.0315\times10^{-1}$, $5.3303\times10^{-2}$, $5.7922\times10^{-2}$, $2.0013\times10^{-3}$}
    D2H6 {$0, \pm 1, \pm 2, \pm 3, \pm 4, \pm 5$} $6.8590\times10^{-1}$ {$3.8694\times10^{-1}$, $2.4178\times10^{-1}$, $5.8922\times10^{-2}$, $5.6153\times10^{-3}$, $2.0652\times10^{-4}$, $3.2745\times10^{-6}$}
    DownLoad: CSV

    表 3  算例设计

    Table 3.  Case design for lattice Boltzmann method

    Case name Step time/s $u_{L{\rm B} ,0}$
    a $1.570796\times10^{-2}$ 0.5
    b $3.141593\times10^{-2}$ 1.0
    c $4.712389\times10^{-2}$ 1.5
    d $6.283185\times10^{-2}$ 2.0
    DownLoad: CSV
    Baidu
  • [1]

    Ammar S, Pernaudat G, Trepanier J 2017 J. Comput. Phys. 343 73Google Scholar

    [2]

    Li Q, Luo K H, Kang Q J, He Y L, Chen Q, Liu Q 2016 Prog. Energ. Combust. 52 62Google Scholar

    [3]

    Shan X W 2010 Phys. Rev. E 81 045701Google Scholar

    [4]

    Shan X W, Chen H D 1993 Phys. Rev. E 47 1815Google Scholar

    [5]

    Sheikholeslami M, Gorji-Bandpy M, Ganji D D 2014 Powder Technol. 254 82Google Scholar

    [6]

    Zhang X L, Xiao L Z, Shan X W, Guo L 2014 Sci. Rep. 4 4843

    [7]

    Shan X W 2011 IMA J. Appl. Math. 76 650Google Scholar

    [8]

    Kang Q J, Zhang D X, Chen S Y, He X Y 2002 Phys. Rev. E 65 036318Google Scholar

    [9]

    Mountrakis L, Lorenz E, Malaspinas O, Alowayyed S, Chopard B, Hoekstra A G 2015 J. Comput. Sci. 9 45Google Scholar

    [10]

    Ladd A J C, Verberg R 2001 J. Stat. Phys. 104 1191Google Scholar

    [11]

    Chopard B, Ouared R, Deutsch A, Hatzikirou H, Wolf-Gladrow D 2010 Acta Biotheor. 58 329Google Scholar

    [12]

    Hirabayashi M, Ohta M, Rjfenacht D A, Chopard B 2003 Phys. Rev. E 68 021918Google Scholar

    [13]

    Latt J, Chopard B 2007 Int. J. Mod. Phys. C 18 619Google Scholar

    [14]

    Yu H D, Luo L S, Girimaji S S 2006 Comput. Fluids 35 957Google Scholar

    [15]

    Yu H D, Girimaji S S, Luo L S 2005 J. Comput. Phys. 209 599Google Scholar

    [16]

    He X Y, Chen S Y, Doolen G D 1998 J. Comput. Phys. 146 282Google Scholar

    [17]

    He X Y, Luo L S 1997 Phys. Rev. E 56 6811Google Scholar

    [18]

    Ye H F, Kuang B, Yang Y H 2019 Chin. Phys. B 28 014701Google Scholar

    [19]

    Dellar P J 2014 J. Comput. Phys. 259 270Google Scholar

    [20]

    Nie X B, Shan X W, Chen H D 2008 EPL-Europhys. Lett. 81 34005Google Scholar

    [21]

    Hazi G, Kavran P 2006 J. Phys. A 39 3127Google Scholar

    [22]

    Shan X W 2016 J. Comput. Sci. 17 475Google Scholar

    [23]

    Shan X W 2010 Phys. Rev. E 81 036702Google Scholar

    [24]

    Shim J W 2013 Phys. Rev. E 87 013312Google Scholar

    [25]

    Ye H F, Gan Z C, Kuang B, Yang Y H 2019 Chin. Phys. B 28 054702Google Scholar

    [26]

    Shan X W, He X Y 1998 Phys. Rev. Lett. 80 65Google Scholar

  • [1] Sun Jia-Kun, Lin Chuan-Dong, Su Xian-Li, Tan Zhi-Cheng, Chen Ya-Lou, Ming Ping-Jian. Solution of the discrete Boltzmann equation: Based on the finite volume method. Acta Physica Sinica, 2024, 73(11): 110504. doi: 10.7498/aps.73.20231984
    [2] Liu Cheng, Liang Hong. Axisymmetric lattice Boltzmann model for three-phase fluids and its application to the Rayleigh-Plateau instability. Acta Physica Sinica, 2023, 72(4): 044701. doi: 10.7498/aps.72.20221967
    [3] Lou Qin, Huang Yi-Fan, Li Ling. Lattice Boltzmann model of gas-liquid two-phase flow of incomprssible power-law fluid and its application in the displacement problem of porous media. Acta Physica Sinica, 2019, 68(21): 214702. doi: 10.7498/aps.68.20190873
    [4] Hu Jia-Yi, Zhang Wen-Huan, Chai Zhen-Hua, Shi Bao-Chang, Wang Yi-Hang. Three-dimensional 12-velocity multiple-relaxation-time lattice Boltzmann model of incompressible flows. Acta Physica Sinica, 2019, 68(23): 234701. doi: 10.7498/aps.68.20190984
    [5] Li Yang, Su Ting, Liang Hong, Xu Jiang-Rong. Phase field lattice Boltzmann model for two-phase flow coupled with additional interfacial force. Acta Physica Sinica, 2018, 67(22): 224701. doi: 10.7498/aps.67.20181230
    [6] Liu Gao-Jie, Guo Zhao-Li, Shi Bao-Chang. A coupled lattice Boltzmann model for fluid flow and diffusion in a porous medium. Acta Physica Sinica, 2016, 65(1): 014702. doi: 10.7498/aps.65.014702
    [7] Huang Wei-Chao, Liu Ding, Jiao Shang-Bin, Zhang Ni. Thermo-fluid coupling of unsteady flow in Czochralski crystal growth. Acta Physica Sinica, 2015, 64(20): 208102. doi: 10.7498/aps.64.208102
    [8] Xie Wen-Jun, Teng Peng-Fei. Study of acoustic levitation by lattice Boltzmann method. Acta Physica Sinica, 2014, 63(16): 164301. doi: 10.7498/aps.63.164301
    [9] Liu Qiu-Zu, Kou Zi-Ming, Jia Yue-Mei, Wu Juan, Han Zhen-Nan, Zhang Qian-Qian. Wettability alteration simulation of modified hydrophobic solid surface by lattice Boltzmann method. Acta Physica Sinica, 2014, 63(10): 104701. doi: 10.7498/aps.63.104701
    [10] Zeng Jian-Bang, Li Long-Jian, Jiang Fang-Ming. Numerical investigation of bubble nucleation process using the lattice Boltzmann method. Acta Physica Sinica, 2013, 62(17): 176401. doi: 10.7498/aps.62.176401
    [11] He Yu-Bo, Lin Xiao-Yan, Dong Xiao-Liang. Use of lattice Boltzmann method to simulate 2-D partial differential equation. Acta Physica Sinica, 2013, 62(19): 194701. doi: 10.7498/aps.62.194701
    [12] Shi Zi-Yuan, Hu Guo-Hui, Zhou Zhe-Wei. Lattice Boltzmann simulation of droplet motion driven by gradient of wettability. Acta Physica Sinica, 2010, 59(4): 2595-2600. doi: 10.7498/aps.59.2595
    [13] Sun Qi-Cheng, Wang Guang-Qian. Force distribution in static granular matter in two dimensions. Acta Physica Sinica, 2008, 57(8): 4667-4674. doi: 10.7498/aps.57.4667
    [14] Feng Shi-De, Zhong Lin-Hao, Gao Shou-Ting, Dong Ping. Equilibrium distribution boundary condition in lattice Boltzmann model and numerical simulation of Darcy-Forcheimer drag for fluid flow across a square cylinder array. Acta Physica Sinica, 2007, 56(3): 1238-1244. doi: 10.7498/aps.56.1238
    [15] Zhao Ying, Ji Zhong-Zhen, Feng Tao. Simulation of thermal convection in a vertical slot using the lattice Boltzmann model. Acta Physica Sinica, 2004, 53(3): 671-675. doi: 10.7498/aps.53.671
    [16] LI HUA-BING, HUANG PING-HUA, LIU MU-REN, KONG LING-JIANG. SIMULATION OF THE MKDV EQUATION WITH LATTICE BOLTZMANN METHOD. Acta Physica Sinica, 2001, 50(5): 837-840. doi: 10.7498/aps.50.837
    [17] FENG SHI-DE, ZHANG QIONG, REN RONG-CAI. SIMULATION OF A FLOW FIELD WITH NONUNIFORM TEMPERATURE BY USING LATTICE BOLTZMANN EQUATION MODEL. Acta Physica Sinica, 2001, 50(7): 1207-1212. doi: 10.7498/aps.50.1207
    [18] FENG SHI-DE, MICHIHISA TSUTAHARA. SIMULATION OF SHOCK WAVES USING A LATTICE BOLTZMANN EQUATION MODEL. Acta Physica Sinica, 2001, 50(6): 1006-1010. doi: 10.7498/aps.50.1006
    [19] Yu Hui-Dan, Zhao Kai-Hua. . Acta Physica Sinica, 2000, 49(4): 816-818. doi: 10.7498/aps.49.816
    [20] YU HUI-DAN, ZHAO KAI-HUA. LATTICE BOLTZMANN MODEL FOR COMPRESSIBLE FLOW SIMULATION. Acta Physica Sinica, 1999, 48(8): 1470-1476. doi: 10.7498/aps.48.1470
Metrics
  • Abstract views:  6323
  • PDF Downloads:  30
  • Cited By: 0
Publishing process
  • Received Date:  26 April 2019
  • Accepted Date:  06 August 2019
  • Available Online:  01 October 2019
  • Published Online:  20 October 2019

/

返回文章
返回
Baidu
map