Wind-blown particulate transport: A review of computational fluid dynamics models

Andrea Lo Giudice; Roberto Nuca; Luigi Preziosi; Nicolas Coste; Andrea Lo Giudice; Roberto Nuca; Luigi Preziosi; Nicolas Coste

doi:10.3934/mine.2019.3.508

Mathematics in Engineering

2019, Volume 1, Issue 3: 508-547. doi: 10.3934/mine.2019.3.508

Previous Article Next Article

Review

Wind-blown particulate transport: A review of computational fluid dynamics models

1.
Optiflow Company, 160 Chemin de la Madrague-Ville, Marseille, 13015, France
2.
Department of Mathematical Sciences Giuseppe Luigi Lagrange, Politecnico di Torino, Corso Duca degli Abruzzi 24, Torino, 10129, Italy

Received: 26 September 2018 Accepted: 26 April 2019 Published: 04 July 2019

The transport of particulate by wind constitutes a relevant phenomenon in environmental sciences and civil engineering, because erosion, transport and deposition of particulate can cause serious problems to human infrastructures. From a mathematical point of view, modeling procedure for this phenomenon requires handling the interaction between different constituents, the transfer of a constituent from the air to the ground and viceversa, and consequently the ground-surface interaction and evolution. Several approaches have been proposed in the literature, according to the specific particulate or application. We here review these contributions focusing in particular on the behavior of sand and snow, which almost share the same mathematical modeling issues, and point out existing links and analogies with wind driven rain. The final aim is then to classify and analyze the different mathematical and computational models in order to facilitate a comparison among them. A first classification of the proposed models can be done distinguishing whether the dispersed phase is treated using a continuous or a particle-based approach, a second one on the basis of the type of equations solved to obtain particulate density and velocity, a third one on the basis of the interaction model between the suspended particles and the transporting fluid.

Keywords:

Citation: Andrea Lo Giudice, Roberto Nuca, Luigi Preziosi, Nicolas Coste. Wind-blown particulate transport: A review of computational fluid dynamics models[J]. Mathematics in Engineering, 2019, 1(3): 508-547. doi: 10.3934/mine.2019.3.508

Related Papers:

[1]	Tianyuan Xu, Shanming Ji, Chunhua Jin, Ming Mei, Jingxue Yin . EARLY AND LATE STAGE PROFILES FOR A CHEMOTAXIS MODEL WITH DENSITY-DEPENDENT JUMP PROBABILITY. Mathematical Biosciences and Engineering, 2018, 15(6): 1345-1385. doi: 10.3934/mbe.2018062
[2]	Wenjie Zhang, Lu Xu, Qiao Xin . Global boundedness of a higher-dimensional chemotaxis system on alopecia areata. Mathematical Biosciences and Engineering, 2023, 20(5): 7922-7942. doi: 10.3934/mbe.2023343
[3]	Sunwoo Hwang, Seongwon Lee, Hyung Ju Hwang . Neural network approach to data-driven estimation of chemotactic sensitivity in the Keller-Segel model. Mathematical Biosciences and Engineering, 2021, 18(6): 8524-8534. doi: 10.3934/mbe.2021421
[4]	Qianhong Zhang, Fubiao Lin, Xiaoying Zhong . On discrete time Beverton-Holt population model with fuzzy environment. Mathematical Biosciences and Engineering, 2019, 16(3): 1471-1488. doi: 10.3934/mbe.2019071
[5]	Chichia Chiu, Jui-Ling Yu . An optimal adaptive time-stepping scheme for solving reaction-diffusion-chemotaxis systems. Mathematical Biosciences and Engineering, 2007, 4(2): 187-203. doi: 10.3934/mbe.2007.4.187
[6]	Xu Song, Jingyu Li . Asymptotic stability of spiky steady states for a singular chemotaxis model with signal-suppressed motility. Mathematical Biosciences and Engineering, 2022, 19(12): 13988-14028. doi: 10.3934/mbe.2022652
[7]	Tingting Yu, Sanling Yuan . Dynamics of a stochastic turbidostat model with sampled and delayed measurements. Mathematical Biosciences and Engineering, 2023, 20(4): 6215-6236. doi: 10.3934/mbe.2023268
[8]	Lin Zhang, Yongbin Ge, Zhi Wang . Positivity-preserving high-order compact difference method for the Keller-Segel chemotaxis model. Mathematical Biosciences and Engineering, 2022, 19(7): 6764-6794. doi: 10.3934/mbe.2022319
[9]	Changwook Yoon, Sewoong Kim, Hyung Ju Hwang . Global well-posedness and pattern formations of the immune system induced by chemotaxis. Mathematical Biosciences and Engineering, 2020, 17(4): 3426-3449. doi: 10.3934/mbe.2020194
[10]	Marcin Choiński, Mariusz Bodzioch, Urszula Foryś . A non-standard discretized SIS model of epidemics. Mathematical Biosciences and Engineering, 2022, 19(1): 115-133. doi: 10.3934/mbe.2022006

Abstract

1. Introduction

Nanomedicine knows an increasing interest due to its estimated importance in efficient medication delivery ^[1]. Crucial medical topics deal with its application in several fields, such as overcoming the blood-brain barrier ^[2], anti-diabetic medication ^[3] and uptake in endothelial cells ^[4]. Nanomedicine deals with the novel chemical, physical or biological properties of medication at nanometric scale, as well as the application of nanotechnology aiming at making a medical diagnosis, treating or preventing the cause of illness. Therefore, nanomedicines have the ability, due to their particular size or surface functionality, to demonstrate physicochemical and/or biological properties that make them interesting for medication delivery, controlled medication release and/or improved medication transport through biological barriers, such as cell membranes, which would otherwise be much more difficult ^[5]. One of the key parameters for successful medication delivery is controlling or improving the permeation of medication through the destined place (e.g. vascular system, cell) in the human body ^[6,7,8]. When speaking about cell permeation and the body uptake efficiency of drugs into the body that it entails, many works take the apparent cell permeability, P_app, through monolayers of human epithelial cells as a widely accepted standard ^[9,10,11]. Values of P_app (often presented in medicine as a volumetric flux per unity of area) of such cells are also often used for the estimation of the blood-brain barrier ^[12,13]. It is useful to understand the apparent cell permeation as a series of parallel and serial resistances (from layer to another). This vision is based on the concept of solubility diffusion ^[14]. The permeability is seen here, therefore, not as an imposed force, but is rather passive in a way that only spontaneous diffusion governs the permeation process. This diffusion process is determined by a concentration gradient. The concentration gradient is classically described by Fick's law or, in case of sizes (within the membranes separating the monolayers) at nanoscale, by a more general evolution equation, describing the extended diffusion law, as will be discussed later. Active transport, as could be caused by transporter proteins, is not considered here. Although, in reality, chemicals/drugs could be transported by proteins (and thus have an effect on P_app), it is not taken into account here. It is well admitted that generally three pathways are possible through a cellular monolayer, i.e. in between the cells (the paracellular way), through the cell (considered here) and along the membrane solely (the lateral way) ^[15,16]. In this case, parallel resistances of the three pathways would apply. Here, we only consider the pathway through the cell. Following the pH partition hypothesis ^[17], it can be assumed that only the neutral form of the drug chemicals can pass the membrane resistances. Therefore, chemicals that are ionized should undergo an acid-base reaction in between the water layer and the membrane before they are able to pass the membrane as a neutral species. According to the aforementioned pH partition hypothesis, such a reaction can be considered to occur instantaneously, meaning that no kinetic hinderness (and therefore no additional resistance) needs to be considered.

