Semi-resolved CFD-DEM study on the influence of ‘surface energy’ factors on the deposition of lignin particles in ceramic membrane pores

Wang, H., Wu, J., Fu, P., Song, X., Li, L., Yang, Y., Xu, W., Li, Z., and Yang, Z. (2026). "Semi-resolved CFD-DEM study on the influence of ‘surface energy’ factors on the deposition of lignin particles in ceramic membrane pores," BioResources 21(3), 6537–6568.

Abstract

This study adopted a semi‑analytical CFD‑DEM coupling method to simulate the movement and deposition of lignin particles in ceramic membrane pores, with the aim of elucidating the microscopic mechanisms of membrane fouling. The movement of lignin particles is predicted to be primarily governed by local hydrodynamic forces and, for sub‑micron particles, by Brownian motion, whereas van der Waals forces determine the strength of particle–particle and particle–surface adhesion during deposition. Because the magnitude of this adhesive interaction was modeled as being controlled by the ‘surface energy’ parameter in the JKR model, calibration of this parameter was essential for reliable simulation results. Accordingly, this study concentrated on systematically analyzing how ‘surface energy’ could influence coordination number, filter‑cake porosity, and deposition morphology during particle sedimentation. The analysis identified a reasonable and physically consistent range for the ‘surface energy’ parameter. The results indicated that setting the particle ‘surface energy’ between 0.2 and 1.0 J/m² yielded deposition behavior that could closely resemble experimental trends reported for lignin filtration, thereby providing a theoretical basis for more accurate prediction and regulation of membrane‑fouling behavior.

Download PDF

Full Article

Semi-Resolved CFD-DEM Study on the Influence of ‘Surface Energy’ Factors on the Deposition of Lignin Particles in Ceramic Membrane Pores

Hao Wang ,^a Junfei Wu,^b,* Ping Fu,^b Xinle Song,^c LongJing Li,^a Yibo Yang,^a Weiwu Xu,^a Zhaolei Li,^a and ZiZhen Yang ^d,*

DOI: 10.15376/biores.21.3.6537-6568

Keywords: Black liquor; Lignin; Ceramic membrane; Semi-resolved CFD-DEM; Surface free energy; Discrete element model; Computational flued dynamics

Contact information: a: College of Automotive Engineering, Zibo Polytechnic University, Zibo, Shandong Province, 255300, China; b: College of Electromechanical Engineering, Qingdao University of Science and Technology, Qingdao, Shandong Province, 266061, China; c: Zibo Polytechnic University, Zibo, Shandong Province, 255300, China; d: Inner Mongolia Metal Material Research Institute, Yantai, Shandong Province, 264003, China; *Corresponding authors: jfw_2002@qust.edu.cn; yzzno52@163.com

INTRODUCTION

Black liquor, which is a byproduct generated during the pulping and cooking stage (Sun et al. 2018), is a high-concentration aqueous organic mixture produced in the treatment process of papermaking wastewater (Ji et al. 2017). Black liquor contains a large amount of lignin and lignin derivatives (Gea et al. 2002; Zhao et al. 2019; Madyaratri et al. 2022), which account for more than 30% of the total solids (Fierro et al. 2006) and is a main by-product of the paper industry. Lignin is an amorphous three-dimensional network polymer with a complex structure and is recalcitrant to biodegradation (Wei 2019; Dominguez et al. 2020). Nevertheless, lignin has excellent physical and chemical properties and can be widely applied in various industrial productions. It is a renewable resource with a very high utilization value and is also the only non-petroleum resource in nature that can provide renewable aryl compounds (Piccitto et al. 2022). Thus, practically speaking, the separation of lignin from black liquor with low loss and high efficiency is of great significance in reducing black-liquor pollution and improving the utilization of waste resources.

The extraction and recovery of lignin have become a research hotspot in biorefinery technology. Due to its excellent chemical stability, thermal stability, strong acid and alkali resistance, high mechanical strength, and long lifespan (Colyar et al. 2008; Wang 2018), the ceramic membrane method is among the most popular and efficient methods for removing lignin from black liquor. However, during the separation process, the problem of ceramic membrane fouling has persisted, which has become the key reason why the ceramic membrane cannot be used for the continuous industrial production of lignin isolation from the black liquor of pulping. The lignin particle motion and deposition processes inside the pores of ceramic membrane are intricate and are important factors affecting the entire filtration process. However, the motion and deposition processes of lignin particles within the ceramic membrane are microscopic processes, and these processes change rapidly. The limitations of existing theories and experimental methods have confined the research on membrane fouling to only the macroscopic level, and the research mainly focuses on the macroscopic operating parameters during the filtration process (Mendes et al. 2021).

Recent experimental studies on ceramic membrane fouling provide valuable context for understanding the microscopic deposition processes relevant to this study. High‑resolution microscopic observations and post‑filtration characterization have shown that particulate suspensions often form heterogeneous, dendritic, or cluster‑based deposits rather than uniform cakes, and these structures evolve dynamically during filtration (Huang et al. 2006; Agbangla et al. 2012; Xiong et al. 2020). Experiments involving filtration of lignocellulosic or particulate systems have further demonstrated that particle aggregates may undergo rotation, bending, partial collapse, and restructuring under hydrodynamic forces during the filtration process (Dingwell et al. 2011; Xie et al. 2019). Moreover, pore‑scale investigations of filtration and particle migration highlight the important roles of local hydrodynamics, interparticle adhesion, and contact mechanics in determining fouling morphology and clogging behavior (Agbangla et al. 2012; Xie et al. 2022). These findings emphasize that membrane fouling is governed by microscopic particle interactions and dynamic deposition processes, motivating the use of numerical approaches, such as computational fluid dynamics simulation with discrete element method (CFD‑DEM), which are capable of resolving particle motion and deposition in membrane pores. With the rapid development of computer technology and the proposal of the discrete element method, the method based on CFD-DEM coupled solution has become an important means for studying the motion and deposition behaviors of lignin particles within the pores of the ceramic membrane. Such approaches can be strengthened by obtaining microscopic information. By such means there is potential to learn above the mechanisms of membrane fouling, because the CFD-DEM coupled solution can characterize the kinetic behavior of particles at the microscopic level.

The CFD-DEM coupling solution approach can microscopically depict the actual particle shape, deposition pattern, particle-fluid interaction, and particle-particle collision process during the filtration process (Fang et al. 2024). However, for micron-sized particles such as lignin with relatively small particle diameters, during the CFD-DEM coupled solution process, their motion states are influenced by forces from multiple aspects. Among these interacting forces, the van der Waals adhesion force plays a decisive role in the motion of the particles (Tao et al. 2018, 2020). Therefore, the effect of van der Waals adhesion on the motion of particles has become the focus and difficulty of current research. As early as 1971, Johnson et al. (1971) developed a JKR model according to DEM approach. The model considered the correlation between the van der Waals adhesion force and ‘surface energy’ between two contacting elastomers. Hence, the JKR model can characterize the van der Waals adhesion of the materials. In addition, an equation was derived for the influence of ‘surface energy’ on the contact size and adhesion force between the surfaces of two lightly-loaded spherical solids. Thus, in the JKR model, the magnitude of the ‘surface energy’ determines the magnitude of the van der Waals adhesion force, making it the most crucial macroscopic characterization parameter reflecting the magnitude of the van der Waals adhesion force. Dominik and Tielens (1997) along with Marshall (2009) modeled the accumulation of small particles in the channel with the JKR model. The outcomes demonstrate that the JKR model has considerable advantages in examining the deposition process of small particles. Qian et al. (2013, 2014) along with Cao et al. (2020) chose the JKR model to simulate the process of particle movement and deposition within the fibrous membrane, and they conducted research on these processes. However, in these studies, the ‘surface energy’ was set as a constant value. As a result, the van der Waals adhesion was substituted by a normal force of constant value, which led to an inability to accurately represent the adhesion between particles and between particles and the membrane (Qian et al. 2014). The approach failed to ensure the accuracy of the simulation of particle motion and deposition under this ‘surface energy’. Tao et al. (2020) examined the influence of the ‘surface energy’ between particles and between particles and fibers on the particle bridging process through using a dual-fiber filtration system. The results revealed that the ‘surface energy’ had a considerable impact on the formation and morphology of particle bridges, demonstrating that the ‘surface energy’ factors strongly affect the particle motion and deposition process. Thus, the alignment of simulation results with reality depends on correct ‘surface energy’ factors.

