A novel semi-resolved CFD-DEM method with two-grid mapping:methodology and validation

The semi-resolved CFD-DEM method has emerged as a prominent tool for modelling particle-ﬂuid interactions in granular materials with high particle size ratios. However, challenges arise from conﬂicting requirements regarding the CFD grid size, which must adequately resolve ﬂuid ﬂow in the pore space while maintaining a physically meaningful porosity ﬁeld. This study addresses these challenges by introducing a two-grid mapping approach. Initially, the porosity ﬁeld associated with ﬁne particles is estimated using a coarse CFD grid, which is then mapped to a dynamically reﬁned grid. To ensure conservation of total solid volume, a volume compensation procedure is implemented. The proposed method has been rigorously validated using benchmark cases, showing its high computational eﬃciency and accurate handling of complex porosity calculations near the surface of coarse particles. Moreover, the previously unreported impact of the empirical drag correlation on ﬂuid-particle force calculations for both coarse and ﬁne particles has been revealed.


Introduction
Fluid-particle flows are common in various industrial sectors, and gaining insight into their behaviour is crucial for developing new downstream processes and formulations [1,2], as well as optimising existing ones.Additionally, comprehending the coupled response of fluidparticle systems is vital in geomechanics applications, particularly in cases such as internal erosion [3,4], where the overall system deformation may be relatively small.Computational Fluid Dynamics coupled with the Discrete Element Method (CFD-DEM) is a numerical tool that can simulate two-phase fluid-particle systems.It is a Lagrangian-Eulerian approach [5,6], in which the fluid is modelled as a continuum phase and the particle phase is treated as a series of discrete elements; the simulation results give particle-scale resolution as the trajectory of each particle is traced.CFD-DEM is usually divided into two categories, namely resolved and unresolved, according to the resolution of the CFD solver.Resolved CFD-DEM treats the particle surfaces as no-slip boundary conditions, the immersed boundary method [7,8] or fictitious domain method [9,10] is used to simulate the fluid flow field, and thus a very fine CFD mesh should be used.In contrast, unresolved DEM-CFD uses the homogenised Navier-Stokes equations, the presence of the particles is represented by the porosity 1 or solids fraction, so that the particle-fluid interface is not resolved, and a minimum cell size of 1.6-3.0times the particle diameter is required in the simulations [11,12].
In the case of bidisperse, bimodal, or gap-graded particle systems with a relatively large size ratio, there are technical challenges associated with use of both unresolved and resolved methods.If the unresolved method is applied, the first challenge is to select a suitable drag force expression to determine the fluid-particle interaction force.While empirical models for mono-disperse particles [13,14] and poly-disperse systems [15] have been used in many studies of the segregation of poly-disperse particle systems in fluidised beds, their reliability is questionable for large size ratios.The second issue in applying the method is determining the porosity (or solids fraction), which is also an important parameter required to calculate the drag force.It is difficult to calculate this accurately, or to obtain a representative value, when the CFD mesh size approaches the size of the largest particle diameter [16].
These problems become more pronounced the larger the particle size ratio of the system.
While resolved CFD-DEM implementations do not require specification of a drag model, this approach has to use a mesh with a grid size that is less than 1/10 of the diameter of the smallest particle in the system [16].This fine resolution requirement significantly increases the computational cost and renders resolved CFD-DEM simulations unsuitable for most industrial applications.Some measures have been applied to address this issue.For example, Tsuji et al. [17] proposed a fictitious particles method in which large particles are made up of smaller fictitious particles; a numerical calibration procedure is required to determine the fictitious particle size and the fictitious volume fraction.
A recently proposed compromise solution is semi-resolved CFD-DEM [18,19], in which the coarse particle and fine particle fractions are simulated using resolved and unresolved schemes, respectively.Note that the term "semi-resolved" CFD-DEM may also be used to refer to a method to smooth out the porosity field as in the works of Wang et al. [20,21], Xie et al. [22].Here, when the term "semi-resolved CFD-DEM" is used we refer to a formulation in which the empirical drag model is applied to determine the fluid-particle interaction force for the fine particles only, while the fluid-particle interaction forces for the coarse particles are calculated analytically based on the fluid flow field.In this way the predictive capacity is improved compared to the unresolved CFD-DEM and the restriction on particle size ratio is removed.However, determining the porosity field for the calculations associated with the finer particles is not straightforward due to the discrepancy between the required CFD mesh sizes for the unresolved and resolved methods.In other words, the fine mesh applied to resolve the flow around the coarse particles may yield non-physical porosity values for the fine particle domain.Currently, a solution for this issue reported in the literature is applying the Gaussian-kernel weighting function to distribute the volume of a single particle to the adjacent CFD cells [19,23].This approach was originally proposed for the unresolved method [24].However, the limitations of such a method are evident.The width of the kernel is arbitrarily determined, and there are difficulties in dealing with porosity calculations nearwalls (and near-coarse particle surfaces), making them inflexible.A recent study by Xie et al. [22] proposed a hybrid CFD-DEM solver that uses fictitious domain and unresolved methods to simulate fluid-particle interactions for grid-to-particle size ratios of < 0.1, 0.1 − 3.0, and > 3.0; this is in fact an improved approach to the semi-resolved method.However, the method also includes a porosity model that employs a Gaussian kernel function.Notably, Xie et al. [22] differed from Yang et al. [19] in that they only applied the kernel-based porosity model to medium sized-particles.This work presents a two-grid semi-resolved CFD-DEM solver, in which a coarse mesh and a fine mesh are overlain.The coarse mesh is introduced to calculate the porosity field for the fine particles, while the fine mesh (refined from the coarse mesh) is used to estimate the porosity field for the coarse particles.This method is inspired by the two-grid approaches that have been applied in unresolved CFD-DEM simulations [25,11,26], which demonstrated the flexibility of a two-grid formulation.The proposed method is validated using existing simulation results obtained using a resolved CFD-DEM method [16].The effects of mesh refinement level and the empirical drag force model adopted on the simulation accuracy are explored.
The paper is organised as follows: in Section 2, the formulation of the semi-resolved CFD-DEM is presented; in Section 3, the porosity calculation method using two mesh grids is shown; the approach to calculate and post-process the fluid-particle interaction forces is given in Section 4 followed by an introduction to the simulation setup in Section 5.The performance of the proposed solver is evaluated in terms of the solid volume conservation, accuracy of the fluid-particle interaction force prediction, etc. in Section 6; Finally, the conclusions and the aspects for future improvement are laid out in Section 8.