It is often suggested that the diffusion through membranes are of non-Fickian nature. We propose that this is either because of non-local effects (the size of the membranes are of the same order of magnitude as the drug molecules or not much bigger) or because the relaxation time becomes important. As controlling and enhancing drug delivery is an important issue in the field of nanomedicine, next to in vitro experiments, modelling is crucial to understand the underlying mechanisms. In nanomedicine, the inertia of the diffusing particles is relevant in addition to the small size of the channels through which the medication passes, being in some cases of the same order of magnitude as the mean free path of the particles. Furthermore, purely microscopic models are often difficult to implement directly to real-life applications. Therefore, further progress in this field can be enhanced by models that are easily applicable, whilst keeping a qualitative efficient and understandable way to link to experiments in a multi-disciplinary way. Therefore, we propose here a model based on thermodynamic principles, where we extend them by incorporating effects, in a rigourous way, that account for the relaxation and size of the drug molecules with respect to the channels through which they diffuse. With the purpose to demonstrate the applicability of our new theory, we treat two cases: drug diffusion through a cell membrane and drug release out of a scaffold. For the first case, let us imagine the existence of two compartments of well-mixed water that are separated by a barrier, i.e. a membrane: the first containing the drug molecules as the donor and the second represents the inner cell or the receptor. The purpose is to follow the evolution of the mass fraction profile through the membrane at the beginning of the diffusion process, where we assume that the drug molecules (in that beginning process) concentration is that high that the mass fraction at the inlet of the membrane is constant. For the second case, we discuss drug release out of a scaffold and follow the concentration profile to assess the rapidity of the drug release. To account for the importance of the relaxation time of the drug molecules and the non-local effects, we propose the formulation of mass transport using Extended Non-Equilibrium Thermodynamics ^[18,19], whose main characteristic is to upgrade the flux and higher order fluxes to the rank of independent variables.

2. Extended equation for the mass flux through a cell membrane

At micro and nanoscales, mass transport is mostly influenced by non-local effects and relaxation times. The classical Fick law

$\mathit{\boldsymbol{J}} = - \rho D\nabla c,$

(1)

relating the mass flux vector J to the mass fraction gradient ∇c, with ρ and D denoting the density and mass diffusion coefficient, respectively, is not applicable at short diffusional time (with respect to the relaxation time, i.e. large relaxation times) and small spatial scales. In order to account for large relaxation times, Fick's law can be generalized by a Cattaneo-like equation (in analogy to the one used for heat transport ^[20]) under the form

${\tau _r}{\partial _t}\mathit{\boldsymbol{J}} + \mathit{\boldsymbol{J}} = - \rho D\nabla c$

(2)

where τ_r designates the relaxation time of the mass flux and ∂_t the partial time derivative. To account furthermore for non-local effects, Extended Non-Equilibrium Thermodynamics (see for underlying developments in previous work ^[18]) proposes a further generalization

${\tau _r}{\partial _t}\mathit{\boldsymbol{J}} + \mathit{\boldsymbol{J}} = - \rho D\nabla c + {\ell ^2}{\nabla ^2}\mathit{\boldsymbol{J}}$

(3)

It is convenient to reformulate (3) in dimensionless form. Therefore, the mass flux J is rescaled with $\frac{{\rho D}}{\delta }$ , the space coordinate with a characteristic length δ (standing also for the membrane thickness), the time t with the characteristic diffusion time ${{\rm{ \mathsf{ τ} }}_d} = \frac{{{\delta ^2}}}{D}$ . Introducing the Knudsen number $Kn \equiv \frac{\ell }{\delta }$ , a dimensional number $M = \frac{{{\tau _r}}}{{{\tau _d}}}$ and using the same notation for the non-dimensional quantities, Eq (3) becomes

$M{\partial _t}\mathit{\boldsymbol{J}} + \mathit{\boldsymbol{J}} = - \nabla c + K{n^2}{\nabla ^2}\mathit{\boldsymbol{J}}$

(4)

It should be noted here that by choosing τ_d as time scale, one implies that any change in M means a change in the relaxation time, which evokes non-local effects. This motivates proposing an adaptation for the expression for M. Indeed, kinetically speaking, one may define a relation $D \equiv \frac{1}{3}\frac{{{\ell ^2}}}{{{\tau _r}}}$ , in analogy to the thermal conductivity ^[1]. Using the definition of τ_d and Kn, one finds a relation $M = \frac{1}{3}K{n^2}$ . However, in suspensions, these kinetic considerations are not always valid and differences could occur. Acknowledging that some partial proportionality is admitted between M and the Kn number (the relaxation time and the mean free path are related in that they both can cause non-local effects), this motivates to define rather $M = {f_{Kn}}\frac{1}{3}K{n^2}$ . Note that using this definition does not change anything for the analysis to be performed (indeed, in principle, any number could be chosen for M in the framework of a dimensionless analysis conditioned by correct interpretation of the results). Moreover, f_Kn is a factor that accounts for the observation that in nanofluids the relaxation time is not necessarily fully proportional to the Kn number. This factor also allows investigating separately the effects of the relaxation time and the mean free path on drug diffusion.

We can see clearly differences between (3) and the conventional constitutional Fick's law for mass diffusion (1), but also some similarities. When we deal with systems where the characteristic size is found to be much larger than the mean free path of the molecules l, Eq (3) reduces to the Cattaneo-like equation (2). This would also correspond to neglect any higher order of mass fluxes. Now, when the changes in the mass flux occur in a time span that is much larger than the relaxation time, we obtain Fick's law. Let us stress that the latter issue depends on the physicochemical relation the medication molecule can have with the carrier fluid or the host environment. Other well-known methods that are used for solving models that deal with nanoscale systems are the phonon-Boltzmann equation, Monte Carlo simulations and Molecular Dynamics. The Monte Carlo simulations and Molecular Dynamics method are numerical methods that generally do not require higher order variables, but instead demand high computing times, whilst being less flexible for modelling various scales at a time, although some efforts have been performed concerning multi-scale molecular dynamics ^[21]. Nonetheless, the flexibility in scale remains less than the present approach. The phonon-Boltzmann equation also demands high computational time for solving a model. The aforementioned aspects are generally less of a problem when it comes to the model that is proposed in this work. It should be noted that the present development could be extended by introducing non-linear mass fluxes or higher-order (higher than the second order) fluxes. However, implementing such an extension, albeit theoretically not an impossible task (see ^[18]), would introduce numerical difficulties and higher computing time, which is an absent aspect of the present approach. Furthermore, it would also loose the ability to propose analytical solutions for some particular cases, offering a pedagogical presentation allowing easier interpretations. Non-linear contributions of the mass flux (at first and higher orders) could be interesting to consider when the length scale becomes much smaller than the mean free path of the molecules. For the purposes of nanomedicine applications, this would in many cases lead to unnecessary complications. Indeed, for many cases a linear mass flux limited to second-order mass fluxes (which leads to the term ∇²J in (3)) is sufficient. In comparison, it already appeared to be sufficient for heat transport applications ^[22] and for predicting nanoflow properties ^[23].

3. Nanomedicine: non-local diffusion

3.1. Membrane permeation

From typical clinical observations that lead to data in the literature, it appears that one of the most important issues to be studied for improving drug delivery are through enhancing the permeability. In porous media, the concentration gradient can be negligible with respect to the flux, which is often the case for nanoporous flow that is pressure-driven by some external pressure. There are some cases in the human body that are comparable from a conceptual point of view, such as oxygen and carbon dioxide transfer from the lungs to the veins within the alveoles. Nonetheless, porous flow within cell membrane pores presents often the difference that, due to the absence of externally imposed pressure, the concentration gradient can be quite important with respect to the flux. In case of passive permeability, one speaks often of diffusional processes. As the diffusion process through the membrane is often the limiting pathway, the flux across the membrane will undergo non-local effects. Flux variation in the perpendicular direction to the flux across the membrane can be neglected, assuming a well-dispersed medication solution at the entrance of the whole membrane. One would then use (3) to describe such a process, reduced to a one-dimensional equation, replacing the ∇ operator by the derivative ∂_y. However, the concentration gradient is not simply a gradient of the concentrations from both sides of the membrane, since concentration jumps can exist at the entry or exit of those membranes due to lipid bilayer interaction or wettability effects. The gradient that counts is the one within the membrane. Therefore, it is important to consider the fact that flux by diffusion across a membrane (even if it were uniform) depends on the solubility. As there are different solubilities at both sides of the membrane and within the membrane as well, we can define partition coefficients ϑ_in and ϑ_out for the entry and exit of the membrane, respectively:

${\vartheta _{in}} = \frac{{{N_{mb, in}}}}{{{N_{in}}}}$

(5)

${\vartheta _{out}} = \frac{{{N_{mb, out}}}}{{{N_{out}}}}$

(6)

where N_{mb, in} and N_{mb, out} are the concentrations within the membrane from the entry and exit sides, respectively. Also, N_in and N_out stand for the same concentrations, but at the medium sides (outside the membrane), respectively. It is more convenient to write (5) and (6) in mass fractions, c, where the subscripts mean the same thing as for the concentrations. Assuming ϑ_in = ϑ_out≡ϑ, we have

$\vartheta = \frac{{{c_{mb, in}}\left( {1 - {c_{in}}} \right)}}{{{c_{in}}\left( {1 - {c_{mb, in}}} \right)}}$

(7)

$\vartheta = \frac{{{c_{mb, out}}\left( {1 - {c_{out}}} \right)}}{{{c_{out}}\left( {1 - {c_{mb, out}}} \right)}}$

(8)

where ϑ is to be determined experimentally as a material property. We can readily see that for a higher partition coefficient, the concentration gradient within the membrane is then larger than on the outer sides of the membrane, which results into a higher membrane solubility and thus a higher permeation. To model the whole process, next to the one-dimensional version of Eq (3), the mass fraction evolution is given by the classical species balance equation, which in dimensionless form gives

${\partial _t}c = - \nabla \cdot \mathit{\boldsymbol{J}}$

(9)

Since the thickness of a cell membrane is much smaller than its size in the other directions, we consider a one-dimensional model through the membrane with thickness δ, assuming that in all directions other than across the membrane the variables remain constant. The entrance of the membrane (at the drug donor side) is denoted by the dimensionless coordinate y = 0 and the exit (at the drug receptor side) by y = 1. The transport of the drug molecules through a cell membrane will then be modeled by the one-dimensional versions of Eqs. (4) and (9)

$M{\partial _t}\mathit{\boldsymbol{J}} + \mathit{\boldsymbol{J}} = - {\partial _y}c + K{n^2}\partial _y^2\mathit{\boldsymbol{J}}$

(10)

${\partial _t}c = - {\partial _y}\mathit{\boldsymbol{J}}$

(11)

Here, we omitted the subscript mb for the membrane and will do so in the following. The drug at the donor side of, but outside, the membrane is given by c_in. The inner-side drug mass fraction of the membrane at y = 0 is then, via the partition theory, given by

${\left. c \right|_{y = 0}} = \frac{{\vartheta {c_{in}}}}{{1 + \left( {\vartheta --1} \right){c_{in}}}}$

(12)

At the cell-side of the membrane, y = 1, we assume a continuity of fluxes so that

$- D\frac{{\partial c}}{{\partial y}} = - {D_{ce}}\frac{{\partial {c_{ce}}}}{{\partial y}}$

(13)

where the subscript ce denotes the inner cell. Assuming that the inner cell is an unstirred well-mixed environment, the drug uptake can be considered to occur under a constant rising concentration. Therefore, the gradient is nearly zero, so that we can say that $\frac{{\partial {c_{ce}}}}{{\partial y}} \approx 0$ . This allows to state approximately that at y = 1

${\left. {\frac{{\partial c}}{{\partial y}}} \right|_{y = 1}} = 0$

(14)

The initial mass fraction within the membrane is given by

${\left. c \right|_{t = 0}} = {c_0}$

(15)

3.2. Medication delivery through a scaffold

A supplementary interest is the continuous delivery of drugs out of a scaffold. In this case, a scaffold, containing the drugs, is introduced into the area of interest (the manner in which this is done is not the point of discussion here). This scaffold is a porous structure, where we are interested to provide a tool for analysing the diffusional flow through this medium. One side of the scaffold contains the drugs, maintained at a certain concentration or mass fraction, c_in. The other side is where the drugs enter the area of interest, where we assume that the uptake is immediate so that the mass fraction is zero there, i.e. c₀ = 0. Let us approximate such a scaffold as cylindrical pores. Contrary to the previous subsection, the dimension of the flux direction (axial direction of the pores) is larger than the radius of such a cylinder (radial direction of the pores). In that case, we can say (and also motivated by previous work in porous systems ^[23]) that the flux in the axial direction of the pores (the direction of drug release) will be approximately constant. However, the flux will be variable over the radial direction of the flux. Indeed, due to the small pore size with respect to the axial dimension of the pores, the wall effects will influence the flux and a flux profile exists in the radial direction. As for the mass fraction gradient (being a slow variable), its radial gradient is negligible with respect to the axial one (influenced by the drug release). As such, Eq (4) gives in dimensionless form (using $\frac{{\rho D}}{R}$ , R and L as scaling for the mass flux, the radial spatial coordinate and the axial spatial coordinate, respectively, where R is the pore radius and L the length of the pore)

$\mathit{\boldsymbol{J}} = - {\rm{ \mathsf{ α} }}\frac{{\partial c}}{{\partial x}} + Kn_r^2\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial r}}} \right),$

(16)

with $K{n_r} \equiv \frac{\ell }{R}$ , r the radial spatial coordinate, x the axial spatial coordinate in the pores and ${\rm{ \mathsf{ α} }} \equiv \frac{R}{L}$ . Assuming maximum flux in the middle of the pore (r = 0), i.e. ${\left. {\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial r}}} \right|_{r = 0}} = 0$ , and first-order slip (defining ${B_s} \equiv \frac{{{\ell _s}}}{R}$ , where ℓ_s stands for the slip length, related to the wall material) at the wall (r = 1), i.e.

$\mathit{\boldsymbol{J}} = - {\left. {{B_s}\frac{{\partial \mathit{\boldsymbol{J}}}}{{\partial r}}} \right|_{r = 1}},$

(17)

the flux can be found as a function of r. Taking subsequently an average value over a cross-section of the cylinder gives the following mass flux

$\mathit{\boldsymbol{J}} = - {\rm{ \mathsf{ α} }}\frac{{\partial c}}{{\partial x}}\left( {1 - \frac{{2Kn_r^2{\cal B}\left[ {1, \frac{1}{{K{n_r}}}} \right]}}{{K{n_r}{\cal B}\left[ {0, \frac{1}{{K{n_r}}}} \right] + {B_s}{\cal B}\left[ {1, \frac{1}{{K{n_r}}}} \right]}}} \right)$

(18)

with ${\cal B}$ the Bessel-I function. It should be noted that Eq (18) depends on how slip is observed and how the mass flux is defined. Other effects can influence the mass flux, but here it is the purpose to show how physical models can be used in the field of nanomedicine. Such an effective flux (taking into account the material of the scaffold (via the slip factor) and non-local effects via the Knudsen number) is now filled in the dimensionless (time scale is taken as $\frac{{{L^2}}}{D}$ ) species balance in the flow direction

${\partial _t}c = - \frac{1}{{\rm{ \mathsf{ α} }}}{\partial _x}\mathit{\boldsymbol{J}}$

(19)

The evolution of the mass fraction is then described as

${\partial _t}c = \frac{{{\partial ^2}c}}{{\partial {x^2}}}\left( {1 - \frac{{2K{n^2}{\cal B}\left[ {1, \frac{1}{{Kn}}} \right]}}{{Kn\;{\cal B}\left[ {0, \frac{1}{{Kn}}} \right] + {B_s}\;{\cal B}\left[ {1, \frac{1}{{Kn}}} \right]}}} \right)$