Accurately modeling micron-particle dynamics with variable ‘surface energy’ is computationally demanding. Unresolved CFD-DEM models, while less computationally intensive, impose mesh size constraints (3 to 5 times particle diameter) that are often incompatible with resolving flow details near membranes (Anderson and Jackson 1967; Tsuji et al. 1993; Zhou et al. 2010; Cheng et al. 2018; Wang et al. 2019, 2020; Penn et al. 2020; Xie et al. 2021). Conversely, resolved CFD-DEM models provide high fidelity by fully resolving flow around particles but require extremely fine meshes (particle diameter 8 to 10 times mesh size), leading to prohibitively high computational costs for micron-scale systems (Hager et al. 2014; Cheng et al. 2018). To overcome these computational bottlenecks for studying lignin deposition in membrane pores, this study employs a semi-resolved CFD-DEM coupling approach (Wang et al. 2021). This method reconstructs background flow information near particles using a kernel function, eliminating strict mesh size dependencies and enabling efficient simulation of micron-particle dynamics in complex geometries such as membrane pores.

Consequently, the motion and deposition processes of lignin particles inside a single pore of ceramic membrane were simulated in this work via the semi-resolved CFD-DEM model, and the impact of ‘surface energy’ factors on these processes was examined. Firstly, the movement and deposition processes of lignin particles within individual pores of ceramic membranes were dynamically characterized from a microscopic perspective. Then, analysis was conducted on the impact of the ‘surface energy’ factors on the filter cake porosity, coordination number together with deposition pattern of lignin particles. Lastly, the calibration of ‘surface energy’ factors was completed considering various factors comprehensively.

EXPERIMENTAL

Implementation of Semi-Resolved CFD-DEM Model

The key distinctions between the semi-resolved CFD-DEM model and the unresolved CFD-DEM model are manifested in two aspects: the way of calculating the void fraction of grid cells and the approach to acquiring the background flow-field information adjacent to particles. Intuitively, during the CFD-DEM coupled solution, the fluid flow area shifts from the original entire grid area to the area of the grid excluding particles. As a result, when solving the Navier-Stokes equation of the fluid, it is necessary to multiply it by the mesh void fraction, which makes the solution of the mesh void fraction extremely important. To improve clarity regarding the modeling approach, Table 1 provides a comparison of the main features, assumptions, and limitations of unresolved, semi‑resolved, and fully resolved CFD‑DEM models.

Table 1. Comparison of Unresolved, Resolved, and Semi‑resolved CFD‑DEM Models

In the unresolved CFD-DEM model, the centroid of every particle is utilized to get the label of its grid unit in the fluid-solid coupling system. For every grid, the void fraction is computed by first subtracting the total volume of all particles within the grid and then dividing the result by the volume of that grid, as shown in Eq. 1 (Sun and Xiao 2015a, b),

(1)

where, V_pi and ΔV are the volume of particle i and the current fluid grid, respectively (m³), and n is the total number of particles contained in the current grid.

When a particle is located at the junction of grid cells, that is, a single particle occupies multiple grid cells simultaneously, this approach cannot reliably represent the percentage of the particle in each grid cell, resulting in a large deviation in the calculation of the void fraction of the grid cell. In addition, when the size of the grid cell is smaller than the particle diameter, the void fraction calculated by this method will be negative, leading to calculation divergence (Xie et al. 2022).

The flow state surrounding each particle is not perfectly resolved by the semi-resolved CFD-DEM model. Rather, it introduces a kernel function to redistribute the particle phase’s volume in space and captures the physical data on the background flow field in the extended area around the particle through this kernel function (as shown in Fig. 1(b)). Then, through local volume fraction weighting, the velocity of the particle and the interaction force between the fluid and the particle are mapped to the local flow field around the particle, thus correcting the interaction force between the fluid and the particle.

Figure 1 shows schematic diagrams of particle volume distribution and background flow field in the semi-resolved along with unresolved CFD-DEM models. As can be seen from Fig. 1, only the data in cell I may be applied as background information in the conventional unresolved CFD-DEM model. When the grid size and particle diameter are not set reasonably, it will be impossible to reasonably acquire the physical information in the background flow field, resulting in inaccurate or even divergent calculation results. In contrast, the semi-resolved CFD-DEM model can reconstruct the background flow field through the kernel function and smoothly distribute the volume of the particle (as shown by the shaded part in Fig. 1(b)), thus eliminating the dependence on grid size in the coupled solution and greatly improving the calculation speed.

Fig. 1. Diagram of particle volume distribution and background flow field. (a) Diagram of unresolved model. (b) Diagram of semi-resolved model (Inspired by Xie et al. 2021)

In the semi-resolved CFD-DEM model, the background flow field information required by particle i can be given by Eq. 2 (Wang et al. 2021),

(2)

where J is the ID of grid cell in smooth area around particle i; is the background flow field information stored in grid cell J; are void fraction and volume in cell-J, respectively. N_cells is the total number of fluid grid cells in smooth area around particles; and is kernel function of particle-i on cell-J.

This model selected a smooth Gaussian function as the kernel function. The Gaussian kernel function of particle-i on grid cell-J is shown in Eq. 3 (Sun and Xiao 2015a, b),

(3)

where T is the transposition of vector, and b is Gaussian smooth distance (m). To meet the best solution accuracy, it is usually 3 to 4 times the particle radius. In this model, the Gaussian smoothing distance b is four times the particle radius.

In this model, a single particle may fill several grid cells, and the grid size may be less than the particle diameter. To avoid the multi-value problem in the solution of void fraction, the void fractions of grid cells within a distance of 4Ri from the center point of the particles was averaged to give the mean void fraction, as shown in Eq. 4 (Xie et al. 2022).

(4)

To provide context for the applicability of the present numerical model, the qualitative trends predicted by the semi‑resolved CFD‑DEM framework have been compared with experimental observations reported in the literature. Previous studies on particle deposition on single fibers and porous substrates (Huang et al. 2006; Xiong et al. 2020; Dingwell et al. 2011) have documented the formation of dendritic particle structures, the transition from initial membrane‑controlled capture to particle‑controlled capture, and the development of anisotropic deposition patterns. These experimentally observed behaviors are consistent with the morphological evolution and particle‑chain growth mechanisms predicted by the current simulation approach. Although the present work focuses on mechanistic modeling rather than reproducing a specific experimental setup, these published results support the physical realism of the deposition processes captured by the numerical model.

Although the JKR contact model provides a well‑established description of adhesive interactions, it is formally derived for contacts between spherical particles. In real black liquor, however, lignin exists not only as nearly spherical precipitated nanoparticles, but also as irregular and fragment‑like structures with highly variable morphology under alkaline conditions. Representing such structural diversity with idealized spheres is therefore a geometric simplification of the present CFD‑DEM framework. This assumption was adopted to ensure numerical stability, tractable collision detection, and compatibility with the semi‑resolved coupling scheme, all of which become significantly more complex for non‑spherical shapes. As a result, the simulated adhesion forces should be interpreted as effective, spherically averaged interactions rather than exact shape‑specific forces. Future research incorporating clumped‑particle or super‑quadric representations will more fully resolve the influence of lignin particle anisotropy on contact mechanics and deposition behavior.