Governing equations for the semi-resolved CFD-DEM
The governing equations of the semi-resolved solver can be derived by combining the governing equations of the resolved and unresolved solvers.First, we briefly introduce the schemes adopted in each of these two solvers.

Resolved solver
The CFD approach used in the resolved solver is known as the fictitious domain method [27,28], which is the general form of the immersed boundary method (IBM).The success of this approach hinges on the ability to enforce a no-slip boundary condition on the particle surface; this can be achieved through various methods such as correcting the velocity field or applying a momentum source term.Additionally, the flow field within the solid object is also solved.The equations for the resolved solver assume incompressible fluid flow.The continuity and momentum equations are where ρ f and u f are the density and velocity of fluid, and g is acceleration due to gravity.
The most important condition is that the following no-slip conditions should be fulfilled at the fluid particle interfaces: where Ω p is the region occupied by the particle or the so-called fictitious domain, and u l is the local particle velocity in a CFD cell, given by: where u p and ω p are the translational and the angular velocity of the particle.
A resolved solver has been implemented in the official release in CFDEM [29], with details provided in [10].The procedure to implement the no-slip boundary is relatively complicated: the fluid velocity in the region occupied by the coarse particles is corrected directly based on the particle velocity field (u l ) and the gradient of a correction factor ψ, from which the resulting velocity field is calculated as: In order to make u f divergence-free , i.e., ∇ • u f = 0, ψ fulfils the following condition: In addition to correcting the velocity field, Hager [10] also suggested expanding the pressure p by the term ∂ψ ∂t to match the momentum equation.However, the algorithms mentioned above are not capable of achieving no-slip boundary conditions.In fact, Eq. 5 shows that there is a noticeable velocity difference between the fluid and particles within the fictitious domain (Ω p ).It has been found that a good solution to overcome this limitation is direct application of the body force term S reso pf , estimated based on the velocity difference between the particle and fluid, to achieve a no-slip boundary condition.This approach has been successfully employed in previous studies, such as [9,30].
The momentum equation (Eq.2) is rewritten as where the S reso pf is calculated by where u f is fluid velocity obtained without considering the presence of the particles, and S reso pf is non-zero only in regions occupied by particles (represented by a porosity ε reso f < 1) and thus eliminates the need for explicit calculation of the Lagrangian multiplier.We validate this method in Section 6.1.

