A.L. Kupershtokh, D.A. Medvedev, & D.I. Karpov (2009). On equations of state in a lattice Boltzmann method. Computers & Mathematics with Applications, 58(5), 965-974. DOI:10.1016/j.camwa.2009.02.024. ↩︎↩︎↩︎
Q. Li, K.H. Luo, Q.J. Kang, Y.L. He, Q. Chen, & Q. Liu (2016). Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Progress in Energy and Combustion Science, 52, 62-105. DOI:10.1016/j.pecs.2015.10.001. ↩︎
Li, Q., Zhou, P., & Yan, H. (2016). Revised Chapman-Enskog analysis for a class of forcing schemes in the lattice Boltzmann method. Physical Review E, 94, 043313. DOI:10.1103/PhysRevE.94.043313. ↩︎↩︎↩︎↩︎
Wagner, A. (2006). Thermodynamic consistency of liquid-gas lattice Boltzmann simulations. Physical Review E, 74, 056703. DOI: 10.1103/PhysRevE.74.056703. ↩︎
""" Function to convert DEM simulation data into coarse graining field. NOTE: This function is used for graunlar system which consists of single type of particles. """ import typing as _T
import numpy as np import numba as nb import numpy.typing as npt
from tqdm import tqdm from scipy.special import erf from scipy.integrate import simpson
@nb.jit(nopython=True) defCG_Gaussian(r: npt.NDArray[np.double], w: float, V_w: float) -> float: """Gaussian Coarse Graining function. Parameters ==== r: numpy.NDArray[double] Relative position vector, whose shape is `(3,)`. w: float Course grain width `w`. V_w: float The coefficient of Gaussian Coarse Graining function. Return === W: float Average weight. """ c = 3 * w # Coarse graining cut-off width c_2 = c**2
# Norm of vector `r` r_norm: float = r[0]**2 + r[1]**2 + r[2]**2 # define weight W = 0.0 if r_norm < c_2: W = np.exp(-r_norm / (2 * w**2)) / V_w return W
@nb.jit("f8[:, :](f8, f8, f8[:], f8[:])", nopython=True, parallel=True, nogil=True) defSIGMA_q_k(mass_i: float, W_i: float, vel_i: npt.NDArray[np.double], U: npt.NDArray[np.double]) -> npt.NDArray[np.double]: """Calculate part of `sigma_qk` contributed by ball `i` Parameters --- mass_i: float Mass of ball i. W_i: float Weight coefficient of ball i. vel_i: numpy.NDArray[double] Velocity vector of ball i, whose shape is (3,). U: numpy.NDArray[double] Coarse grained velocity vector, whose shape is (3,). Return --- res: numpy.NDArray[double] 3x3 Matrix of `sigma_qk`caused by ball i. """ res = np.zeros((3,3), dtype=np.double) for i in nb.prange(3): for j in nb.prange(3): res[i, j] += mass_i * W_i * \ (vel_i[i] - U[i]) * (vel_i[j] - U[j]) return res
@nb.jit("f8[:, :](f8[:], f8[:], f8)", nopython=True, parallel=True, nogil=True) defSIGMA_q_c(F_ij: npt.NDArray[np.double], l_ij: npt.NDArray[np.double], W_ij: float) -> npt.NDArray[np.double]: """Calculate part of `sigma_qk` contributed by contact from i to j. Parameters --- F_ij: numpy.NDArray[double] Contact force vector, pointing from i to j. l_ij: numpy.NDArray[double] Contact branch vector, pointing from i to j. W_ij: float Integrated weight coefficient, from i to j. Return --- res: numpy.NDArray[double] 3x3 Matrix of `sigma_qc`caused by contact from i to j. """ res = np.zeros((3,3), dtype=np.double) for i in nb.prange(3): for j in nb.prange(3): res[i,j] = F_ij[i] * l_ij[j] * W_ij return res
classCG_mono_3D: def__init__(self, rho: float, pos: npt.NDArray[np.double], vel: npt.NDArray[np.double], radius: npt.NDArray[np.double], ball_id: npt.NDArray[np.int_], contact_force: npt.NDArray[np.double], contact_branch: npt.NDArray[np.double], contact_end_ball: npt.NDArray[np.double], in_flow: _T.Optional[npt.NDArray[np.bool_]] = None, dp: _T.Optional[float] = None, *, progress_bar: bool = False ) -> None: """ Parameters ---- rho: float Density of all particles. pos: numpy.NDArray[double] Position vector, whose shape is `[particle_num, 3]`. vel: numpy.NDArray[double] Velocity vector, whose shape is `[particle_num, 3]`. radius: numpy.NDArray[double] Radius list, whose shape is `(particle_num,)`. ball_id: numpy.NDArray[int] ID number of all balls. contact_force: numpy.NDArray[double] Contact force list, whose shape is `[contact_num, 3]`. contact_branch: numpy.NDArray[double] List of contact branch, whose shape is `[contact_num, 3]`. Each line `i` means the distance vector from `contact_end_ball[i, 0]` to `contact_end_ball[i, 1]`. contact_end_ball: numpy.NDArray[double] List of contacted balls tuple, whose shape is `[contact_num, 2]`. in_flow: numpy.NDArray[np.bool_], optional Mask array of ball-in-flow, whose shape is `(particle_num,)`. dp: float, optional Mean diameter of all particles. progress_bar: bool, optional, keyword argument only. Switch to control whether show the process. Default is False. Note ---- 1. Take the `i`-th ball for example, its ball ID is `ball_id[i]`. And its radius and velocity are `radius[i]` and `vel[i,:]`, respectively. Note that `ball_id[i]` is not equal to `i`. 2. For the `j`-th contact, The ID of two contacted balls are `b1 = contact_end_ball[j, 0]`, and `b2 = contact_end_ball[j ,1]`. """ self.particle_num = pos.shape[0] self.contact_num = contact_force.shape[0]
self._id = ball_id # Ball id self.pos = pos # Position vector self.vel = vel # Velocity vector self.radius = radius # Ball radius
# Ball density # |=> NOTE: Since this function is used for mono-disperse system, # The density of particles are constant. self.rho = rho
# => Ball-ID position array. # For example, `self._c_end_id[j,0]` corresponds to the first end # ball of `j`-th contact (whose ID is `self.c_ends[j,0]`). # And we have `self.c_ends[j,0] = self._id[self._c_end_id[j,0]]`. self._c_end_id = [ [np.where(ball_id == contact_end_ball[i,0])[0][0], np.where(ball_id == contact_end_ball[i,1])[0][0]] for i inrange(self.contact_num) ] self._c_end_id = np.array(self._c_end_id, dtype=np.intc)
# Unit normal of each contact self.contact_norm = []
# Real contact branch vector, pointing from contact point to balls' center self.branch_vector_list = self._get_branch_vector()
# => STATUS of domain config self._domain_init = False
# => Max search width self._MAX_width = 7
# => Option to show progress bar self._use_tqdm = progress_bar pass
defset_domain(self, xlim: _T.Tuple[float, float], ylim: _T.Tuple[float, float], zlim: _T.Tuple[float, float], periodic_flag: _T.Tuple[bool, bool, bool]): """Function to set cuboid domain of coarse graining. Parameters ---- xlim: (float, float) Range of x-coordinate, contains with `(xMin, xMax)`. ylim: (float, float) Range of y-coordinate, contains with `(yMin, yMax)`. zlim: (float, float) Range of z-coordinate, contains with `(zMin, zMax)`. periodic_flag: (bool, bool, bool) Flag tuple to determinate which axis is periodic. """ xlim = (min(*xlim), max(*xlim)) ylim = (min(*ylim), max(*ylim)) zlim = (min(*zlim), max(*zlim)) self._lims = np.array([xlim, ylim, zlim])
defget_CG_sigma_qc(self, R: npt.NDArray[np.double], CG_width: float, V_w: float) -> npt.NDArray[np.double]: """ Get stress tensor contributed by contact. Parameter --- R: numpy.NDArray[double] Center of coarse graining node. CG_width: float Width of coarse graining. V_w: float The coefficient of Gaussian Coarse Graining function. Return --- sigma_qc: numpy.NDArray[double] Stress tensor contributed by contact. Note --- $\sigma_{\\alpha \\beta}^{q(\mathrm{ck})} = \sum_{i \in q} \sum\limits_{j \in \bar{Q}}^{j \neq i} F_{ij,\\alpha} l_{ij,\\beta} \int_{0}^{1} W(\\boldsymbol{r} - \\boldsymbol{r}_i + s \\boldsymbol{l}_{ij}) \mathrm{d}s$ where: * $v_{i \\alpha}^{'} = u_{\\alpha}^{q}(\\boldsymbol{r}, t) - u_{i,\\alpha}$ is the $\\alpha$ component of the fluctuation velocity of particle i. * $\\bar{Q}$ is the union of all particles classes, $q$. * $\\boldsymbol{l}_{ij}$ is as the branch vector between two particles, i and j. (from i to j) * $\\boldsymbol{F}_{ij}$ is contact force vector between two particles, i and j. (from i to j) """ MAX_width = self._MAX_width * CG_width
# Result array sigma_qc = np.zeros((3,3), dtype=np.double)
# ========================== # # === Auxiliary function === # # ========================== # def_W_ij_4_sigma_qc(self, dR, l_ij, CG_width, V_w) -> float: """ Parameters --- dR: numpy.NDArray[double] Relative position vector, whose shape is `(3,)`. l_ij: numpy.NDArray[double] Contact branch vector, whose shape is `(3,)`. CG_width: float Width of coarse graining. V_w: float The coefficient of Gaussian Coarse Graining function. Return ---- W_ij: float Integrated weight coefficient used in the calculation of `sigma_qc`. """ _func = lambda s: CG_Gaussian(dR + s * l_ij, CG_width, V_w) _x = np.linspace(0, 1, 201) _y = np.array([_func(s) for s in _x]) W_ij = simpson(_y, x=_x) return W_ij
def_get_distance_vec(self, R: npt.NDArray[np.double], r_i: npt.NDArray[np.double]) -> npt.NDArray[np.double]: """Get distance between center of coarse graining node `R` and loaction of particle's center `r_i`. Parameter --- R: numpy.NDArray[double] Center of coarse graining node. r_i: numpy.NDArray[double] Loaction of particle's center. Return --- dR: numpy.NDArray[double] Distance between `R` and `r_i` (i.e. `R - r_i`), with periodic boundary considered. """ dR = get_distance(R, r_i, self._lims, self._hasPeriodic, self._isPeriodic) return dR
@property defprogress_bar(self): """Bool flag to control whether show the process by tqdm.""" returnself._use_tqdm
@progress_bar.setter defprogress_bar(self, value: bool): self._use_tqdm = bool(value) @property defsearch_width_coef(self): """Coefficient of search radius, and the search radius is `self.search_width_coef * CG_width`. """ returnself._MAX_width @search_width_coef.setter defsearch_width_coef(self, value: float): if value <= 0: raise ValueError("`search_width_coef` should be positive") elif value <= 3: raise ValueError("`search_width_coef` should be larger than cut off width coefficient (3.0)") else: self._MAX_width = value defget_CG_factors(self, w: float) -> float: """Calculate the coefficient of Gaussian Coarse Graining function. Parameter ---- w: float Course grain width `w`. Return ---- V_w: float The coefficient of Gaussian Coarse Graining function. """ c = 3 * w V_w = np.sqrt(8) * np.pi**1.5 * w**3 * erf(c/w/np.sqrt(2)) - \ 4 * c * w**2 * np.pi * np.exp(-c**2 / w**2 / 2) return V_w
def_get_branch_vector(self) -> list[list[np.ndarray]]: """Calculate the real contact branch vector, because `self.c_branch` only means the distance of the 2 balls' spherical centers on the same contact. """ branch_vector_list = [[], []] self.contact_norm = []
for i inrange(self.contact_num): id_i = self._c_end_id[i, 0] id_j = self._c_end_id[i, 1] r_i = self.radius[id_i] r_j = self.radius[id_j]
# Distance of the 2 balls' spherical centers (from i to j) l_ij = self.c_branch[i, :] l_norm = np.linalg.norm(l_ij) N = l_ij / l_norm
离散单元法(Discrete Element Method,简称DEM)是颗粒流数值模拟的重要数值方法,被用于各种各样的颗粒流定量研究中。颗粒流的粗粒化(Coarse Graining),便是将颗粒流流场的大量数据转化为连续场数据的数值方法。
在开始介绍之前,先简要列出主要的参考资料:
Breard, E. C. P., Dufek, J., Fullard, L., & Carrara, A. (2020). The Basal Friction Coefficient of Granular Flows With and Without Excess Pore Pressure: Implications for Pyroclastic Density Currents, Water‐Rich Debris Flows, and Rock and Submarine Avalanches. Journal of Geophysical Research: Solid Earth, 125(12), e2020JB020203. DOI:10.1029/2020JB020203.
Lacaze, L., & Kerswell, R. R. (2009). Axisymmetric Granular Collapse: A Transient 3D Flow Test of Viscoplasticity. Physical Review Letters, 102(10), 108305. DOI:10.1103/PhysRevLett.102.108305.
Weinhart, T., Labra, C., Luding, S., & Ooi, J. Y. (2016). Influence of coarse‐graining parameters on the analysis of DEM simulations of silo flow. Powder Technology, 293, 138–148. DOI:10.1016/j.powtec.2015.11.052
单种颗粒的粗粒化
为简单起见,先从仅存在单种球形颗粒的颗粒流开始说起,并且只考虑获取空间内一点 r 在 t 时刻的粗粒化数据。
粗粒化函数
在这里,引入高斯函数
W(r)={Vw1exp(−2w∣r∣2),0,∣r∣<cotherwise
其中: w 为粗粒化宽度; c=3w 为粗粒化截断长度; Vw 是确保其密度积分等于总质量的系数
Vw=22⋅w3π23erf(2wc2)−4cw2πexp(2w2−c2)
Weinhart 等建议 w 的取值范围应在 [0.75dmean,1.25dmean] 之间,其中 dmean 为平均颗粒直径。
密度和动量的计算
在 t 时刻,单种颗粒的集合为 q 。对于其中的第 i 个球,其位置矢量为 ri ,质量为 mi ,速度为 ui。
Chen, Z., Shu, C., Wang, Y., Yang, L. M., & Tan, D. (2017). A Simplified Lattice Boltzmann Method without Evolution of Distribution Function. Advances in Applied Mathematics and Mechanics, 9(1), 1–22. DOI:10.4208/aamm.OA-2016-0029 ↩︎
Z. Chen, C. Shu, D. Tan; Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows. Physics of Fluids. 1 May 2017; 29 (5): 053601. DOI:10.1063/1.4983339 ↩︎↩︎↩︎
Chen Z, Shu C, Tan D, Wu C. On improvements of simplified and highly stable lattice Boltzmann method: Formulations, boundary treatment, and stability analysis. Int J Numer Meth Fluids. 2018; 87: 161–179. DOI:10.1002/fld.4485 ↩︎
Chen, Z., & Shu, C. (2020). Simplified and Highly Stable Lattice Boltzmann Method. WORLD SCIENTIFIC. DOI:10.1142/12047 ↩︎
C. Shu, Y. Wang, C. Teo, and J. Wu. Development of lattice Boltzmann flux solver for simulation of incompressible flows. Advances In Applied Mathematics And Mechanics. 6, 436 (2014). ↩︎
GRAD H. Note on N-Dimensional Hermite Polynomials[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 325–330. http://doi.org/10.1002/cpa.3160020402↩︎
GRAD H. On the Kinetic Theory of Rarefied Gases[J]. Communications on Pure and Applied Mathematics, 1949, 2(4): 331–407. http://doi.org/10.1002/cpa.3160020403↩︎
SHAN X, YUAN X-F, CHEN H. Kinetic Theory Representation of Hydrodynamics: A Way beyond the Navier–Stokes Equation[J]. Journal of Fluid Mechanics, 2006, 550(1): 413. http://doi.org/10.1017/S0022112005008153↩︎↩︎
CAO W. Investigation of the Applicability of the Lattice Boltzmann Method to Free-Surface Hydrodynamic Problems in Marine Engineering[D/OL]. Laboratoire de recherche en Hydrodynamique, Énergétique et Environnement Atmosphérique (LHEEA): École centrale de Nantes, 2019. https://tel.archives-ouvertes.fr/tel-02383174↩︎