Although the semi‑resolved CFD‑DEM framework improves the accuracy of local void fraction estimation and background flow‑field reconstruction, it remains a numerical model whose predictions must be interpreted in light of available experimental observations. Experimental studies of particle deposition on fibers and ceramic membranes (Huang et al. 2006; Dingwell et al. 2011; Xie et al. 2019) have consistently reported phenomena such as dendritic chain growth, particle rotation or bending under hydrodynamic loading, and partial collapse or restructuring of particle aggregates. These behaviors broadly agree with the predicted evolution of dendritic structures in the present simulations, particularly the transition from membrane‑dominated capture to secondary interception by growing particle chains. However, experiments also revealed additional effects, such as partial particle coalescence, local compaction of aggregates, and progressive restructuring in densely packed regions, which cannot be fully represented in the current model because (i) the JKR formulation treats adhesion at the particle scale without explicitly modelling material deformation or coalescence, and (ii) the semi‑resolved CFD‑DEM approach locally averages fluid stresses rather than resolving micro‑scale shear variations within porous deposits. Consequently, minor discrepancies between simulated and experimentally observed deposition thickness, local coordination number, or fine‑scale pore blockage patterns are expected. These differences do not contradict the major qualitative trends captured by the model but highlight the need for future coupling with higher‑resolution interparticle constitutive descriptions or multiscale frameworks to more closely reproduce experimentally observed deposit consolidation mechanisms.

Control Equation for Fluid–Solid Two-Phase Flow

The coupled solution process of CFD-DEM is illustrated in Fig. 2. The whole coupled solution process is performed on the basis of the Euler-Lagrange framework. The motions of the particles and the fluid are solved via the Newton’s second law and Navier-Stokes (N-S) equations, separately, and the coupling is realized via Newton’s third law.

Fig. 2. Flowchart of the CFD-DEM coupling

Fluid phase

The governing equations for the fluid phase adhere strictly to the laws of conservation of mass and conservation of momentum for locally averaged variables (Qian et al. 2014). Furthermore, in order to characterize the interaction forces between the fluid and the particles in the surrounding local flow field as well as the particle velocities, the void fraction is adopted to characterize the N – S equations. The governing equations for the fluid phase are displayed in Eqs. 5 through 7 (Chu et al. 2009),

(5)

(6)

(7)

where µ is viscosity of fluid (Pa·s), u is velocity of fluid (m/s), ƿ is density of fluid (kg/m³), S represents the total fluid drag forces F_D operating on the grid cell volume (N), and ΔV denotes the grid cell volume (m³).

Discrete mode

When particles are in motion, collisions occur not only among the particles themselves but also between the particles and the membrane. Therefore, these collisions have a significant impact on the movement of particles within the entire system, causing the particle motion during the separation process to exhibit a contact-dominated flow characteristic. Therefore, precisely capturing this process is essential for accurately characterizing the particle deposition pattern.

Among numerous models for describing particle collisions, the hard-sphere model and the soft-sphere model are two of the most fundamental ones (Tsuji et al. 1993). The hard-sphere model highly simplifies the particle contact process, regarding it as a collision behavior that can be completed instantaneously. In contrast, the soft-sphere model uses a spring-damper system to represent the dynamic process of contact between two colliding particles, as displayed in Fig. 3(a) (Cheng et al. 2022). Additionally, the contact force is determined using the soft-sphere model with the tangential displacement along with normal overlap between particles. This enables it to capture the process in which a specific particle collides with multiple surrounding particles simultaneously, thereby more realistically depicting the collisions between moving particles. Therefore, in this study, the non-linear soft-sphere model was selected to describe particle collisions.

In the framework of the soft-sphere model, the particle discrete element method determines the velocity and position of particles at every specific instant based on the external forces acting on the particles and Newton’s second law (Wang et al. 2010). The particle moving in the flow field is subjected to external forces including collision force, particle-fluid force and gravity, as depicted in Fig. 3(b). Additionally, since the particle size in this research is 1 μm, the influence of Van der Waals adhesion cannot be overlooked. The governing equations for the particle phase are shown in Eqs. 8 and 9 (El-Emam et al. 2023),

(8)

(9)

The interaction forces between particles in colloidal systems may also include electrostatic double‑layer forces arising from surface charges in the liquid medium. In black liquor systems, lignin particles typically carry surface charges due to dissociated phenolic and carboxylic acid groups, which can generate electrostatic repulsion between particles. However, black liquor is characterized by a relatively high ionic strength, which compresses the electrical double layer and significantly weakens long‑range electrostatic interactions. Therefore, in the present model the dominant interparticle attractive interaction is assumed to be van der Waals adhesion, which is represented by the Hertz–Mindlin with JKR contact model. The electrostatic double‑layer force is not explicitly included in the current simulations, and its potential influence on particle aggregation and deposition behavior should be considered in future research.

Fig. 3. (a) Diagram of soft sphere model; (b) diagram of the force in the particle’s collision process (Reprinted with permission from Elsevier, Qian et al. 2014)

The adhesion force in contact area was described through the Hertz-Mindlin with JKR model, which is a cohesive contact model that estimates the van der Waals adhesion at equilibrium in the contact region (Li et al. 2011) for ideal, smooth surfaces. Its expression is formulated in Eq. 10 as follows:

(10)

where is ‘surface energy’ (J/m²), E^* is the equivalent Young’s modulus (Pa), α denotes the normal overlap (m), and R^* represents the relative radius (m).

F_fp represents the force between the fluid and particles. Compared with the inertia of the particles themselves, the magnitudes of the buoyancy force, virtual mass force, and Saffman lift force are relatively small. Therefore, when considering the interaction forces between the liquid-solid two-phase, these forces are omitted, and only the drag force F_drag and pressure gradient force F_p, are considered. Therefore, the interaction force between the fluid and the particle is given by Eq. 11. F_drag stands for the fluid resistance that acts on the particle (Qiu et al. 2023), which is presented in Eq. 12,

(11)

(12)

where A is the particle’s projected area (m²), C_D represents the drag coefficient, and its expression is formulated in Eq. 13.

(13)

In Eq. 13, denotes the force of pressure gradient (Napolitano et al. 2022), which can be calculated by Eq. 14,

(14)

where F_c stands for the contact force between particles (Mindlin et al. 1949; Madrid et al. 2022), as presented in Eq. 15,

(15)

where F_cn is normal contact force acting on particles after particles collision (N), and F_ct denotes tangential contact force that acts on the particles upon collision (N).

Although both static and rolling friction coefficients are included in Table 2, their use in the present DEM framework does not indicate that lignin particles are assumed to undergo hydrodynamic sliding along surfaces. For spherical particles at the microscale, fluid shear primarily generates a torque rather than a net tangential force, and rolling is therefore expected to be the dominant mode of motion. In the soft‑sphere Hertz–Mindlin contact formulation, the frictional parameters regulate the tangential compliance during particle–particle and particle–membrane contact, ensuring numerical stability and allowing the model to capture partial rolling, resisted rotation, and small amounts of microslip that occur when adhesive contacts deform under load. Thus, the friction coefficients should be interpreted as effective parameters governing rotational resistance in the JKR–Hertz–Mindlin contact law, rather than as an assumption that particles translate by sliding under hydrodynamic forcing.

Although the JKR formulation captures normal adhesive interactions between deformable spherical particles, it does not explicitly model the development of a continuous lignin phase at contact points, as would be expected during multi-layer or dendritic deposition. In such cases, a local ‘healing effect’ may occur, in which adjacent particles partially merge or form a more continuous cohesive bridge. This phenomenon is beyond the microscopic resolution of the present CFD–DEM framework. Based on the simplified approach used here to represent the colloidal forces and energies, the JKR model is essentially an equilibrium theory, more suitable for describing adhesive behavior after particles have come to rest. Its use in this study to simulate dynamic deposition involves simplifying assumptions that enable the simulations but also limit the applicability of the model to some extent. Therefore, the results presented here should be interpreted as qualitative trends.

Through the ‘surface energy’ defined herein, the model can still implicitly capture the stability trends of particle bridges and dendritic structures. It should be noted that this ‘surface energy’ is actually a computed interaction energy term based on simplified assumptions, not a true material surface free energy; for consistency with the existing literature. Therefore, for studies of dynamic deposition processes, future readers are advised to adopt a combined CFD–DEM and DLVO approach rather than directly using the JKR model, as this would allow a more comprehensive description of particle interaction forces before contact equilibrium and potentially lead to more reliable predictions.