Unresolved solver
The equations describing the behaviour of the fluid phase in the unresolved solver use the volume averaged Navier-Stokes equations, which were originally derived by Anderson & Jackson [31].The presence of the particle phase in the fluid is represented by a porosity field (ε unreso f ), and empirical drag models are employed to calculate the fluid-particle force.The continuity equation and momentum equations [5,6] are written as The momentum source term, S unreso pf , in Eq. 10 is estimated as where G pf is the momentum source coefficient.G pf is calculated as the sum of the drag forces (obtained through empirical drag models) of the particles within a CFD cell, as detailed below.
The porosity field exists in both the resolved and unresolved solvers and it is defined as the volume fraction of the fluid (gas or liquid) in a CFD cell as follows: where V s and V cell are the volume of solid material and the CFD grid cell, respectively.Fig. 1 is a schematic diagram of the porosity fields (ε reso f and ε unreso f ).The value of ε reso f is 0 (inside the particle) or 1 (outside the particle).The porosity is calculated for each cell in the CFD mesh.As one CFD cell contains multiple particles, the value of ε unreso f of the mono-disperse particles is estimated to be in the range of [0.37, 1.0] (based upon the data in [32]).
Along with the calculation of the ε unreso f in the unresolved solver, the contribution of a particle i to the total solid volume in a specific CFD cell j is also recorded and expressed as a weight term W i,j = V p,i,j V s,tot,j , which will be used to interpolate the particle proprieties to the CFD cells.This is the reason that the porosity calculation is also called "averaging" or "interpolation" [18].The interpolation scheme that maps the solid particle velocities to the CFD cells is given by The momentum source coefficient G pf,j associated with the particles in Eq. 11 is obtained by interpolation to give:

Semi-resolved solver
The idea behind the semi-resolved solver is to extend the governing equation of the unresolved solver with a source term, S coarse pf , as in the resolved solver, to give it the functionality of a resolved solver.The resulting equations may then be transformed into resolved or unresolved formats, depending on whether coarse or fine particles occupy the regions.This is feasible since coarse and fine particles cannot physically coexist in the same area.Referring to the above equations,the governing equations for the fluid phase in a semi-resolved solver (where flow around the coarse particles is resolved and and flow around the fine particles is not resolved) can be expressed as follows: The momentum contribution from the fine particles (S f ine pf ) is determined from the drag forces on the particles: G pf is a coefficient determined by the drag force (F d,i ) on the particles, which is estimated by the empirical correlations.The contribution of the coarse particles to the momentum, S coarse pf , is calculated according to [9,33] as where u f is an intermediate velocity field determined without considering the immersed (coarse) particle.
The governing equations for the DEM solver are based on Newton's second law of motion and involve calculating the net force acting on each particle in the system, taking into account both contact and non-contact (fluid-particle interaction) forces.The translational and angular velocity of the particle is calculated by where u p,i is the velocity of particle i, N pp is the number of adjacent particles, F ∇p,i , F ∇•τ,i , and F d,i are the pressure gradient force, viscous force and the drag force, respectively, F c,ij is the inter-particle contact force which is calculated by the Hertz-Mindlin model [34,35].