(20)

The boundary condition on one side of the scaffold is given by

${\left. {\frac{{\partial c}}{{\partial x}}} \right|_{x = 0}} = 0$

(21)

at x = 0, which stands for assuming that the drug only goes on one direction and

${\left. c \right|_{x = 1}} = {c_0}{e^{ - t}}$

(22)

at x = 1, which allows for a drug uptake relaxation time of the same order of magnitude as the characteristic diffusion time (which is reasonable, considering that transport in cells is mainly of diffusive character), describing the mass fraction going from c₀ to 0 (complete drug uptake). This depends of course of the system in question and the area in the body. Nevertheless, studying the effect of non-locality and slip factor on drug uptake, the choice of a relaxation time is not of importance here and does not influence our analysis. Eq (22) also shows the initial condition by evaluating it at t = 0.

3.3. Numerical method

The model is solved using a semi-implicit difference method for the spatial discretization and the forward Euler method for the temporal discretization for each time step t. At the beginning of the simulation, the initial conditions are imposed. At a certain time t, the equations for the mass flux and mass fractions are written in matrix form A*V = B, where A is a matrix of dimension (n-2)×(n-2) (n being the number of nodes in the discretized spatial domain and two nodes being dedicated to the boundary conditions, the bulk equations are thus discretized on (n-2) nodes) containing the coefficients in the equations for the mass fractions and mass fluxes at time t+1, V a vector of dimension (n-2) containing the unknown mass fraction and mass flux values at the spatial nodes at time t+1 and B a vector of dimension (n-2) containing the known mass fraction and mass flulx values at time t. Due to the spatial gradients in the bulk equations, the values at nodes 2 and n-1 depend on the boundary values at nodes 1 and n, respectively. The vector V is calculated as V = A^-1B. The boundary condition for at the membrane entrance is imposed and held at that value at each time step. The boundary condition at the exit of the membrane at time t+1 is then obtained through the bulk values at time t+1. This procedure is repeated until the relative errors of the bulk and boundary values of both the mass fraction and mass flux between the previous and present loop become smaller than 10^-4. The obtained values at time t+1 are stocked in the output matrix and used as the known values for the next time step t+2.

4. Results and discussion

4.1. Membrane diffusion

We use our model, given by Eqs. (10)–(12), (14)–(15), and solve it by means of a semi-implicite finite difference algorithm, explained in section 4.1. We take c_in = 0.8, c₀ = 0 and ϑ = 1 (we neglect the effect of the partition coefficient for the moment). The results are presented in Figures 1-5. Figures 1-3 show the evolution of the mass fraction across the membrane for Kn numbers Kn = 0.1, Kn = 1 and Kn = 10, respectively. The different Kn numbers stand for different sizes of the membrane with respect to the molecular mean free path of the drug molecules in the fluid present in the membrane (which can be taken as a watery solution). Furthermore, for each Kn number different f_Kn values are used: f_Kn = 0.1, f_Kn = 1 and f_Kn = 10.

Figure 1. Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 0.1.

DownLoad: Full-Size Img PowerPoint

Figure 2. Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 1.

DownLoad: Full-Size Img PowerPoint

Figure 3. Mass fraction evolution across the membrane layer for different values of f_Kn for Kn = 10.

DownLoad: Full-Size Img PowerPoint

Figure 4. Mass fraction as a function of time at y = 1 for different values of f_Kn and Kn = 0.1 (left), 1 (middle) and 10 (right).

DownLoad: Full-Size Img PowerPoint

Figure 5. A zoom around y = 1 for several f_Kn values for Kn = 1 (left) and for several Kn values for f_Kn = 1 (right) at y = 0.3.

DownLoad: Full-Size Img PowerPoint

We can see that for a given Kn number a higher f_Kn number (higher relaxation time with respect to a fixed diffusion characteristic time), the evolution of the mass flux is significantly slower, with a large gradient at the membrane entrance. This implies that the mass flux is very high at the membrane entrance but very low at the membrane exit. This effect is precisely referred to as non-locality. It is important here that by choosing the diffusion characteristic time as the scale time, any change in f_Kn implies a change in the relaxation time. This means that one compares systems with different relaxation times, which illustrates the effect of medication diffusion as a function of different medication-carrier fluid systems. A higher Kn number has the same non-local effect, albeit on the geometrical level (for the same system). This latter effect is not so clearly visible in Figures 1–3, but can be more appreciated looking to the evolution of the mass fraction at the membrane exit y = 1 in Figure 4.

Effectively, Figure 4 shows that a higher Kn number, for all f_Kn values, results into a lower mass fraction, which is representative for a lower drug diffusion/solubility in the membrane. The non-local effect at high Kn numbers is responsible for a non-Fickian mass fraction profile, which is indicated by a mass fraction at the exit of the membrane that is less than the initial mass fraction at steady state (this in comparison with Fickian behaviour, where for a zero-flux boundary condition, the profile should become flat). This means that for drugs that have large relaxation times or mean free paths, the drug uptake could be troublesome. It is interesting to note that the values shown in Figure 4 are exactly at the membrane exit, but that a mass fraction jump exists just before that exit. Not being well discernible in Figures 1-3, a zoom of some of those curves has been plotted in Figure 5.

This mass fraction jump can be seen clearly in Figure 5 and stems from the non-local behaviour from both the mass flux time derivative (∂_tJ) and mass flux diffusion term (∇²J). Now, if the partition coefficient was not at unity, we would have not only a jump due to a non-locality at the membrane side of the interface, but also at the cell side, due the partition coefficient. Figure 6 shows the influence of the partition coefficient on the mass fraction evolution.

Figure 6. Mass fraction across the membrane at t = 3 for f_Kn = Kn = 1 (left) and as function of time at y = 1 for f_Kn = Kn = 1 (right).

DownLoad: Full-Size Img PowerPoint

It can be seen that a higher partition coefficient results into a higher mass fraction and a higher drug uptake (higher mass fraction at the membrane exit). Since a zero-flux boundary condition (standing for a quick drug diffusion in the cell and not for a wall-like boundary, as is discussed under Eq (13)) is maintained at the membrane exit, the partition coefficient does not affect the mass fraction here. However, it does effect the quantitative value of the mass fraction just outside the membrane (in the cell). The higher the partition coefficient, the lower the mass fraction will be here. Also, as Figure 6 seems to show, a higher partition coefficient does not exhibit a different behaviour of the mass fraction at the membrane exit. That is, a higher partition coefficient gives a higher mass fraction of the drug at the membrane entrance, which translates itself into a higher mass fraction at the membrane exit. Nonetheless, this increase is cancelled by the partition coefficient itself at the membrane exit. This motivates using drugs that have a higher partition coefficient at the membrane entrance, but a lower one at the membrane exit. This would, of course, not lead to a higher total uptake (mass conservation implies that the total quantity of the uptake is equal to the administered quantity), but may lead to a higher efficiency if a higher mass fraction is needed for a better healing in a shorter time period. Having discussed drug uptake by means of drug diffusion through a cell membrane, it is also interesting to provide tools to drug release by scaffolds, applicable in the field of nanomedicine. This is treated in the next subsection.

4.2. Scaffold drug release

We use our model, given by Eqs (20)–(22), and solve it in the same way as for the non-local drug diffusion model. We take c₀ = 0.8 and present the results in Figures 7 and 8 for different Kn numbers and slip factors B_s. When we take B_s = 0, this stands for a zero slip length of the mass flux. This means that the mass flux at the wall of the pores in the scaffold is equal to zero. However, when B_s→∞, this means that there is full slip at the pore walls and that the mass flux is everywhere equal in the radial direction. We take a value B_s = 10 to be quite large, standing for "large slip" and compare its results to "no-slip". The slip behaviour of the mass flux is not clearly defined in the literature, but is rather often used in the framework of fluid flow, related to its velocity. The very definition of a mass flux in a porous medium being ^[23]