Physical Model

In practical situations, some ceramic ball-shaped particles accumulate in the ceramic membrane during its preparation and sintering to form filtration channels. This porous ceramic membrane structure can be close to the actual filtration process in the simulation. Given that the authors have successfully achieved simplified modeling of single-pore ceramic membranes in previous investigations into the influence of the rolling friction coefficient on the deposition behavior of lignin particles within these membranes (Wang et al. 2023), this study continues to adopt the methodology established by Tao et al. (2020) and Xiong et al. (2020), which simplifies multi-channel fiber membranes to single-channel fiber membranes, to construct the single-pore ceramic membrane structure. Figure 4 shows the structure of a single ceramic membrane pore. A ceramic particle sphere in 10 μm diameter was utilized to create the pore of ceramic membrane in accordance with the suspended-particle approach for the fabrication of porous ceramic membranes (Xie et al. 2019; Wang et al. 2023). The calculation range was – 10 μm < X < 10 μm, – 10 μm < Y < 10 μm, and – 30 μm < Z < 30 μm.

Fig. 4. Structure of the pores of the ceramic membrane and the size of the computational domain

Boundary Conditions and Parameter Settings

Table 2 lists the primary parameters utilized in the CFD-DEM coupled solution procedure. During the solution process, the boundary conditions of inlet velocity and outlet pressure were adopted. The inlet velocity was set to 0.5 m/s, and the outlet pressure was equal to atmospheric pressure. The particle size distribution of lignin particles in black liquor is extensive. Because simulating the motion of very small particles using the CFD–DEM coupling method is computationally expensive, the diameter of the lignin particles was set to 1 μm in the simulations. This enables a reduction in calculation time and the representation of micron-sized particles with a smaller particle size. The particle size distribution of lignin particles in black liquor is typically broad, ranging from submicron to several micrometers according to experimental observations. In the present simulation, a representative particle diameter of 1 µm was selected to reduce the computational cost associated with CFD‑DEM simulations involving a large number of very small particles. This simplified assumption allows the study to concentrate on the fundamental mechanisms governing particle transport, adhesion, and dendritic deposition behavior in membrane pores. Nevertheless, it is noteworthy that the actual deposition behavior may be influenced by the particle size distribution. Variations in particle size can affect drag forces, collision frequency, and adhesion interactions, which may further influence the morphology and stability of the dendritic deposition structures. Therefore, the results presented here represent the behavior of typical micron‑scale lignin particles, while simulations incorporating realistic particle size distributions can provide further insight in future studies. In addition, lignin particles present in black liquor exhibit significant heterogeneity in chemical composition and molecular structure, which leads to variations in surface properties, such as ‘surface energy’, charge density, and hydrophobicity. The lignin present in black liquor is chemically heterogeneous and consists of a mixture of fragments with different molecular weights, functional groups, and surface properties.

Previous studies have indicated that membrane fouling is often dominated by the fraction of the lignin mixture that exhibits the strongest tendency to adsorb or aggregate on surfaces, rather than by the average behavior of all lignin components (Colyar et al. 2008; Dingwell et al. 2011; Mendes et al. 2021). In the present study, this complexity was represented in a simplified manner by using an effective ‘surface energy’ parameter in the JKR adhesion model (Johnson et al. 1971; Li et al. 2011; Tao et al. 2018, 2020). Accordingly, lignin particles were treated as having uniform material properties, and the variability in composition was indirectly represented by varying the ‘surface energy’ parameter. The selected range of ‘surface energy’ therefore reflects the adhesion characteristics of the lignin fractions with the highest deposition propensity, allowing the model to capture the dominant fouling mechanisms while maintaining computational tractability.

Furthermore, lignin particles in practical black liquor suspensions often exhibit irregular and non‑spherical morphologies. For computational efficiency and to simplify the contact‑detection algorithm in DEM simulations, the particles are assumed to be spherical in the present model. This assumption is commonly adopted in CFD‑DEM studies and enables efficient calculation of collision and hydrodynamic interactions. Nevertheless, particle shape can influence drag forces, contact mechanics, and aggregation behavior, and future work may incorporate non‑spherical particle models to further improve the representation of realistic lignin particle morphology.

In actual kraft black liquor, lignin-derived colloidal particles are electrostatically stabilized at high pH due to the dissociation of phenolic hydroxyl groups (pKa ≈ 10) and, upon further acidification, carboxylic acid groups (pKa ≈ 3.5). These dissociated groups impart substantial negative charge, generating electrostatic repulsion that suppresses particle aggregation until the pH is lowered. This mechanism explains why progressive acidification promotes lignin precipitation and aggregation, as reviewed in Hubbe et al. (2019). In the present numerical model, such electrostatic stabilization was not explicitly resolved; instead, its macroscopic effect on aggregation propensity was represented through the ‘surface energy’ parameter in the JKR adhesion model. Accordingly, the selected range of ‘surface energy’ values should be understood as effective adhesion strengths corresponding to lignin fractions or solution conditions, in which electrostatic repulsion has been sufficiently reduced for deposition to occur.

Table 2. Simulation Parameters

Given that the numerical model employed in this study is consistent with the model previously established for investigating the influence of the rolling friction coefficient on the deposition behavior of lignin particles in single-pore ceramic membranes (Wang et al. 2023), the mesh generation scheme was adopted in the current work to ensure the reliability and comparability of the numerical simulation results. It used the hexahedral grid to divide the fluid domain. The mesh resolution ranged from 0.2 to 1 μm (46,400 total cells), giving a particle‑to‑grid size ratio between 1 and 5. The mesh quality ranged from 0.7 to 1, indicating acceptable element shape and resolution for the simulations.

RESULTS AND DISCUSSION

Analysis of the Motion and Deposition Processes of Lignin Particles

In previous CFD-DEM studies of the fiber filtration process, a ‘freezing’ model was adopted to describe the particle deposition process on fibers (Li et al. 2007). This model assumes that particles no longer move after depositing on the fibers and are ‘frozen’ in place, ignoring the movement of particles after deposition. In fact, the ‘freezing’ model was only applied to the initial phase of filtering. When the deposited particles form a dendritic structure, on the one hand, the dendritic structure continues to capture particles for growth. On the other hand, the fluid force acting on the dendritic structure can cause the dendritic structure to rotate, bend, fold, or break (Huang et al. 2006), until the dendritic structure reaches an equilibrium state. This not only seriously influences the deposition pattern of particles on the membrane, but also, due to the interlacing of dendritic structures, it significantly influences the movement of particles within the membrane pores. Therefore, in this section, the kinetic behavior of the movement together with deposition of lignin particles on ceramic membranes was investigated, encompassing the process of initial lignin particle deposition on ceramic membranes, the development of the dendritic structure, and the final movement process towards equilibrium.

Movement and deposition of lignin particles in the ceramic membrane pores were simulated under the conditions that the ‘surface energy’ between lignin particles is assumed to be 0.6 J/m² and the ‘surface energy’ between lignin particles and the ceramic membrane is assumed to be 0.5 J/m².

Figure 5 presents the dynamic variation of the motion and deposition pattern of lignin particles over time under these conditions. The arrow in the figure indicates the flow direction of the fluid. As time goes on, the quantity of lignin particles that have been deposited keeps growing. In the initial stage of filtration (t = 0 to 0.5 ms), the lignin particles gradually approach the ceramic membrane, carried by the black liquor. It is assumed that, during this period, only a few lignin particles that come into direct contact with the ceramic membrane surface are captured. Most of the remaining lignin particles exit the computational area with the black fluid after immediately passing through the membrane pores. As the filtration process advances (t = 0.5 to 1.0 ms), the lignin particles that have already been formed on the ceramic membrane surface gradually begin to be crucial to the particle-capturing process. The lignin particles are captured not just by the ceramic membrane but also by the lignin particles that were previously deposited. Under the influence of this dual-capture mechanism, during this phase, the morphology of the deposited lignin particles alters dramatically. The deposited lignin particles begin to extend outwards to form a dendritic structure, as shown in the black-framed area in Fig 6. As the filtration process proceeds further (t = 1.0 to 2.0 ms), the dendritic structure becomes more and more distinct and independent. On the one hand, the dendritic structure continuously captures lignin particles to keep growing. On the other hand, under the action of fluid forces, phenomena, including rotation or bending start to occur, elevating the contact area between the moving lignin particles and dendritic structure. As a result, the dendritic structure currently captures the majority of the lignin particles. After the injection of lignin particles ceases (t = 2.0 to 2.2 ms), the lignin particles on the dendritic structure are not replenished. At this time, the dendritic structure continues to experience bending, rotation, and other phenomena under the action of fluid forces until it reaches the final equilibrium state.