Dynamic mesh refinement and porosity calculation
As mentioned earlier, difficulties still exist in calculating ε f ine f using semi-resolved CFD-DEM.To address this issue and improve the flexibility of the porosity calculation, a two-grid method was developed.This method utilises two CFD grids to calculate the porosity fields; the coarse (and static) grid is used to calculate ε f ine f .The coarse grid mesh size is 1.6 times greater than the fine particle diameter, which is in line with most methods proposed to determine the porosity such as the particle centre divided method (DM) [29].The second grid, namely the dynamically refined grid used for calculating ε coarse f and solving the governing equations, is generated using the dynamicRefineFvMesh feature in OpenFOAM.As shown in Fig. 2, the strategy is to refine the mesh close to the coarse particle surface and at the same time allow as little mesh distortion as possible.It does not change the shape of the mesh cells, rather it performs topological refinements at the region where 0 The mesh refinement level (RL) is defined relative to the starting background mesh 2 .For example, if the background mesh size is 2 mm, the mesh sizes at one and two levels of refinement will be 1 mm and 0.5 mm respectively.
The advantage of the dynamic refinement is that it achieves a high resolution where necessary around the particles while maintaining high computational efficiency.The ε coarse f field on the fine mesh grid can be estimated by a couple of methods [9,36,10] and we adopted the method proposed by Kempe et al. [36], which is based on the signed-distance level set function of the particle surface and is easy to implement.Due to the dynamic refinement process, the size of the refined CFD grid may be smaller than the fine particles in the corresponding regions.Hence, it is not possible to estimate ε f ine f directly on this mesh.
Therefore, a field mapping procedure is conducted to transfer the variables (ε f ine f , u s , G f ine pf ) from the coarse grid to the fine grid, and is introduced in the following section.

Field mapping with solid volume compensation
The mapFields tool 3 already available in OpenFOAM is used to perform the field mapping.Since the two grids have identical boundary conditions, and the refined mesh cells are obtained by splitting cells on the coarse mesh, the mapping option used is Nearest.This method searches for the cell in the coarse grid that is closest to the target cell in the fine grid and uses the value of the this closest cell directly in the target cell.
A difficulty in mapping the porosity field between two grids is that the volume of the 2 https://www.openfoam.com/documentation/guides/latest/doc/guide-meshing-snappyhexmeshcastellation.html 3 https://openfoamwiki.net/index.php/MapFieldsparticles may not be conserved, meaning that the particle volume stored in the CFD mesh may be underestimated or overestimated.Figure 5 illustrates the field mapping procedures for ε f ine f from the coarse mesh to the fine mesh.Since ε f ine f was initially estimated in the coarse mesh (without any refinement), the values outside a coarse particle may expand into its occupied region (see Fig. 5e ), where the value of ε f ine f should be one.To avoid such errors, a solid volume compensation procedure is implemented to correct the ε f ine f values.
The basic idea of the method is to set all the cell values in the internal area of the coarse particle to one and then compensate for the "lost" volume of fine particles by decreasing the porosity of the cells at the boundary of that coarse particle (see Fig. 5f-g ).
Compared to existing porosity models for semi-resolved CFD-DEM [19], the proposed two-grid method has two key advantages.Firstly, it is highly flexible as any existing porosity calculation algorithms can be used to calculate the porosity field on the coarse grid without the need for additional treatment, and the same procedure can be applied for field mapping.
Secondly, the proposed method does not require additional arbitrary input parameters such as the bandwidth in the Gaussian kernel method [19,24] or the size of the porous sphere [37] or porous cube [38].
Along with the field mapping of ε f ine f , the momentum source field G f ine pf , and the solid velocity u s are also mapped.Specifically, G f ine pf was calculated using the same weight function, based on proportion of the particle volume in each CFD cell, that was adopted to calculate ε f ine f .As a result, a correction must be applied to G f ine pf , as illustrated in Fig. 5, to ensure conservation of the total momentum.

Two force calculation schemes
The fluid-particle interaction forces acting on the coarse particles can be calculated using two different methods: the direct method [9,33] and the Shirgaonkar method [39].The direct method uses the value of the force (or momentum source term) applied in Equation 16, which is given by In the Shirgaonkar method, the total fluid-particle interaction force is estimated by integrating the force over the domain occupied by the coarse particle, given by As the fluid field around the boundary of a coarse particle is resolved in the solver, according to the fluid field with a high resolution, the resulting torque due to the fluidparticle interaction can be calculated by