Figure 7. Mass fraction across the scaffold pores for several times t for Kn = 0.1 (left), 1 (middle) and 10 (right) with a no-slip condition at the pore walls, B_s = 0.

DownLoad: Full-Size Img PowerPoint

Figure 8. Mass fraction across the scaffold pores for several times t for Kn = 0.1 (left), 1 (middle) and 10 (right) with a high-slip condition at the pore walls, B_s = 10.

DownLoad: Full-Size Img PowerPoint

$\mathit{\boldsymbol{J}} = \phi {\rho _f}\left( {{\mathit{\boldsymbol{v}}_f} - \mathit{\boldsymbol{v}}} \right)$

(23)

where ϕ is the porosity of the porous medium, ρ_f the density of the fluid flowing through it, v_f the fluid velocity and v the barycentric velocity, we see that the mass flux is directly related to the fluid velocity. If we note that in a two-component system the barycentric velocity is expressed as ρv = ϕρ_fv_f+(1-ϕ)ρ_sv_s, wherein the total mass density ρ is ρ = ϕρ_f+(1-ϕ)ρ_s with ρ_s and v_s the density and velocity of the solid component, respectively, and that in case of a scaffold v_s = 0, we can quickly see that in (23) the only variable the mass flux is related directly to is the fluid velocity v_f. As such, the discussion on the slip factor becomes more physical and understandable.

Figures 7 and 8 show that the Kn number tends to delay the drug release. This suggests that for pores of sizes of the order of magnitude of the mean free path of the drug molecules the drug release is slowed down. On the other hand, in some applications, it is necessary to have pores of such a size. Figure 8 shows that the slower drug release can be compensated by a higher slip factor, which accelerates the release of medication, illustrated by mass fraction profiles that evolve faster with respect to Figure 7. Figure 9 shows the mass fraction evolution for different Kn numbers for "no-slip" (B_s = 0) and "high-slip" (B_s = 10) at the membrane exit (y = 1). It can be seen that indeed, when B_s is increased, the reduction of the mass fraction due to an increase in the Kn number (from Kn = 0.1 to Kn = 1) is compensated, resulting into a similar profile. This confirms that by increasing the slip factor, which amounts to say altering the material properties of the scaffold pores to properties that reduce the interaction with the drug molecule solutions, compensates the reduction of drug release in nanopores.

Figure 9. Mass flux as a function of time at x = 1 for different values of Kn and B_s = 0 (left) and 10 (right).

DownLoad: Full-Size Img PowerPoint

5. Conclusions

This work presents new developments in Extended Non-Equilibrium Thermodynamics, proposing a model for studying drug uptake through cell membranes and drug release out of scaffolds. The theory is based on elevating fluxes to the same independent status of the usual variables. In this particular case, the mass flux has the same status as the mass fraction or concentration of the medication solution. This means that the mass flux follows its own evolution next to that of the mass fraction, being connected on an equal level via a species balance. By taking into account the relative importance of the relaxation time with respect to the diffusion characteristic time, and of the mean free path with respect to the system characteristic size, non-local effects are included in this way.

As for the diffusion through a cell membrane, it appeared that for sizes and relaxation times of the order of magnitude of or smaller than the mean free path (Kn = O(1)) and diffusion characteristic time (M = O(1)), respectively, leads to a lower drug uptake. This is due to the non-local effect, mentioned before. Non-Fickian profiles result from this, where it can be seen that at high non-locality a large mass fraction gradient is present close to the membrane entrance, but a nearly horizontal one near the exit. The latter contributes to the reduction of drug uptake. Furthermore, it has also been shown that at the membrane exit, a mass fraction jump (within the membrane) occurs that is not due to the boundary condition, but also due the non-locality effect. This is an effect that has not only been observed for heat transport ^[25], where it is already an established phenomenon, but also for mass transport ^[26]. The effect of the partition coefficient on the mass fraction is also studied, since they are typically introduced to account for the typical mass fraction jump (outside the membrane) observed in membrane diffusion studies, though still presenting experimental uncertainty to its precise value ^[26,27,28]. No particular effect has been observed other than that a higher partition coefficient causes a higher mass fraction jump at the membrane boundaries. If the medication could be chosen in such a way that the partition coefficient is high at the entrance and low at the exit, the efficiency of the medication could be increased if a higher concentration would lead to a better healing. Of course, mass conservation implies that this would not influence the total quantity of the uptake, but rather the speed of uptake, which could influence the kinetics of healing.

As for the drug release from the scaffold, the non-locality effect also showed to have a reducing impact on the release of drugs from the scaffold. It appeared that such a reduction could be countered by choosing the material of the scaffold in such a way that (besides biocompatibility, which is not a complicated task on its own) the interaction between the scaffold pore walls and the medication solution is weakened. This results into a higher slip factor enhancing drug release.

Future developments of this work include experimental demonstration, using the proposed model in this work. As quantitative agreement between experiments and models are difficult to find in the literature, comparison with experimental data should first be performed by nanomaterials of which material properties are well known and documented. As example, we can mention the delivery of nucleic acids with DNA nanostructures as new means of drug delivery by nanocarriers ^[29,30]. In the case of the nanocarriers, which can also be seen as scaffolds, the diffusion of the carriers through the cell membranes could be described by the set of Eqs (10)–(12) and (14)–(15). Subsequently, the delivery of the drugs out of the carriers into the cells could be described as release from a scaffold by Eqs (20)–(22). The difficulty lies in determining the interaction the medication can have with the nanocarrier, such as the DNA nanostructures ^[29]. In other words, it is important to know the slip length or slipping coefficient in Eq (17). Enhancement of drug release does not depend only on the release of the drug out of the carrier, but also the solubility of the carrier itself through the cell membranes. It is reported that DNA nanostructures can overcome the usually large circulation times of carriers ^[29,30]. The aforementioned discussion implies a necessity of measuring some material properties in order to use the presented model in this paper. In some cases, information about material properties could be obtained by in vitro experiments with a controlled experimental setup in a similar way performed by ^[31] for high-precision regulation of cell behaviour. This brings the intermediate step between modelling and full experimental demonstration to be identified by bio-mimicking cell membranes or scaffolds. Finally, such information can complete an experimental demonstration of our model.

Acknowledgements

Financial support from BelSPo is acknowledged.

Conflict of interest

The author declares there is no conflict of interest.

References