Fig. 5. Process of lignin particle deposition morphology changing with time. (, )

The changes over time of number of deposited lignin particles along with cumulative number of lignin particles passing through the membrane pores were statistically examined in order to better understand the influence of the creation of the dendritic structure on the lignin particle deposition. Figure 6 displays the findings. The slope of the curve showing the total number of lignin particles that pass through the pores of membrane steadily diminishes with time, while the slope of the curve for the number of deposited lignin particles in the membrane pores gradually increases. This occurs because, during the early stage of filtration, only a small fraction of lignin particles are directly captured by the membrane surface, while the open pores allow most particles to pass through and exit the computational domain. However, with the formation of the dendritic structure, it provides a secondary interception of lignin particles. Thus, more and more lignin particles are intercepted, leading to a gradual decrease in the number of lignin particles passing through. Due to the relatively small total number of lignin particles generated throughout the filtration process, the entire filtration process precisely lies within the clean-filtration and transitional-filtration stages. Therefore, the dendritic structure cannot fully develop to the extent of completely blocking the membrane pores, and thus the phenomenon of the membrane pores being completely blocked and particles being unable to pass through the pores cannot be observed. Furthermore, throughout the entire filtration process, a total of 286 lignin particles were trapped. However, only 46 of them were directly caught by the membrane, constituting just 16% of the overall captured particles. Additionally, as the filtration process advances, this proportion will further decrease. To sum up, in the initial stage of filtration, the ceramic membrane serves as the main way to capture lignin particles, while in the subsequent stage, it is the dendritic structure that mainly captures these particles.

Fig. 6. Curve of cumulative penetrating particle number and deposited particle number with time

Influence of ‘Surface Energy’ on Particle Deposition Pattern

During the deposition process, lignin particles form particle chains with a dendritic structure. Under the action of fluid forces, the particle chains may experience various phenomena, including toppling, bending, rotation, fracture, and collapse. This phenomenon persists until all particles reach a state of force equilibrium again or until the end of the filtration (Xiong et al. 2020). Because the setting of ‘surface energy’ directly affects the force of van der Waals adhesion between the ceramic membrane and lignin particles as well as between lignin particles, it changes the ability of the particle chains to resist fluid forces, thus influencing the deposition process of lignin particles and ultimately affecting the deposition pattern of lignin particles on the ceramic membrane. Therefore, in this section, a quantitative description of the impact of ‘surface energy’ on the lignin particle deposition pattern is provided.

Under the condition that the ‘surface energy’ between the membrane and lignin particles is constantly 0.5 J/m² , Fig. 7 displays the findings of a study on the impact of varying ‘surface energy’ among lignin particles on their deposition pattern. The figure exhibits that the deposition pattern of lignin particles is considerably influenced by the ‘surface energy’ between them. In the JKR model, the adhesion force and the torque caused by the adhesion force are proportional to (Liu 2011). Therefore, the greater the ‘surface energy’ between lignin particles, the greater the adhesion force between them. When the ‘surface energy’ between lignin particles is relatively small although lignin particles can deposit on the surface of membrane, owing to the low ‘surface energy’, the adhesion force generated between lignin particles is insufficient to resist the driving force generated by fluid flow. As a result, the dendritic structure bends and topples under the action of fluid force, and it finally reaches an equilibrium state on the downstream side of the ceramic membrane, presenting a morphology in which the dendritic structure grows along the fluid flow direction. As the ‘surface energy’ between lignin particles continues to increase the adhesion force between lignin particles also keeps strengthening. The ability of the dendritic structure to resist bending and toppling is continuously enhanced, and independent dendritic structures begin to grow on the upstream side of the ceramic membrane. In addition, the greater the ‘surface energy’ between lignin particles, the greater the contact stiffness between them. It becomes more difficult for particles to rotate around the point of contact, and the dendritic structure is also more difficult to bend. This leads to a relatively fixed position of lignin particles within the dendritic structure, and the independence of the dendritic structure is further enhanced, resulting in a depositional morphology of the whole lignin particles similar to a ‘forest’. This is also consistent with the deposition pattern of particles on a single fiber observed experimentally by Huang et al. (2006).

Fig. 7. Influence of ‘surface energy’ between lignin particles on particle deposition morphology

In addition, to further investigate the influence of the ‘surface energy’ between lignin particles on the deposition pattern of lignin particles, the variation of the number of deposited lignin particles over time under different surface energies between lignin particles was statistically analyzed, and the results are shown in Fig. 8. It can be observed that the curves under various working conditions exhibit similar trends. The reasons for such trends have been explained earlier. Here, the focus is on studying the impact of the ‘surface energy’ between lignin particles on the number of deposited lignin particles. In the first step of filtration, since the membrane is still primarily responsible for capturing lignin particles, there is no discernible difference in the number of deposited lignin particles under different working conditions. As the filtration process progresses, the number of deposited lignin particles begins to vary with different surface energies. The number of deposited lignin particles does not increase continuously with the increase of the ‘surface energy’ between lignin particles. When the ‘surface energy’ between lignin particles is less than 0.2 , the number of deposited lignin particles increases as the ‘surface energy’ between lignin particles increases. However, when the ‘surface energy’ between lignin particles is greater than 0.2 , the number of deposited lignin particles decreases as the ‘surface energy’ increases. This is caused by the two-sided effect of ‘surface energy’ on deposition. Taking into account the lignin particles’ deposition pattern, when the ‘surface energy’ between the lignin particles is small, due to the weak adhesion between lignin particles, the dendritic structure is more likely to topple towards the membrane surface under the action of fluid force. The dendritic structure is closer to the membrane surface, and its structural stability is better, thus trapping more lignin particles under a lower ‘surface energy’. When the ‘surface energy’ between lignin particles is large, the dendritic structure is in a state of independent growth. The deposition pattern of lignin particles indicated that the number of lignin particles in the roots of the dendritic structures near the membrane surface was relatively small. In contrast, the lignin particles are carried away from the flow field by the fluid force.

Fig. 8. Influence of ‘surface energy’ between lignin particles on the number of particle deposits with time

Figure 9 demonstrates the effect of ‘surface energy’ between lignin particles and membrane (0.05, 0.1, 0.5, and 1.0 J/m²) on the deposition pattern of the lignin particles at a constant ‘surface energy’ of 0.5 J/m² between lignin particles. As Fig. 9 indicates, the ‘surface energy’ between lignin particles and the membrane has an insignificant impact on the deposition pattern of the lignin particles. The degree of particle aggregation, the length of the dendritic structure, and the counter-current growth situation are similar. This is because, as described in earlier studies, the lignin particles captured directly by the membrane are too small in number.

Fig. 9. Influence of ‘surface energy’ between lignin particles and membrane on particle deposition morphology

Figure 10 illustrates how ‘surface energy’ affects the quantity of lignin particles that are only captured by the ceramic membrane. The results suggest that the number of lignin particles that the membrane is able to catch is independent of the ‘surface energy’ between the particles. Nevertheless, when the ‘surface energy’ between the membrane and the lignin particles grows, so does the amount of lignin particles that the membrane is able to trap. However, the quantity of lignin particles that the membrane is able to catch is insufficient when compared to the total number of lignin particles that were collected. Therefore, although the ‘surface energy’ between lignin particles and the membrane has a certain influence on the number of lignin particles captured by the membrane, it has no effect on the overall deposition pattern of lignin particles. Therefore, one of the main factors affecting the lignin particle deposition pattern is the ‘surface energy’ between the particles.