Fluid-particle interaction force normalisation
In the following text,F pf,i is chosen for comparison with the fully resolved simulation because it is easily obtained from both solvers and is directly related to particle motion.
To facilitate better visualization of the fluid-particle interaction force, all the particle-fluid interaction forces discussed in the following are normalised by the Stokes force as

Simulation setup
The semi-resolved solver was implemented in the CFD-DEM open-source code CFDEM [29], which couples the CFD solver OpenFOAM [40] and the DEM solver LIGGGHTS [29].
Fig. 3 shows the flow chart of the implemented algorithm.
The semi-resolved solver outlined here is suitable for bi-modal or gap-graded particle systems with a large size ratio.In order to rigorously evaluate the accuracy and efficiency of the solver, previously obtained resolved simulation results from Knight [16] were used for a bi-modal particle assembly with a particle size ratio of 4.0.Table 1 displays the overall configuration of Knight's simulation cases, where particle assemblies with a fine particle volume fraction (f f ine ) ranging from 0.11-0.51were selected.The coupling interval, i.e., the period in which the two solvers exchange information, was set to 100∆t DEM .
The procedures for performing the immersed boundary simulations have been outlined in previous contributions [16,41].Here, we provide only a high-level description.Initially, particles were randomly placed within cubic periodic boundaries to create the samples.Then, the sample was subjected to increasing isotropic compression up to an effective stress of 100 kPa in the DEM solver using servo-controlled periodic boundaries [42].During laminar flow, the fluid-particle interactions were determined by IBM simulations using the Multiflow code [43,8,44].
Fig. 4 shows the particle assembly considered, with a particle size ratio of 4. In the CFD solver, the boundaries in the y and z direction are periodic.In the x -direction, the inlet boundary has a fixed fluid velocity and zero pressure gradient, and the outlet has a fixed pressure of zero.Periodic boundaries are also applied in the DEM solver.The sample is fixed in space and the particles do not move during the simulation.Table 2 lists the particle properties and numerical settings in the semi-resolved CFD-DEM simulations.As the "divided" porosity model is used to estimate ε f ine f , the coarse mesh size should be greater than the fine particle diameter [29].Thus the mesh-particle size ratio (SR) of 1.6 was chosen to generate the background coarse mesh.In order to capture the pore fluid flow around the coarse particle, the mesh near the boundary of the coarse particles is dynamically refined with a refinement level of 2 to 4. The simulations were terminated when the pressure drop across the particle assembly reached a stable level; then the total fluid-particle interaction force (F pf,i ) was extracted from the particle data directly.

Validation of the resolved solver: flow through ordered packings
As mentioned in section 2.1, the resolved solver now includes a new FD (or IBM) method to improve accuracy.To ensure the accuracy of the semi-resolved method based on IBM, this new method must be fully validated.Validation cases involve ordered packings of monodisperse spheres at low Reynolds numbers.Analytical solutions for these packings, including periodic simple cubic (SC), body-centred cubic (BCC), and face-centred cubic (FCC) arrays, have been provided by Zick and Homsy [45], and these data have been used by several researchers for benchmark validation purposes [46,47,48].Fig. 6 shows the variation of F pf,i with ϕ for SC, BCC and FCC arrays predicted by analytical results [45], IBM simulations (taken from [41] with the radius retraction parameter of 0.2 and D/∆x = 64) and the current work (D/∆x = 64), where D is the particle diameter and ∆x is the grid size.Both sets of simulation data almost coincide with the analytical results, with a slight discrepancy in cases with high solid fraction.Knight et al. [41] discussed the suitability of the IBM for predicting particle-fluid interaction of ordered packings.Although the accuracy of the IBM method is somewhat lower than some of the more refined calculation methods such as unstructured mesh methods, it strikes a good balance between computational cost and accuracy.