[1]	Abadie MO, Mendes N (2008) Numerical assessment of turbulence effect on the evaluation of wind-driven rain specific catch ratio. Int Commun Heat Mass Transfer 35: 1253–1261. doi: 10.1016/j.icheatmasstransfer.2008.08.013
[2]	Al-Hajraf S, Rubini P (2001) Three-dimensional homogeneous two-phase flow modelling of drifting sand around an open gate. WIT Trans Eng Sci 30: 309–325.
[3]	Alhajraf S (2000) Numerical simulation of drifting sand, PhD thesis, Cranfield University.
[4]	Anderson RS, Haff PK (1988) Simulation of eolian saltation. Science 241: 820–823. doi: 10.1126/science.241.4867.820
[5]	Anderson RS, Hallet B (1986) Sediment transport by wind: Toward a general model. Geol Soc Am Bull 97: 523–535. doi: 10.1130/0016-7606(1986)97<523:STBWTA>2.0.CO;2
[6]	Andreotti B, Claudin P, Douady S (2002a) Selection of dune shapes and velocities part 1: Dynamics of sand, wind and barchans. Eur Phys J B 28: 321–339.
[7]	Andreotti B, Claudin P, Douady S (2002b) Selection of dune shapes and velocities part 2: A two-dimensional modelling. Eur Phys J B 28: 341–352.
[8]	Andrews M, O'rourke P (1996) The multiphase particle-in-cell (MP-PIC) method for dense particulate flows. Int J Multiphase Flow 22: 379–402. doi: 10.1016/0301-9322(95)00072-0
[9]	Arastoopour H, Pakdel P, Adewumi M (1990) Hydrodynamic analysis of dilute gas-solids flow in a vertical pipe. Powder Technol 62: 163–170. doi: 10.1016/0032-5910(90)80080-I
[10]	Balachandar S, Eaton J (2010) Turbulent dispersed multiphase flow. Annu Rev Fluid Mech 42: 111–133. doi: 10.1146/annurev.fluid.010908.165243
[11]	Bang B, Nielsen A, Sundsbø P, et al. (1994) Computer simulation of wind speed, wind pressure and snow accumulation around buildings (SNOW-SIM). Energy Build 21: 235–243. doi: 10.1016/0378-7788(94)90039-6
[12]	Barnea E, Mizrahi J (1973) A generalized approach to the fluid dynamics of particulate systems: Part 1. General correlation for fluidization and sedimentation in solid multiparticle systems. Chem Eng J 5: 171–189.
[13]	Benyahia S, Syamlal M, O'Brien TJ (2006) Extension of Hill-Koch-Ladd drag correlation over all ranges of Reynolds number and solids volume fraction. Powder Technol 162: 166–174. doi: 10.1016/j.powtec.2005.12.014
[14]	Beyers J, Sundsbø P (2004) Numerical simulation of three-dimensional, transient snow drifting around a cube. J Wind Eng Ind Aerod 92: 725–747. doi: 10.1016/j.jweia.2004.03.011
[15]	Beyers M, Waechter B (2008) Modeling transient snowdrift development around complex three-dimensional structures. J Wind Eng Ind Aerod 96: 1603–1615. doi: 10.1016/j.jweia.2008.02.032
[16]	Blocken B (2014) 50 years of computational wind engineering: past, present and future. J Wind Eng Ind Aerod 129: 69–102. doi: 10.1016/j.jweia.2014.03.008
[17]	Blocken B, Carmeliet J (2010) Overview of three state-of-the-art wind-driven rain assessment models and comparison based on model theory. Build Environ 45: 691–703. doi: 10.1016/j.buildenv.2009.08.007
[18]	Blocken B, Stathopoulos T, Carmeliet J (2007) CFD simulation of the atmospheric boundary layer: Wall function problems. Atmos Environ 41: 238–252. doi: 10.1016/j.atmosenv.2006.08.019
[19]	Boutanios Z, Jasak H (2017) Two-way coupled Eulerian-Eulerian simulations of drifting snow with viscous treatment of the snow phase. J Wind Eng Ind Aerod 169: 67–76. doi: 10.1016/j.jweia.2017.07.006
[20]	Businger J, Wyngaard J, Izumi Y, et al. (1971) Flux profile relationships in the atmospheric surface layer. J Atmos Sci 28: 181–189. doi: 10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2
[21]	Canuto C, Lo Giudice A (2018) A multi-timestep robin-robin domain decomposition method for time dependent advection-diffusion problems. App Math Comput (In press).
[22]	Chapman S, Cowling T (1970) The Mathematical Theory of Non-uniform Gases, Cambridge Mathematical Library.
[23]	Chen C, Wood P (1985) A urbulence closure model for dilute gas particle flows. Can J Chem Eng 63: 349–360. Available from: https://doi.org/10.1002/cjce.5450630301.
[24]	Choi E (1992) Simulation of wind-driven-rain around a building. J Wind Eng 52: 60–65.
[25]	Choi E (1997) Numerical modelling of gust effect on wind-driven rain. J Wind Eng Ind Aerod 72: 107–116. doi: 10.1016/S0167-6105(97)00246-8
[26]	Creyssels M, Dupont P, Moctar AOE, et al. (2009) Saltating particles in a turbulent boundary layer: Experiment and theory. J Fluid Mech 625: 47–74. doi: 10.1017/S0022112008005491
[27]	Deen N, Annaland MVS, der Hoef MV, et al. (2007) Review of discrete particle modeling of fluidized beds. Chem Eng Sci 62: 28–44. doi: 10.1016/j.ces.2006.08.014
[28]	Di Felice R (1994) The voidage function for fluid-particle interaction systems. Int J Multiphase Flow 20: 153–159. doi: 10.1016/0301-9322(94)90011-6
[29]	Durán O, Parteli EJ, Herrmann HJ (2010) A continuous model for sand dunes: Review, new developments and application to barchan dunes and barchan dune fields. Earth Surf Processes Landforms 35: 1591–1600. doi: 10.1002/esp.2070
[30]	Dyer J (1974) A review of flux profile relationships. Boundary-Lay Meteorol 7: 363–372. doi: 10.1007/BF00240838
[31]	Elghobashi S (1991) Particle-laden turbulent flows: Direct simulation and closure models, In: Oliemans, R.V.A. Editor, Computational Fluid Dynamics for the Petrochemical Process Industry, Dordrecht: Springer, 91–104.
[32]	Elghobashi S (1994) On predicting particle-laden turbulent flows. Appl Sci Res 52: 309–329. doi: 10.1007/BF00936835
[33]	Ergun S (1952) Fluid flow through packed columns. J Chem Eng Prog 48: 89–94.
[34]	Faber TE (1995) Fluid Dynamics for Physicists, Cambridge University Press.
[35]	Farimani AB, Ferreira AD, Sousa AC (2011) Computational modeling of the wind erosion on a sinusoidal pile using a moving boundary method. Geomorphology 130: 299–311. doi: 10.1016/j.geomorph.2011.04.012
[36]	Gauer P (1999) Blowing and drifting snow in alpine terrain: a physically-based numerical model and related field measurements. PhD thesis, ETH Zurich, Switzerland.
[37]	Germano M, Piomelli U, Moin P, et al. (1991) A dynamic subgrid-scale eddy viscosity model. Phys Fluids A 3: 1760–1765. doi: 10.1063/1.857955
[38]	Gidaspow D (1994) Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions, Academic press.
[39]	Gosman A, Lekakou C, Politis S, et al. (1992) Multidimensional modeling of turbulent two-phase flows in stirred vessels. AIChE J 38: 1946–1956. doi: 10.1002/aic.690381210
[40]	Hangan H (1999) Wind-driven rain studies. A C-FD-E approach. J Wind Eng Ind Aerod 81: 323–331. doi: 10.1016/S0167-6105(99)00027-6
[41]	Ho TD, Dupont P, Ould El Moctar A, et al. (2012) Particle velocity distribution in saltation transport. Phys Rev E 85: 052301. Available from: https://doi.org/10.1103/PhysRevE. 85.052301.
[42]	Ho TD, Valance A, Dupont P, et al. (2011) Scaling laws in aeolian sand transport. Phys Rev Lett 106: 094501. Available from: https://doi.org/10.1103/PhysRevLett.106.094501. doi: 10.1103/PhysRevLett.106.094501
[43]	Ho TD, Valance A, Dupont P, et al. (2014) Aeolian sand transport: Length and height distributions of saltation trajectories. Aeolian Res 12: 65–74. doi: 10.1016/j.aeolia.2013.11.004
[44]	Hrenya CM, Sinclair JL (1997) Effects of particle-phase turbulence in gas-solid flows. AlChE J 43: 853–869. Available from: https://aiche.onlinelibrary.wiley.com/doi/abs/10. 1002/aic.690430402. doi: 10.1002/aic.690430402
[45]	Hsu TJ, Jenkins JT, Liu PL (2004) On two-phase sediment transport: Sheet flow of massive particles. Proc R Soc A 460: 2223–2250. Available from: https://doi.org/10.1098/rspa. 2003.1273. doi: 10.1098/rspa.2003.1273
[46]	Huang S, Li Q (2010a) A new dynamic one-equation subgrid-scale model for large eddy simulations. Int J Numer Methods Eng 81: 835–865.
[47]	Huang S, Li Q (2010b) Numerical simulations of wind-driven rain on building envelopes based on Eulerian multiphase model. J Wind Eng Ind Aerod 98: 843–857.
[48]	Huang S, Li Q (2011) Large eddy simulations of wind-driven rain on tall building facades. J Struct Eng 138: 967–983.
[49]	Iversen J, Greeley R, White BR, et al. (1975) Eolian erosion of the martian surface, part 1: Erosion rate similitude. Icarus 26: 321–331. doi: 10.1016/0019-1035(75)90175-X
[50]	Iversen J, Rasmussen K (1999) The effect of wind speed and bed slope on sand transport. Sedimentology 46: 723–731. Available from: https://doi.org/10.1046/j.1365-3091. 1999.00245.x. doi: 10.1046/j.1365-3091.1999.00245.x
[51]	Jenkins JT, Hanes DM (1998) Collisional sheet flows of sediment driven by a turbulent fluid. J Fluid Mech 370: 29–52. Available from: https://doi.org/10.1017/S0022112098001840. doi: 10.1017/S0022112098001840
[52]	Ji S, Gerber A, Sousa A (2004) A convection-diffusion CFD model for aeolian particle transport. Int J Numer Methods Fluids 45: 797–817. Available from: https://doi.org/10.1002/fld. 724. doi: 10.1002/fld.724
[53]	Jiang H, Dun H, Tong D, et al. (2017) Sand transportation and reverse patterns over leeward face of sand dune. Geomorphology 283: 41–47. doi: 10.1016/j.geomorph.2016.12.030
[54]	Jiang H, Huang N, Zhu Y (2014) Analysis of wind-blown sand movement over transverse dunes. Sci Rep 4: 7114.
[55]	Jones W, Launder B (1972) The prediction of laminarization with a two-equation model of turbulence. Int J Heat Mass Transfer 15: 301–314. doi: 10.1016/0017-9310(72)90076-2
[56]	Kang L (2012) Discrete particle model of aeolian sand transport: Comparison of 2D and 2.5D simulations. Geomorphology 139: 536–544.
[57]	Kang L, Guo L (2006) Eulerian-Lagrangian simulation of aeolian sand transport. Powder Technol 162: 111–120. doi: 10.1016/j.powtec.2005.12.002
[58]	Kang L, Liu D (2010) Numerical investigation of particle velocity distributions in aeolian sand transport. Geomorphology 115: 156–171. doi: 10.1016/j.geomorph.2009.10.001
[59]	Kang L, Zou X (2011) Vertical distribution of wind-sand interaction forces in aeolian sand transport. Geomorphology 125: 361–373. doi: 10.1016/j.geomorph.2010.09.025
[60]	Kato M, Launder B (1993) The modeling of turbulent flow around stationary and vibratingsquare cylinders, In: Ninth Symposium on Turbulent Shear Flows, American Society of Mechanical Engineers, 1–6.
[61]	Kobayashi H (2005) The subgrid-scale models based on coherent structures for rotating homogeneous turbulence and turbulent channel flow. Phys Fluids 17: 045104. doi: 10.1063/1.1874212
[62]	Kobayashi H, Ham F, Wu X (2008) Application of a local sgs model based on coherent structures to complex geometries. Int J Heat Fluid Flow 29: 640–653. doi: 10.1016/j.ijheatfluidflow.2008.02.008
[63]	Kok JF, Parteli EJ, Michaels TI, et al. (2012) The physics of wind-blown sand and dust. Rep Prog Phys 75: 106901. doi: 10.1088/0034-4885/75/10/106901
[64]	Kolmogorov AN (1941) The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl Akad Nauk SSSR 30: 299–303.
[65]	Kubilay A, Carmeliet J, Derome D (2017) Computational fluid dynamics simulations of wind-driven rain on a mid-rise residential building with various types of facade details. J Build Perform Simul 10: 125–143. doi: 10.1080/19401493.2016.1152304
[66]	Kubilay A, Derome D, Blocken B, et al. (2013) CFD simulation and validation of wind-driven rain on a building facade with an Eulerian multiphase model. Build Environ 61: 69–81. doi: 10.1016/j.buildenv.2012.12.005
[67]	Kubilay A, Derome D, Blocken B, et al. (2014) Numerical simulations of wind-driven rain on an array of low-rise cubic buildings and validation by field measurements. Build Environ 81: 283–295. doi: 10.1016/j.buildenv.2014.07.008
[68]	Kubilay A, Derome D, Blocken B, et al. (2015a) Numerical modeling of turbulent dispersion for wind-driven rain on building facades. Environ Fluid Mech 15: 109–133.
[69]	Kubilay A, Derome D, Blocken B, et al. (2015b) Wind-driven rain on two parallel wide buildings: field measurements and CFD simulations. J Wind Eng Ind Aerod 146: 11–28.
[70]	Lakehal D, Mestayer P, Edson J, et al. (1995) Eulero-Lagrangian simulation of raindrop trajectories and impacts within the urban canopy. Atmos Environ 29: 3501–3517. doi: 10.1016/1352-2310(95)00202-A
[71]	Li Z, Wang Y, Zhang Y, et al. (2014) A numerical study of particle motion and two-phase interaction in aeolian sand transport using a coupled large eddy simulation-discrete element method. Sedimentology 61: 319–332. Available from: https://doi.org/10.1111/sed. 12057. doi: 10.1111/sed.12057
[72]	Liston G, Brown R, Dent J (1993) A two-dimensional computational model of turbulent atmospheric surface flows with drifting snow. Ann Glaciol 18: 281–286. doi: 10.3189/S0260305500011654
[73]	Lopes A, Oliveira L, Ferreira AD, et al. (2013) Numerical simulation of sand dune erosion. Environ Fluid Mech 13: 145–168. Available from: https://doi.org/10.1007/s10652-012-9263-2. doi: 10.1007/s10652-012-9263-2
[74]	Lugo J, Rojas-Solorzano L, Curtis J (2012) Numerical simulation of aeolian saltation within the sediment transport layer using granular kinetic theory. Rev Fac Ing 27: 80–96.
[75]	Lun CKK, Savage SB, Jeffrey DJ, et al. (1984) Kinetic theories for granular flow: Inelastic particles in Couette flow and slightly inelastic particles in a general flowfield. J Fluid Mech 140: 223–256. Available from: https://doi.org/10.1017/S0022112084000586. doi: 10.1017/S0022112084000586
[76]	Marval JP, Rojas-Solórzano LR, Curtis JS (2007) Two-dimensional numerical simulation of saltating particles using granular kinetic theory, In: ASME/JSME 2007 5th Joint Fluids Engineering Conference, 929–939.
[77]	Menter FR (1994) Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 32: 1598–1605. doi: 10.2514/3.12149
[78]	Menter FR, Kuntz M, Langtry R (2003) Ten years of industrial experience with the SST turbulence model. Turbul heat mass transfer 4: 625–632.
[79]	Mochida A, Lun I (2008) Prediction of wind environment and thermal comfort at pedestrian level in urban area. J Wind Eng Ind Aerod 96: 1498–1527. doi: 10.1016/j.jweia.2008.02.033
[80]	Monin A, Obukhov A (1954) Basic laws of turbulent mixing in the surface layer of the atmosphere. Contrib Geophys Inst Acad Sci 151: 163–187.
[81]	Moore I (1995) Numerical modelling of blowing snow around buildings. PhD thesis, University of Leeds.
[82]	Naaim M, Florence N, Martinez H (1998) Numerical simulation of drifting snow: Erosion and deposition models. Ann Glaciol 26: 191–196. doi: 10.3189/1998AoG26-1-191-196
[83]	Nalpanis P,Hunt J,Barrett C (1993) Saltating particles over flat beds. J Fluid Mech 251: 661–685. doi: 10.1017/S0022112093003568
[84]	Okaze T, Niiya H, Nishimura K (2018) Development of a large-eddy simulation coupled with Lagrangian snow transport model. J Wind Eng Ind Aerod 183: 35–43. doi: 10.1016/j.jweia.2018.09.027
[85]	Okaze T, Takano Y, Mochida A, et al. (2015) Development of a new κ–? model to reproduce the aerodynamic effects of snow particles on a flow field. J Wind Eng Ind Aerod 144: 118–124.
[86]	Parteli E, Schwämmle V, Herrmann H, et al. (2006) Profile measurement and simulation of a transverse dune field in the lençóis maranhenses. Geomorphology 81: 29–42. doi: 10.1016/j.geomorph.2006.02.015
[87]	Pasini JM, Jenkins JT (2005) Aeolian transport with collisional suspension. Philos Trans R Soc A 363: 1625–1646. Available from: https://doi.org/10.1098/rsta.2005.1598. doi: 10.1098/rsta.2005.1598
[88]	Patankar N, Joseph D (2001) Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach. Int J Multiphase Flow 27: 1659–1684. doi: 10.1016/S0301-9322(01)00021-0
[89]	Pettersson K, Krajnovic S, Kalagasidis A, et al. (2016) Simulating wind-driven rain on building facades using eulerian multiphase with rain phase turbulence model. Build Environ 106: 1–9. doi: 10.1016/j.buildenv.2016.06.012
[90]	Pischiutta M, Formaggia L, Nobile F (2011) Mathematical modelling for the evolution of aeolian dunes formed by a mixture of sands: Entrainment-deposition formulation. Commun Appl Ind Math 2: 0003777.
[91]	Pomeroy J, Gray D (1990) Saltation of snow. Water Resour Res 26: 1583–1594. doi: 10.1029/WR026i007p01583
[92]	Preziosi L, Fransos D, Bruno L (2015) A multiphase first order model for non-equilibrium sand erosion, transport and sedimentation. Appl Math Lett 45: 69–75. doi: 10.1016/j.aml.2015.01.011
[93]	Sagaut P (2006) Large eddy simulation for incompressible flows: An introduction. Springer Science & Business Media.
[94]	Sato T, Uematsu T, Nakata T, et al. (1993) Three dimensional numerical simulation of snowdrift, In: Murakami, S. Editor, Computational Wind Engineering 1, Elsevier, Oxford, 741–746. Available from: https://doi.org/10.1016/B978-0-444-81688-7.50082-6.
[95]	Sauermann G, Andrade J, Maia L, et al. (2003) Wind velocity and sand transport on a barchan dune. Geomorphology 54: 245–255. doi: 10.1016/S0169-555X(02)00359-8
[96]	Sauermann G, Kroy K, Herrmann HJ (2001) Continuum saltation model for sand dunes. Phys Rev E 64: 031305. doi: 10.1103/PhysRevE.64.031305
[97]	Schlichting H, Gersten K (2016) Boundary-Layer Theory, Springer.
[98]	Shi X, Xi P, Wu J (2015) A lattice Boltzmann-Saltation model and its simulation of aeolian saltation at porous fences. Theor Comput Fluid Dyn 29: 1–20.
[99]	Shih TH, Liou WW, Shabbir A, et al. (1995) A new k-? eddy viscosity model for high Reynolds number turbulent flows. Comput Fluids 24: 227–238.
[100]	Smagorinsky J (1963) General circulation experiments with the primitive equations: I. the basic experiment. Mon Weather Rev 91: 99–164. doi: 10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
[101]	Sun Q, Wang G, Xu Y (2001) DEM applications to aeolian sediment transport and impact process in saltation. Partl Sci Technol 19: 339–353. doi: 10.1080/02726350290057877
[102]	Sun X, He R, Wu Y (2018) Numerical simulation of snowdrift on a membrane roof and the mechanical performance under snow loads. Cold Reg Sci Technol 150: 15–24. doi: 10.1016/j.coldregions.2017.09.007
[103]	Sundsbø P (1998) Numerical simulations of wind deflection fins to control snow accumulation in building steps. J Wind Eng Ind Aerod 74: 543–552.
[104]	Surry D, Inculet D, Skerlj P, et al. (1994) Wind, rain and the building envelope: A status report of ongoing research at the university of western ontario. J Wind Eng Ind Aerod 53: 19–36. doi: 10.1016/0167-6105(94)90016-7
[105]	Syamlal M, O' Brien T (1987) The derivation of a drag coefficient formula from velocity-voidage correlations. Tech Note, US Department of Energy, Office of Fossil Energy, NETL, Morgantown, WV.
[106]	Thiis TK (2000) A comparison of numerical simulations and full-scale measurements of snowdrifts around buildings. Wind Struct 3: 73–81. doi: 10.12989/was.2000.3.2.073
[107]	Tominaga Y, Okaze T, Mochida A (2011) CFD modeling of snowdrift around a building: An overview of models and evaluation of a new approach. Build Environ 46: 899–910. doi: 10.1016/j.buildenv.2010.10.020
[108]	Tong D, Huang N (2012) Numerical simulation of saltating particles in atmospheric boundary layer over flat bed and sand ripples. J Geophys Res Atmos 117. Available from: https://doi. org/10.1029/2011JD017424.
[109]	Uematsu T, Nakata T, Takeuchi K, et al. (1991) Three-dimensional numerical simulation of snowdrift. Cold Reg Sci Technol 20: 65–73. doi: 10.1016/0165-232X(91)90057-N
[110]	Vinkovic I, Aguirre C, Ayrault M, et al. (2006) Large-eddy simulation of the dispersion of solid particles in a turbulent boundary layer. Boundary-layer Meteorol 121: 283–311. doi: 10.1007/s10546-006-9072-6
[111]	Wang H, Hou X, Deng Y (2015) Numerical simulations of wind-driven rain on building facades under various oblique winds based on Eulerian multiphase model. J Wind Eng Ind Aerod 142: 82–92. doi: 10.1016/j.jweia.2015.02.006
[112]	Wen C, Yu Y (1966) A generalized method for predicting the minimum fluidization velocity. AIChE J 12: 610–612. doi: 10.1002/aic.690120343
[113]	Werner B (1990) A Steady-State model of Wind-Blown sand transport. J Geol 98: 1–17. Available from: https://doi.org/10.1086/629371. doi: 10.1086/629371
[114]	Wilcox D (2008) Formulation of the k-ω turbulence model revisited. AIAA J 46: 2823–2838. Available from: https://doi.org/10.2514/1.36541. doi: 10.2514/1.36541
[115]	Wilcox DC (1988) Reassessment of the scale-determining equation for advanced turbulence models. AIAA J 26: 1299–1310. doi: 10.2514/3.10041
[116]	Wilcox DC (1998) Turbulence Modeling for CFD, 2nd Edition, DCW industries La Canada, California.
[117]	Xiao F, Guo L, Li D, et al. (2012) Discrete particle simulation of mixed sand transport. Particuology 10: 221–228. doi: 10.1016/j.partic.2011.10.004
[118]	Yakhot V, Orszag S (1986) Renormalization group analysis of turbulence. i. basic theory. J Sci Comput 1: 3–51. Available from: https://doi.org/10.1007/BF01061452.
[119]	Zhou X,Kang L,Gu M,et al. (2016a) Numerical simulation and wind tunnel test for redistribution of snow on a flat roof. J Wind Eng Ind Aerod 153: 92–105.
[120]	Zhou X, Zhang Y, Wang Y, et al. (2016b) 3D numerical simulation of the evolutionary process of aeolian downsized crescent-shaped dunes. Aeolian Res 21: 45–52.

This article has been cited by:

Lu Xu, Chunlai Mu, Qiao Xin, Global boundedness and asymptotic behavior of solutions for a quasilinear chemotaxis model of multiple sclerosis with nonlinear signal secretion, 2023, 28, 1531-3492, 1215, 10.3934/dcdsb.2022118

Reader Comments

Your name:*

Email:*
© 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)