Fig. 10. Influence of ‘surface energy’ on the number of particles captured by the membrane. (a) Influence of ‘surface energy’ between lignin particles on the number of particles captured by the membrane. (b) Influence of ‘surface energy’ between lignin particles and membrane on the number of particles captured by the membrane

In addition, to deeply explore the influence of the ‘surface energy’ between lignin particles on the deposition pattern of the lignin particles, the distribution of the deposition angles of lignin particles under different surface energies between lignin particles was statistically analyzed. Figure 11 (a) presents the definition criteria for the deposition angles of lignin particles. Based on these criteria, the distribution range of the deposition angles of lignin particles was obtained, and the findings are displayed in Fig. 11 (b).

Fig. 11. Particle deposition angle distribution. (a) Schematic showing the angle criterion. (b) Influence of ‘surface energy’ between lignin particles on deposition angle distribution of lignin particles

Figure 11(b) shows the distribution of the particle deposition angle when the ‘surface energy’ between the lignin particles and the film is , with the ‘surface energy’ between the lignin particles being 0.02, 0.1, 0.5, and 1.0 J/m² respectively. As can be seen from the figure, regardless of the magnitude of the ‘surface energy’ between lignin particles, the deposition of lignin particles on the membrane is approximately symmetrically distributed. There are about 75.8% of the particles deposited in the range of β = -90 ° to β = 90 ° as . This indicates that the majority of lignin particles are deposited on the upstream side of the membrane, with only a few located on the downstream side. However, as the ‘surface energy’ decreases, the number of lignin particles that were deposited on the upstream side of the membrane progressively reduced. When the lignin particles deposited in the range of β = -90 ° to β = 90 ° decrease to 73.4%, and the lignin particles deposited in the range of β = -90 ° to β = 90 ° are only 61.1 % and 62.2% as By considering the deposition pattern of lignin particles, it is not difficult to realize that the growth process of the particle chain greatly affects the deposition angle of lignin particles. When the ‘surface energy’ among lignin particles is high, the particle chains’ ability to resist bending and rotation is also strong. The counter-current growth of the particle chains can lead to most particles being on an upstream side of membrane. When the ‘surface energy’ among lignin particles is low, the particle chain is unable to resist the action of the fluid force. It begins to bend or move along the fluid flow direction, causing a rising number of particles to be on the downstream side of the membrane. In addition, due to the characteristics of fluid flow, there are not many lignin particles deposited in the range of β = -30 ° to β = 30 °, and most of the lignin particles are in the range of β = -30 °~ -90 °and β = 30 ° ~ 90 °.

Influence of ‘Surface Energy’ on Coordination Number of Particles

The coordination number, i.e., the number of particles that are in contact with the center particle, serves as the most sensitive parameter reflecting the local microstructure of the particle (Yang et al. 2000). Generally, the minimum coordination number plays a vital role in keeping the particle chain (the force chain) continuous. For 1 μm particles, the most appropriate minimum coordination number is 2 (Yang et al. 2000), as can be seen from the last particles of the particle chain within the black dashed box in Fig. 5. At this moment, a particle, while supporting other particles, is itself also supported by other particles. Nevertheless, the local coordination number will be more prominent at the positions where particle chains are interwoven and particles overlap. Consequently, within this section, the effect of ‘surface energy’ on both the average coordination number and the distribution of the coordination number of lignin particles was explored.

Figure 12 presents the curve depicting the impact of ‘surface energy’ on the mean coordination number of the lignin particle deposition structure. As the ‘surface energy’ between lignin particles continuously increases, the average coordination number also increases, with the rate of increase becoming increasingly slow. When the ‘surface energy’ between lignin particles rises from 0.01 to 1.0 J/m², the average coordination number increases from 2.03 to 3.69. Although an increase in ‘surface energy’ increases the inter‑particle adhesion, leading to the formation of more independent branches in the dendritic structure – which may reduce the average coordination number – it also increases the normal overlap between lignin particles (see Table 3), which tends to increase the coordination number. The net effect is that the average coordination number gradually increases with ‘surface energy’. The acceptable minimum mean coordination number for particles with a particle size of 1 μm by gravity alone is 2.13, according to a study on the distribution of the coordination number of the deposited particle packing structure under the sole action of gravity (Yang et al. 2000). Therefore, solely in terms of meeting the requirement of the minimum average coordination number, the acceptable minimum ‘surface energy’ between lignin particles is 0.02 J/m². In addition, when the ‘surface energy’ between lignin particles remains constant, the average coordination number fluctuates irregularly with changes in the ‘surface energy’ between lignin particles and the membrane. This indicates that the ‘surface energy’ between lignin particles and the membrane has no significant impact on the average coordination number. The reason, again, is that the number of lignin particles captured by the membrane is too small.

Fig. 12. Relation curve of average coordination number with ‘surface energy’

Table 3. Influence of the ‘Surface Energy’ between Lignin Particles on the Normal Overlap Between Particles

Figure 13 shows the influence of ‘surface energy’ on the distribution of coordination numbers. The ‘surface energy’ between the membrane and lignin particles has no significant impact on the distribution of coordination numbers. The primary element influencing the distribution of coordination numbers is still the ‘surface energy’ between lignin particles. As the ‘surface energy’ between lignin particles gradually increases, the normal overlap between lignin particles also continuously increases, thus causing the distribution range of coordination numbers to gradually expand.

Fig. 13. Influence of ‘surface energy’ on coordination number distribution. (a) Influence of ‘surface energy’ between lignin particles on coordination number distribution. (b) Influence of ‘surface energy’ between lignin particles and membrane on coordination number distribution

Influence of ‘Surface Energy’ on Porosity

Porosity () is a significant parameter utilized to characterize the properties of the particle deposition structure. The porosity of the lignin particle deposition structure is also influenced by the ‘surface energy’ setting. The porosity magnitude can be computed from Eq. 16 (Xiong et al. 2020),

(16)

The porosity profile of the lignin-particle deposition structure calculated by Eq. 16 changing with ‘surface energy’ is shown in Fig. 14. The ‘surface energy’ between lignin particles and the membrane has no effect on the porosity. The porosity only gradually decreases as the ‘surface energy’ between lignin particles increases. The porosity drops from 0.91 to 0.66 as the ‘surface energy’ among the lignin particles rises from 0.01 to 1.0 J/m². This is due to the fact that as the ‘surface energy’ among the lignin particles rises, the normal overlap between the lignin particles also rises, which leads to a reduction in the porosity. The porosity of the filter cake formed by lignin particles is between 0.61 and 0.71 (Dingwell et al. 2011). When the ‘surface energy’ between lignin particles is in the range of 0.2 to 1.0 J/m², the porosity is between 0.66 and 0.71. Therefore, it can be determined that during the simulation process, setting the ‘surface energy’ between lignin particles in the range of 0.2 to 1.0 J/m² is more in line with the actual porosity situation.

Fig. 14. The variation curve of porosity with ‘surface energy’

CONCLUSIONS