Solid volume conservation
The field mapping approach with particle volume compensation was evaluated by comparing the total solid volume of the fine fraction in the CFD mesh before and after applying the compensation.The results are shown in Fig. 7, where the y-axis represents the ratio between the total solid volume in the CFD mesh and the real solid volume of the particle assembly.There is a relatively high error ( 13%) in the solid volume for the cases that do not use the solid volume compensation algorithm.However, after applying the modification, these errors are reduced to lower than 0.1%.This indicates that the solid volume compensation algorithm is necessary and accurate for ensuring the accuracy of the field mapping approach.
Fig. 8 shows the distributions of ε f ine f and ε coarse f with different refinement levels on a slice through the domain in the y-z plane at a distance of 0.33 cm from the inlet.Figs. 9 and 10 show the variation in ε coarse f and ε f ine f along a straight line (as indicated in Fig. 8).
As Fig. 10 indicates, due to the low resolution in the coarse CFD grids, the ε f ine f field cannot represent the fine particle distribution correctly without the dynamic mesh refinement.As the dynamic cell refinement level increases, the porosity fields at the boundary region of the coarse particle sharpen and the curves become closer, in other words the local porosity at the edge of the coarse particles is captured more accurately as RF increases.The difference between them becomes less pronounced when RL ≥ 3.

Particle-fluid interaction force
Fig. 11 compares the total fluid-particle interaction forces calculated from the resolved and the semi-resolved methods for refinement levels (RL) of 2, 3 and 4 considering both the fine and coarse particles.Table 3 shows the corresponding Pearson correlation (P C) between the fluid-particle interaction forces calculated from these two methods.As RL increases from 2 to 3, the P C coarse increases by 0.144 to 0.746; as RL further increases to 4, P C coarse increases by approximately 0.04 to 0.781, which is a minor change compared with the former case.In contrast, the change in the P C f ine is not pronounced under different RL values.Therefore, mesh refinement can improve the accuracy of F coarse pf,i until a specific RL has been reached, but the accuracy of the F f ine pf,i is not affected in this process.This result is not surprising for two reasons.Firstly, the refinement of the CFD meshes helps to improve the accuracy of the fictitious domain method, as it is recommended that the mesh size should be smaller than 1/10 of the particle diameter for the resolved solver; secondly, the F f ine pf,i is mainly determined by the empirical drag model adopted, and the mesh density is less important.

Effect of coarse particle force model
The F coarse pf,i calculated using the above mentioned two force calculation methods (see Equations 21 and 22) is compared in Fig. 12.Though Shirgaonkar's method results in slightly lower force values, the differences are minor.The two force models are derived using different ideas and show close results, which again reflects the reliability of the present solver.

Effect of the empirical drag model
The proposed semi-resolved solver calculates the drag force on the fine particles using an empirical model.However, the calculated drag forces impact the flow field and so the sensitivity of both F coarse pf,i and F f ine pf,i to the drag model must be considered.In order to evaluate the effect of the drag expression adopted on the simulation accuracy, the empirical models of Ergun [49], Di Felice [50], Tang [51] and Tenneti [52] were applied in the simulation cases with bi-modal samples.These drag expressions were chosen as they have previously been applied to similar systems [16,53].Fig. 13 compares the calculated fluid particle interaction force with values from the resolved solvers; Fig. 14 shows the mean and the standard deviation (SD) values of F pf,i , and the Pearson correlation coefficient (P C) of the corresponding case is also listed in Table 4.For each case, the data for the fine and coarse particles are shown separately.After a careful analysis of those data, the following conclusions emerged.
(i) For cases with f f ine = 11%, 25% and 50%, adopting the Di Felice, Ergun and Tang models led to the best match for F coarse pf,i data, and Ergun, Tang and Tang models provided the best match for the F f ine pf,i data, respectively.(ii) In all cases, the values of P C f ine are relatively low (< 0.62) and show more scatter compared to the IBM data (with lower SD), which indicates that the stability of this approach to accurately predict F f ine pf,i is low.This result is not surprising as the empirical correlations were used to calculated the drag force (and thus the F pf,i ) for the fine particles in each case.Previous studies [41,53] have indicated that though the empirical drag correlations (for a mono-disperse particle system) can provide a good approximation to the overall drag force, they deliver a poor performance in the fluid-particle interaction force prediction for individual particles compared to the IBM simulations.In contrast, the fluid flow is resolved around the surface of the coarse particles, the F coarse pf,i was calculated directly and thus P C coarse are higher than P C f ine .We emphasise that the semi-resolved solver is mainly proposed to improve the feasibility and accuracy of the interaction between the coarse particle and the fluid flow.

Computational efficiency
From a computational standpoint, an essential criterion for any semi-resolved CFD-DEM implementation is its ability to reduce the computational cost compared to fully resolved cases.Resolved simulations for bi-modal particle assemblies, for example, typically require 10 6 to 10 7 fluid cells and 12 to 48 hours to reach steady-state flow conditions on 100 to 1000 processing units, as the mesh cell size should be smaller than one-tenth of the smallest particle diameter in the flow system.Such a high computational cost confirms that resolved solvers are not suitable for real industrial simulations.In contrast, all the semi-resolved simulations presented in this paper were run using 10 5 to 10 6 fluid cells (under mesh refinement levels of 3 to 4, see Table 3), requiring 2 to 10 hours to reach steady-state flow conditions on 4 to 12 processing units.The computational cost can be roughly estimated based on the number of CFD cells, as the calculation loads are mainly spent on the CFD solver side.The mesh cell number is only 1% to 10% of that in the resolved method.Accordingly, the semi-resolved solver roughly reduces the computational cost by more than one-tenth.The increase in computational efficiency is mainly due to the flexible algorithm for the particle resolution as well as the dynamic mesh refinement, which significantly reduces the number of CFD mesh cells.

Data Availability and Reproducibility Statement
The numerical data from Figures 6-7 and 9-14 have been tabulated in the Supplementary Material.The source code of the fully-resolved CFD-DEM solver, cfdemSolverIBPICI, along with the validation cases, is provided as a .zipfile in the Supplementary Material.

Conclusion
A two-grid semi-resolved CFD-DEM approach has been proposed, which is suitable for particle-fluid flow scenarios with a high particle size ratio, such as bi-modal particle systems.
The semi-resolved CFD-DEM method provides a trade-off between computational efficiency and accuracy, with its range of applicability falling somewhere between fully and unresolved solvers.
1. Using two CFD grids is an effective solution to address the challenge of balancing the requirements of resolving the fluid flow in the pore space and maintaining a realistic porosity field for the fine particles.This approach is easy to implement and provides great flexibility.Additionally, it has been found that a solid volume compensation procedure can ensure total solid volume conservation in the system.
2. In the case of bi-modal distributed particle systems, which is the main application scenario of this semi-resolved solver, it was found that while the proposed solver provided a highly accurate estimation of the fluid-particle interaction force for the coarse particles, the accuracy of the force prediction for each individual fine particle was relatively low.This was due to the limitations of the empirical drag models and this is a fundamental issue with any unresolved approach to CFD-DEM coupling.However, the accuracy can still be improved by applying the empirical model which performs best for the particular range of volume fraction of fine particles in the system.
3. While the computational cost of the semi-resolved CFD-DEM method depends on various factors such as the volume fraction of coarse particles and the size of the computational domain, our simulations have demonstrated that for a typical scenario involving a particle assembly with a size ratio of 4, the use of two sets of CFD grids can lead to a reduction of the total computational cost by about one order of magnitude.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :
Figure 1: Schematic diagram of the porosity fields in resolved and unresolved CFD-DEM solver

Figure 13 :Figure 14 :
Figure 13: Comparison of the F pf,i from fully resolved and semi-resolved solvers illustrating sensitivity to the drag model adopted for the finer particles (inset gives magnified image of data for the finer particles)

Table 2 :
Parameters of the particle assemblies for the semi-resolved CFD-DEM validation (coarse-fine particle size ratio was 4).

Table 1 :
Parameters of the particle assemblies for the semi-resolved CFD-DEM validation; note that the coarse-fine particle size ratio was 4 in all cases.