In this study, a semi-resolved CFD-DEM coupling model was used to simulate the movement of lignin particles in a single ceramic membrane pore, with all simulations based on the simplified approach used here to represent colloidal forces and energies. It should be noted that the JKR model is essentially an equilibrium theory strictly applicable to adhesive behavior after particles have come to rest; its use for simulating dynamic deposition therefore has theoretical limitations. Nevertheless, to obtain exploratory insights into trends such as dendritic deposition morphology under the available computational constraints, the model was adopted as a simplified tool. The ‘surface energy’ defined herein is not a true material surface free energy, but rather a computed interaction energy term based on simplifying assumptions. Under this premise, the effects of ‘surface energy’ on filter cake porosity, particle coordination number, and deposition morphology were systematically examined. The results indicate that this simplified approach provides usable effective parameters, and the simulated trends qualitatively approach the real filtration process.
The ‘surface energy’ between particles is a critical factor influencing the deposition morphology and deposition angle. The larger the ‘surface energy’, the more independent the particle chains become, resulting in a more pronounced ‘forest’ morphology and an increased number of particles depositing upstream.
The ‘surface energy’ between particles is also a critical factor influencing the porosity and coordination number of the deposition structure. In contrast, the ‘surface energy’ between particles and the membrane has a relatively negligible impact. As the ‘surface energy’ between particles increases, the mean coordination number rises from 2.03 to 3.69, while the porosity decreases from 0.91 to 0.66.
Considering the effects of the inter-particle ‘surface energy’ on porosity, coordination number, and deposition morphology, it is recommended to set the ‘surface energy’ in the range of 0.2 to 1.0 J/m² in similar simulations, which makes the simulated trends closer to the actual lignin filtration process. Future work should attempt to combine CFD-DEM with the DLVO model, which can take into consideration the interaction forces before the particles come into equilibration upon contact, thereby providing a more mechanistic basis for predicting membrane fouling rather than relying solely on an empirical ‘surface energy’ parameter.

ACKNOWLEDGMENTS

The authors gratefully acknowledge Prof. Wu Junfei and Prof. Fu Ping from Qingdao University of Science and Technology for their generous support of this study. This study was funded by Shandong Provincial Natural Science Foundation (Grant No. ZR2023MB086) and the Science and Technology Support Plan of Youth Innovation Team of Shandong Higher School (Grant No. 2024KJH150).

REFERENCES CITED

Agbangla, G. C., Climent, É., and Bacchin, P. (2012). “Experimental investigation of pore clogging by microparticles: Evidence for a critical flux density of particle yielding arches and deposits,” Separation and Purification Technology 101, 42-48. https://doi.org/10.1016/j.seppur.2012.09.011

Anderson, T. B., and Jackson, R. O. Y. (1967). “Fluid mechanical description of fluidized beds. Equations of motion,” Industrial & Engineering Chemistry Fundamentals 6(4), 527-539. https://doi.org/10.1021/i160024a007

Cao, B., Wang, S., Dong, W., Zhu, J., Qian, F., Lu, J., and Han, Y. (2020). “Investigation of the filtration performance for fibrous media: Coupling of a semi-analytical model with CFD on Voronoi-based microstructure,” Separation and Purification Technology 251, article 117364. https://doi.org/10.1016/j.seppur.2020.117364

Cheng, K., Wang, Y., and Yang, Q. (2018). “A semi-resolved CFD-DEM model for seepage-induced fine particle migration in gap-graded soils,” Computers and Geotechnics 100, 30-51. https://doi.org/10.1016/j.compgeo.2018.04.004

Cheng, K., Zhu, J., Qian, F., Cao, B., Lu, J., and Han, Y. (2022). “CFD–DEM simulation of particle deposition characteristics of pleated air filter media based on porous media model,” Particuology 72, 37-48. https://doi.org/10.1016/j.partic.2022.02.003

Chu, K. W., Wang, B., Yu, A. B., and Vince, A. (2009). “CFD-DEM modelling of multiphase flow in dense medium cyclones,” Powder Technology 193(3), 235-247. https://doi.org/10.1016/j.powtec.2009.03.015

Colyar, K. R., Pellegrino, J., and Kadam, K. (2008). “Fractionation of pre-hydrolysis products from lignocellulosic biomass by an ultrafiltration ceramic tubular membrane,” Separation Science and Technology 43(3), 447-476. https://doi.org/10.1080/01496390701812517

Dingwell, K., Sedin, M., and Theliander, H. (2011). “Filtration of lignin from a lignocellulose-based ethanol pilot plant,” Environmental Engineering Science 28(11), 775-779. https://doi.org/10.1089/ees.2010.0397

Domínguez, E., del Río, P. G., Romaní, A., Garrote, G., Gullón, P., and de Vega, A. (2020). “Formosolv pretreatment to fractionate paulownia wood following a biorefinery approach: Isolation and characterization of the lignin fraction,” Agronomy 10(8), article 1205. https://doi.org/10.3390/agronomy10081205

Dominik, C., and Tielens, A. G. G. M. (1997). “The physics of dust coagulation and the structure of dust aggregates in space,” The Astrophysical Journal 480(2), article 647. https://doi.org/10.1086/303996

El-Emam, M. A., Zhou, L., and Omara, A. I. (2023). “Predicting the performance of aero-type cyclone separators with different spiral inlets under macroscopic bio-granular flow using CFD–DEM modelling,” Biosystems Engineering 233, 125-150. https://doi.org/10.1016/j.biosystemseng.2023.08.003

Fang, W., Wang, X., Zhu, C., Han, D., Zang, N., and Chen, X. (2024). “Analysis of film unloading mechanism and parameter optimization of air suction-type cotton plough residual film recovery machine based on CFD–DEM coupling,” Agriculture 14(7), article 1021. https://doi.org/10.3390/agriculture14071021

Fierro, V., Torné-Fernández, V., and Celzard, A. (2006). “Kraft lignin as a precursor for microporous activated carbons prepared by impregnation with ortho-phosphoric acid: Synthesis and textural characterisation,” Microporous and Mesoporous Materials 92(1-3), 243-250. https://doi.org/10.1016/j.micromeso.2006.01.013

Gea, G., Murillo, M. B., and Arauzo, J. (2002). “Thermal degradation of alkaline black liquor from straw: Thermogravimetric study,” Industrial & Engineering Chemistry Research 41(19), 4714-4721. https://doi.org/10.1021/ie020283z

Hager, A., Kloss, C., Pirker, S., and Goniva, C. (2014). “Parallel resolved open source CFD-DEM: Method, validation and application,” The Journal of Computational Multiphase Flows 6(1), 13-27. https://doi.org/10.1260/1757-482X.6.1.13

Huang, B., Yao, Q., Li, S. Q., Zhao, H. L., Song, Q., and You, C. F. (2006). “Experimental investigation on the particle capture by a single fiber using microscopic image technique,” Powder Technology 163(3), 125-133. https://doi.org/10.1016/j.powtec.2006.01.014

Hubbe, M., Alén, R., Paleologou, M., Kannangara, M., and Kihlman, J. (2019). “Lignin recovery from spent alkaline pulping liquors using acidification, membrane separation, and related processing steps: A review,” BioResources 14(1), 2300-2351. https://doi.org/10.15376/biores.14.1.2300-2351

Ji, X. (2017). Experimental Study and Numerical Simulation on Straw Pulp Black Liquor Combustion in Fluidized Bed for Direct Causticization, Ph.D. Dissertation, Harbin Institute of Technology, Harbin, China.

Johnson, K. L., Kendall, K., and Roberts, A. A. D. (1971). “ ‘Surface energy’ and the contact of elastic solids,” Proceedings of the Royal Society of London A: Mathematical and Physical Sciences 324(1558), 301-313. https://doi.org/10.1098/rspa.1971.0141

Li, S. Q., and Marshall, J. S. (2007). “Discrete element simulation of micro-particle deposition on a cylindrical fiber in an array,” Journal of Aerosol Science 38(10), 1031-1046. https://doi.org/10.1016/j.jaerosci.2007.08.004

Li, S., Marshall, J. S., Liu, G., and Yao, Q. (2011). “Adhesive particulate flow: The discrete-element method and its application in energy and environmental engineering,” Progress in Energy and Combustion Science 37(6), 633-668. https://doi.org/10.1016/j.pecs.2011.02.001

Liu, G. (2011). Discrete Element Methods of Fine Particle Dynamics in Presence of van der Waals and Electrostatic Forces, Ph.D. Dissertation, Tsinghua University, Beijing, China.

Madrid, M. A., Fuentes, J. M., Ayuga, F., and Gallego, E. (2022). “Determination of the angle of repose and coefficient of rolling friction for wood pellets,” Agronomy 12(2), article 424. https://doi.org/10.3390/agronomy12020424

Madyaratri, E. W., Iswanto, A. H., Nawawi, D. S., Lee, S. H., and Fatriasari, W. (2022). “Improvement of thermal behavior of rattan by lignosulphonate impregnation treatment,” Forests 13(11), article 1773. https://doi.org/10.3390/f13111773

Marshall, J. S. (2009). “Discrete-element modeling of particulate aerosol flows,” Journal of Computational Physics 228, 1541-1561. https://doi.org/10.1016/j.jcp.2008.10.035

Mendes, S. F., Rodrigues, J. S., de Lima, V. H., Botaro, V. R., Cardoso, V. L., and Reis, M. H. (2021). “Forward black liquor acid precipitation: Lignin fractionation by ultrafiltration,” Applied Biochemistry and Biotechnology 193(10), 3079-3097. https://doi.org/10.1007/s12010-021-03580-2

Mindlin, R. D. (1949). “Compliance of elastic bodies in contact,” in: Mathematical Theory of Elasticity, Springer. https://doi.org/10.1007/978-1-4613-8865-4_24

Napolitano, E. S., Renzo, A. D., and Maio, F. P. D. (2022). “Coarse-grain CFD-DEM modelling of dense particle flow in gas-solid cyclone,” Separation and Purification Technology 287, article 120591. https://doi.org/10.1016/j.seppur.2022.120591

Penn, A., Padash, A., Lehnert, M., Pruessmann, K. P., Müller, C. R., and Boyce, C. M. (2020). “Asynchronous bubble pinch-off pattern arising in fluidized beds due to jet interaction: A magnetic resonance imaging and computational modeling study,” Physical Review Fluids 5(9), article 094303.

Piccitto, A., Scordia, D., Corinzia, S. A., Cosentino, S. L., and Testa, G. (2022). “Advanced biomethane production from biologically pretreated giant reed under different harvest times,” Agronomy 12(3), article 712. https://doi.org/10.3390/agronomy12030712

Qian, F., Huang, N., Zhu, X., and Lu, J. (2013). “Numerical study of the gas-solid flow characteristic of fibrous media based on SEM using CFD-DEM,” Powder Technology 249, 63-70. https://doi.org/10.1016/j.powtec.2013.07.030

Qian, F., Huang, N., Lu, J., and Han, Y. (2014). “CFD-DEM simulation of the filtration performance for fibrous media based on the mimic structure,” Computers & Chemical Engineering 71, 478-488. https://doi.org/10.1016/j.compchemeng.2014.09.018

Qiu, M., Chen, Z., Jiang, L., Liu, R., Tang, Y., and Liu, M. (2023). “Numerical simulation of uranium tetrafluoride fluorination in a multistage spouted bed using the improved CFD-DEM chemical reaction model,” Particuology 75, 119-136. https://doi.org/10.1016/j.partic.2022.07.010

Sun, M., Wang, Y., Shi, L., and Klemeš, J. J. (2018). “Uncovering energy use, carbon emissions and environmental burdens of pulp and paper industry: A systematic review and meta-analysis,” Renewable and Sustainable Energy Reviews 92, 823-833. https://doi.org/10.1016/j.rser.2018.04.036

Sun, R., and Xiao, H. (2015a). “Diffusion-based coarse graining in hybrid continuum–discrete solvers: Theoretical formulation and a priori tests,” International Journal of Multiphase Flow 77, 142-157. https://doi.org/10.1016/j.ijmultiphaseflow.2015.08.014

Sun, R., and Xiao, H. (2015b). “Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD-DEM,” International Journal of Multiphase Flow 72, 233-247. https://doi.org/10.1016/j.ijmultiphaseflow.2015.02.014

Tao, R., Yang, M., and Li, S. (2020). “Effect of adhesion on clogging of microparticles in fiber filtration by DEM-CFD simulation,” Powder Technology 360, 289-300. https://doi.org/10.1016/j.powtec.2019.09.083

Tao, R., Yang, M. M., and Li, S. Q. (2018). “Filtration of micro-particles within multi-fiber arrays by adhesive DEM-CFD simulation,” Journal of Zhejiang University–SCIENCE A 19(1), 34-44. https://doi.org/10.1631/jzus.A1700156

Tsuji, Y., Kawaguchi, T., and Tanaka, T. (1993). “Discrete particle simulation of two-dimensional fluidized bed,” Powder Technology 77(1), 79-87. https://doi.org/10.1016/0032-5910(93)85010-7

Wang, G., Hao, W., and Wang, J. (2010). Discrete Element Method and Its Practice on EDEM, Northwestern Polytechnical University Press, Xi’an, China.

Wang, H., Wang, X., Wu, Y., Wang, S., Wu, J., Fu, P., and Li, Y. (2023). “Study of CFD-DEM on the impact of the rolling friction coefficient on deposition of lignin particles in a single ceramic membrane pore,” Membranes 13(4), article 382. https://doi.org/10.3390/membranes13040382

Wang, W. P. (2018). The Preparation of Inorganic Microporous Membrane and Application in Lignin Separation, Master’s Thesis, South China University of Technology, Guangzhou, China.

Wang, Z., and Liu, M. (2020). “Semi-resolved CFD-DEM for thermal particulate flows with applications to fluidized beds,” International Journal of Heat and Mass Transfer 159, article 120150. https://doi.org/10.1016/j.ijheatmasstransfer.2020.120150

Wang, Z. K., and Liu, M. B. (2021). “Whole-process cross-scale modelling of laser direct deposition with semi-resolved VOF-DEM coupling,” Chinese Journal of Theoretical and Applied Mechanics 53(12), 3228-3239.

Wang, Z., Teng, Y., and Liu, M. (2019). “A semi-resolved CFD-DEM approach for particulate flows with kernel-based approximation and Hilbert curve-based searching strategy,” Journal of Computational Physics 384, 151-169. https://doi.org/10.1016/j.jcp.2019.01.017

Wei, Z. M. (2019). Study on the Correlation between Lignin Characteristic Structure and COD, Master’s Dissertation, South China University of Technology, Guangzhou, China. https://doi.org/10.27151/d.cnki.ghnlu.2019.001111

Xie, C., Ma, H., and Zhao, Y. (2019). “Investigation of modeling non-spherical particles by using spherical discrete element model with rolling friction,” Engineering Analysis with Boundary Elements 105, 207-220. https://doi.org/10.1016/j.enganabound.2019.04.013

Xie, Z., Wang, S., and Shen, Y. (2021). “CFD-DEM modelling of the migration of fines in suspension flow through a solid packed bed,” Chemical Engineering Science 231, article 116261. https://doi.org/10.1016/j.ces.2020.116261

Xie, Z., Wang, S., and Shen, Y. (2022). “Particle-scale modelling of rapid granular filtration in a dual-media filter,” Separation and Purification Technology 302, article 122076. https://doi.org/10.1016/j.seppur.2022.122076

Xiong, G., Gao, Z., Hong, C., Qiu, B., and Li, S. (2020). “Effect of the rolling friction coefficient on particles’ deposition morphology on single fibre,” Computers and Geotechnics 121, article 103450. https://doi.org/10.1016/j.compgeo.2020.103450

Yang, R. Y., Zou, R. P., and Yu, A. B. (2000). “Computer simulation of the packing of fine particles,” Physical Review E 62(3), article 3900. https://doi.org/10.1103/PhysRevE.62.3900

Zhao, Q. J., Ding, S., Wen, Z. M., and Toppinen, A. (2019). “Energy flows and carbon footprint in the forestry-pulp and paper industry,” Forests 10, article 725. https://doi.org/10.3390/f10090725

Zhou, Z. Y., Kuang, S. B., Chu, K. W., and Yu, A. (2010). “Discrete particle simulation of particle-fluid flow: Model formulations and their applicability,” Journal of Fluid Mechanics 661, 482-510. https://doi.org/10.1017/S002211201000306X

Article submitted: December 8, 2025; Peer review completed: March 21, 2026; Revisions accepted: May 16, 2026; Published: May 29, 2026.

DOI: 10.15376/biores.21.3.6537-6568