Analytical Modeling of Transverse and Longitudinal Motion of Single Particles in a Horizontal Boundary Layer with Cross-Flow Velocity Pulsations

1. Introduction

The study of particle-laden flows in both turbulent and laminar regimes has long attracted attention due to its significance in industrial processes, environmental engineering, and energy systems. Classical fluid dynamics provides the essential framework for understanding continuous media. Early works on fluid and gas mechanics [1], hydrodynamics [2], and mathematical physics [3] established fundamental principles. These theories are crucial for describing particle motion within fluids and the interactions between dispersed and continuous phases [4].

Particle behaviour in nonuniform and turbulent flows is governed by complex mechanisms, including drag, lift, and inertial effects. Pioneering studies of small rigid spheres in nonuniform flows highlighted the critical role of particle size and density relative to the surrounding medium [5]. Contemporary research has extended these insights to examine particle transport in laminar boundary layers [6] and to develop analytical methods for suspended particle dynamics [7]. Detailed experimental and numerical studies have revealed transversal and longitudinal particle motion, highlighting resonance phenomena that influence deposition patterns and trajectory characteristics [8,9,10].

Advanced modelling approaches now employ matrix and vector formulations. These methods provide a comprehensive representation of multidimensional particle dynamics and enable the derivation of closed-form solutions through harmonic decomposition [11,12,13,14,15]. They have been applied to predict deposition, re-suspension, and agglomeration in wall-bounded turbulent flows [16], and to refine Large-Eddy Simulation models for particle-laden systems [17,18].

The impact of particle and thermal inertia on transport and heat transfer in non-isothermal turbulent flows has been emphasised in recent computational studies [19,20]. Moreover, the spatial distribution and deposition of cylindrical, spherical, and charged inertial particles have been investigated to assess the effects of turbulence intensity, boundary layer characteristics, and particle-particle interactions [21,22,23]. Such studies are particularly relevant for process intensification strategies in chemical and biomass reactors, where optimal particle behaviour directly affects process efficiency and product yield [24,25,26].

Particle deposition and resuspension have been thoroughly investigated in turbulent channel and pipe flows using both Lagrangian and Eulerian [27,28] approaches. Boelens et al. developed Lagrangian models to analyse particle deposition and re-suspension under turbulent conditions, providing quantitative predictions of particle-wall interactions. Balachandar and Eaton [29] reviewed dispersed multiphase flows, highlighting key mechanisms governing particle transport and turbulence modulation. Hetsroni et al. [30] examined particle dispersion in wall-bounded turbulent flows, elucidating the effects of boundary layer structures on particle trajectories.

Historical studies remain critical references for understanding particle motion. Morsi and Alexander [31] offered seminal insights into particle paths in two-phase flow systems, establishing baseline equations for trajectory prediction in both laminar and turbulent regimes. Elghobashi [32] extended these concepts to particle-laden turbulent flows, emphasising the role of particle-fluid interactions in transport dynamics. Crowe and Sommerfeld [33] provided a detailed overview of particle-laden flows, integrating classical theory with practical modelling approaches, demonstrating the engineering relevance of multiphase flow comprehension. Boivin et al. [34] applied direct numerical simulations to analyse detailed particle behaviour in turbulent channels, quantifying deposition, lift-off, and near-wall interactions. Basset and Mindlin [35] introduced the concept of history forces on spherical particles, revealing cumulative effects of fluid acceleration on particle motion. Van Dop and Beij [36] investigated particle deposition and lift-off in turbulent boundary layers, demonstrating subtle interactions between particle inertia and near-wall turbulence structures.

Building on these historical foundations, the stability and qualitative behaviour of particle–fluid dynamical systems—particularly under impulsive, harmonic, or weakly stochastic excitation—can be rigorously analysed within a matrix-based dynamical framework. When the governing equations are expressed in compact matrix–vector form, they allow systematic assessment of transient and steady-state regimes through eigenvalue spectra, state-transition matrices, and transform-based operator methods. The classical Lyapunov stability theory, originating from Lyapunov’s foundational monograph [37] and further developed in modern nonlinear dynamics and matrix differential systems [39,40,41,42], provides rigorous criteria for boundedness, asymptotic stability, and resonance amplification in multidimensional particle motions. These analytical tools are especially powerful for systems subjected to transverse velocity pulsations, where small perturbations in the state vector may yield sustained oscillatory responses or transient divergence. Matrix formulations, combined with Laplace-transform techniques [43,44] and criteria for impulsive and periodically driven systems [45,46,47], therefore offer an indispensable framework for identifying stability boundaries, resonance conditions, and long-term particle behaviour—phenomena that often remain inaccessible to direct numerical integration alone. Additional extensions of these methods include impulsive differential models, coupled matrix differential equations, and stochastic operator approaches [48,49,50,51,52], which provide versatile tools for capturing complex particle–fluid interactions in laminar and turbulent boundary layers.

These studies provide a robust foundation for analysing particle transport in wall-bounded flows and motivate the use of compact state-space and transform-based methods to separate transient relaxation, forced oscillations and long-time drift. Building on this background, the present work develops a coupled analytical model for a single inertial particle in a horizontal laminar boundary layer subjected to harmonic transverse velocity pulsations. The governing equations are formulated in a matrix–vector form and solved in closed form using Laplace-transform techniques, enabling explicit identification of decay rates, phase shifts and resonance conditions. The analytical results are evaluated numerically for representative light and heavy particles, with full datasets provided in the Supplementary Materials to ensure reproducibility.

2. Theoretical Framework, Physical Model and Governing Equations

When single solid impurities move in the boundary layer of a plane horizontal flow of an incompressible fluid, numerous forces act upon them: aerodynamic drag force; added-mass inertial force; Basset force; Saffman force; Magnus force; pressure-gradient force; turbophoretic force; photophoretic force; electrophoretic force; and the body (mass) force. The magnitude and direction of each of these forces are determined primarily through empirical relations, and for every specific problem, only those forces whose influence is dominant are taken into account in the equations describing the particle velocity and trajectory.

The flow under consideration is a viscous fluid moving in a horizontal plane channel with half-width b. The coordinate system is defined as follows: the X-axis is oriented along the flow direction and coincides with the plane where the velocity is zero (the lower bounding wall); the Y-axis is directed vertically upward, perpendicular to the streamlines; and the Z-axis is directed perpendicular to the plane of the schematic.

The velocity field of the carrier phase consists of one constant longitudinal component U_g and one transverse pulsation component varying according to a harmonic law v′(t):

• along the x-axis: U_g= U_g (у),
• along the y-axis: v_g′(t)=v₀cos(ωt),

where:

- v₀ and ω are the amplitude and the frequency of the pulsations, identical for all points in the flow plane within the boundary layer.

The longitudinal carrier-phase velocity is described by a wall-bounded boundary-layer profile:

$$U_g\left(y\right)=U_{\infty }\left(1-e^{-{\displaystyle \frac{y}{\delta }}}\right), у\ge 0,$$

(1)

where:

- δ [m]- denotes the characteristic boundary-layer thickness,

If extended linearly for y < 0 using the wall shear rate, the longitudinal velocity profile becomes:

$$U_g\left(y\right)={\left.{\displaystyle \frac{d{U}_g}{dy}}\right|}_{y=0} y={\displaystyle \frac{U_{\infty }}{\delta }}y={\gamma }_{\omega }.y$$

(2)

Accordingly, the wall shear rate γ_ω is defined as:

$${\gamma }_{\omega }={\left.{\displaystyle \frac{d{U}_g}{dy}}\right|}_{y=0}={\displaystyle \frac{U_{\infty }}{\delta }}, \left[{\gamma }_{\omega }\right]=s^{-1}$$

(3)

The physical model, coordinate system, velocity components, and forces acting on a single particle in the horizontal channel boundary layer are shown schematically in Figure 1.

An important observation must be noted regarding the projections of the Saffman force: the inclusion of the transverse velocity pulsation changes the direction of the projection of the force on the x-axis. The Saffman force Fs along the y-axis has a negative sign. This follows from the algebraic operations within the brackets in the vector equation for this force:

$$F_s=6.44d_p^2\sqrt{\mu {\rho }_p}{\left|\nabla {\times V}_g\right|}^{-0.5}\left(V_g-V_p\right)\times \left(\nabla {\times V}_g\right)$$

(4)

For the wall-bounded velocity profile U_g(y), the carrier-phase vorticity is dominated by the wall-normal gradient of the longitudinal velocity. Consequently, the magnitude of the vorticity entering the Saffman force can be expressed as:

$$\left(\nabla {\times V}_g\right)\approx \left|{\displaystyle \frac{d{U}_g}{dy}}\right|={\gamma }_{\omega }$$

(5)

establishing a direct physical link between the Saffman lift mechanism and the wall shear rate.

where:

- d_p — particle diameter;

- μ — dynamic viscosity of the carrier fluid;

- ρ_p — density of the particle material;

- V_g = (U_g(y),v_g′(t),0)— velocity vector of the carrier fluid;

- V_p = (u_p,v_p,0)— velocity vector of the particle;

- ∇×V_g — vorticity vector of the carrier flow;

- ∣∇×V_g∣ — magnitude of the flow vorticity;

- (V_g−V_p) — relative velocity between fluid and particle.

Taking into account the above considerations when formulating the equations of motion, the following system is obtained for the particle velocities: u_p(t) - longitudinal and v_p(t) -transverse:

$${\displaystyle \frac{d{u}_p}{dt}}=A\left(U_g\left(y_p\right)-u_p\right)-B\left(v^{\text{'}}-v_p\right)$$

(6)

$${\displaystyle \frac{d{v}_p}{dt}}=A\left(v^{\text{'}}-v_p\right)-B\left(U_g\left(y_p\right)-u_p\right)-g\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)$$

(7)

where:

- A — drag/inertia coefficient, representing the influence of the carrier fluid along the motion direction:

$$A={\displaystyle \frac{18\mu }{{\rho }_p{d}_p^2}}$$

(8)

- B – Saffman/vorticity coefficient, representing the coupling between longitudinal and transverse motion due to fluid shear:

$$B={6.44d}_p^2\sqrt{\mu {\rho }_p}{\left|\nabla \times V_g\right|}^{-0.5}\approx {6.44d}_p^2\sqrt{\mu {\rho }_p}{\gamma }_{\omega }^{-0.5}$$

(9)

- ε_p – particle-to-fluid density ratio:

$${\epsilon }_p={\displaystyle \frac{{\rho }_p}{{\rho }_f}}$$

(10)

- ρ_p – particle density;

- ρ_f – density of the carrier fluid.

- g – acceleration due to gravity - g(1−1/ε_p) accounts for the effective body force on the particle in the vertical direction, i.e., the net gravitational acceleration reduced by the buoyant force (Archimedes’ principle).

The initial conditions corresponding to the physical process and required for solving the system are:

$$t=0;y_p\left(0\right)=y_0;u_p\left(0\right)=U_g\left(y_0\right);v_{p0}=0;$$

2.1. Transverse Motion: Governing Equations and Analytical Solution

The solution begins with the definition of the function v_p(t). Equation (7) is solved according to the difference (U_g-u_p), whereby transverse dynamics can be expressed as:

The closed-form analytical solutions derived in this section are exact for constant U_g and remain valid as a local approximation under the near-wall linearisation U_g(y)≈γ_ωy. For the full boundary-layer profile U_g(y_p(t)), the coupled first-order system is integrated numerically.

The solution begins with the definition of the function v_p(t). Equation (7) is solved with respect to the difference (U_g−u_p), yielding:

$$-B\left(U_g-u_p\right)={\displaystyle \frac{d{v}_p}{dt}}+\mathrm{A}\left(v^{\text{'}}-v_p\right)-g\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)$$

(11)

Differentiating (7) with respect to time gives:

$$B{\displaystyle \frac{d{u}_p}{dt}}=-{\displaystyle \frac{d^2{v}_p}{d{t}^2}}+\mathrm{A}\left({\displaystyle \frac{d{v}^{\text{'}}}{dt}}-{\displaystyle \frac{{dv}_p}{dt}}\right)$$

(12)

Substituting into equation (6) yields the second-order differential equation for transverse motion:

$${\displaystyle \frac{d^2{v}_p}{d{t}^2}}+2\mathrm{A}{\displaystyle \frac{{dv}_p}{dt}}+\left(A^2+B^2\right)v_p+\mathrm{A}\mathrm{g}\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)$$

This equation describes the translational and oscillatory motion of a solid particle in a viscous medium under the action of forced external disturbances in the form of harmonic velocity pulsations of the surrounding fluid. The solution of Equation (13) is composed of the homogeneous solution corresponding to damped natural oscillations (left-hand side) and a particular solution that accounts for the influence of the imposed external pulsations on the particle motion.

The homogeneous solution of (13) is:

$$v_p^{\left(h\right)}\left(\mathrm{t}\right)=\left(C_1\cos \left(Bt\right)+C_2\sin \left(Bt\right)\right)e^{-At}-{\displaystyle \frac{Ag}{A^2+B^2}}\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)$$

Given the harmonic form of the imposed velocity pulsations, the particular solution of the governing differential equation is sought in a harmonic form with the same angular frequency. Due to the linearity of the system and the sinusoidal nature of the forcing term, the steady-state response of the particle is also oscillatory, differing from the excitation by its amplitude and phase shift. Therefore, the particular solution is assumed in the form:

$$\eta =H_{1}cos\left(\omega t+{\delta }_1\right)$$

(15)

H₁ and δ₁ are constants determined by the system parameters and the characteristics of the external forcing. The constant H₁, while δ₁ denotes the phase shift between the particle motion and the imposed harmonic pulsations.

where:

- H₁ - represents the amplitude of the steady-state oscillatory response:

$$H_1={\displaystyle \frac{\left(A^2+B^2\right)v_0}{\left(A^2+B^2-{\omega }^2\right)cos\left({\delta }_1\right)-2.A\omega sin\left({\delta }_1\right)}}$$

(16)

- δ₁ - represents the phase shift accounting for the lag or lead of the particle response relative to the imposed harmonic excitation:

$${\delta }_1=arctan\left({\displaystyle \frac{-A\omega \left(A^2+B^2+{\omega }^2\right)}{{\left({A^2+B}^2\right)}^2+\left(A^2-B^2\right){\omega }^2}}\right)$$

(17)

The complete solution for the transverse velocity is obtained by superposition of the homogeneous solution (Equation (14)) and the particular solution (Equation (15)).

$$v_p=\left(C_{1}cos\left(Bt\right)+C_{2}sin\left(Bt\right)\right)e^{-At}-{\displaystyle \frac{Ag}{A^2+B^2}}\left(1-{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\epsilon }_p$}\right.}\right)+H_{1}cos\left(\omega .t+{\delta }_1\right).$$

(18)

The integration constants C₁ and C₂ are determined from the following initial conditions:

$$t=0; v_p=0; {\displaystyle \frac{d{v}_p}{dt}}=0$$

fter performing the necessary algebraic manipulations, the following expression is obtained:

$$C_1={\displaystyle \frac{Ag}{A^2+B^2}}\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)-H_{1}cos\left({\delta }_1\right)$$

(19)

$$C_2={\displaystyle \frac{A}{B}}{C}_1+{\displaystyle \frac{H_1\omega }{B}}sin{\delta }_1$$

(20)

In conclusion, Equation (17) fully characterizes the transverse particle displacement by capturing transient decay, forced oscillations, initial-condition effects, and gravitational contributions through the density ratio ε_p:

$$\begin{array}{l}y\left(t\right)=\int v_p\left(t\right)dt={\displaystyle \frac{e^{-At}}{A^2{+B}^2}}\left[\left(B{C}_1-A{C}_2\right)\sin \left(Bt\right)-\left(BA{C}_1+B{C}_2\right)c\mathit{os}\left(Bt\right)\right]\\ +\frac{H_1}{\omega }\left[sin\left(\omega t+{\delta }_1\right)-sin{\delta }_1\right]+\frac{A{C}_1-B{C}_2}{A^2{+B}^2}-\frac{A_g}{A^2{+B}^2}\left(1-\frac{1}{{\epsilon }_p}\right)t\end{array}$$

(21)

2.2. Longitudinal Motion: Governing Equations and Analytical Solution

The closed-form solution below is derived under the locally frozen approximation U_g(y_p)≈U_g(y₀). In the numerical evaluation, U_g is updated at each time step as U_g(y_p(t)).

The differential equation governing the longitudinal particle velocity as a function of time, u_p(t), is derived from the coupled system of Equations (6) and (7). By rearranging the first equation of the system and solving it with respect to the term B(v′−v_p), the following relation is obtained:

$$B\left(v^{\text{'}}-v_p\right)=-{\displaystyle \frac{d{u}_p}{dt}}+A\left(U_g-u_p\right)$$

(22)

Differentiating Equation (18) with respect to time yields:

$$\left({\displaystyle \frac{d{v}^{\text{'}}}{dt}}-{\displaystyle \frac{d{v}_p}{dt}}\right)=-{\displaystyle \frac{d^2{u}_p}{{dt}^2}}-A{\displaystyle \frac{d{u}_p}{dt}}$$

(23)

Substituting into equation (3) yields the second-order differential equation for longitudinal motion:

$$\begin{array}{l}{\displaystyle \frac{d^2{u}_p}{d{t}^2}}+2\mathrm{A}{\displaystyle \frac{{du}_p}{dt}}+\left(A^2+B^2\right)u_p-\left(A^2+B^2\right)U_g\left(y_p\right)\\ + \mathrm{B}\mathrm{g}\left(1-\frac{1}{{\epsilon }_p}\right)=B{v}_0\sin \left(\mathsf{\omega }\mathrm{t}\right)\end{array}$$

(24)

The solution of the resulting equation is given as the sum of the solution corresponding to the damped natural oscillations (left-hand side) and a particular solution accounting for the influence of the forced velocity pulsations in the fluid surrounding the particle (right-hand side). The homogeneous solution of Equation (20) is given by:

$$u_p^{\left(h\right)}\left(\mathrm{t}\right)=U_g\left(y_p\right)+\left(C_3\cos \left(Bt\right)+C_4\sin \left(Bt\right)\right)e^{-At}-{\displaystyle \frac{Bg}{A^2+B^2}}\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)$$

(25)

The homogeneous solution of (21) is:

$${\eta }_2=H_{2}cos\left(\omega t+{\delta }_2\right)$$

(26)

where:

-

η₂ - denotes the particular (forced) oscillatory component of the particle displacement associated with the imposed harmonic velocity pulsations.

The constants δ₂ and H₂ are determined following the procedure described above. By equating and summing the coefficients of the cos(ωt) terms in the expression for δ₂, the particular solution is obtained as:

$${\delta }_2={\displaystyle \frac{A^2+B^2-{\omega }^2}{2A\omega }}$$

(27)

$$H_2={\displaystyle \frac{B{v}_0\omega }{\left(A^2-B^2-{\omega }^2\right)sin{\delta }_2+2A\omega cos{\delta }_2}}$$

(28)

The complete solution for longitudinal velocity is:

$$u_p=\left(C_{3}cos\left(Bt\right)+C_{4}sin\left(Bt\right)\right)e^{-At}-{\displaystyle \frac{Bg}{A^2+B^2}}\left(1-{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\epsilon }_p$}\right.}\right)+H_{2}cos\left(\omega t+{\delta }_2\right).$$

(29)

The integration constants are determined from the initial conditions:

$$C_3={\displaystyle \frac{Bg}{A^2+B^2}}\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)-H_{2}cos\left({\delta }_2\right)$$

(30)

$$C_4={\displaystyle \frac{A}{B}}{C}_3+{\displaystyle \frac{H_2\omega }{B}}sin{\delta }_2$$

(31)

The longitudinal displacement is:

This section establishes a consistent physical and mathematical model for the motion of a single solid particle in a viscous boundary-layer flow subjected to harmonic transverse velocity pulsations. Starting from the governing force balance, coupled differential equations for the transverse and longitudinal velocity components are derived and solved analytically. The resulting closed-form solutions capture both transient and steady-state particle dynamics, accounting for viscous damping, Saffman-induced coupling, gravitational effects, density ratio influence, and externally imposed oscillatory forcing, and provide a solid basis for further parametric and numerical analysis.

3. Advanced Analytical Extensions of the Governing Particle Motion Model

An extended analytical framework is developed to describe the motion of individual solid particles within a horizontal boundary layer subjected to transverse velocity fluctuations. Building on the linear particle dynamics introduced earlier, the formulation incorporates impulsive contributions to represent sudden, localized changes in particle velocity arising from micro-vortical structures, particle–wall interactions, or transient turbulent accelerations.

The governing equations are expressed in matrix–vector form, enabling spectral analysis and the application of Laplace transforms to capture both transient and steady-state responses. This approach facilitates the investigation of resonance phenomena, stability via eigenvalue analysis, and phase-space dynamics. By combining smooth harmonic fluctuations with discrete impulsive effects, the framework provides a rigorous link between idealized particle motion and the complex behaviour observed in boundary-layer flows.

3.1. Impulsive Differential Extension of the Governing Dynamics

Real boundary layers are inherently non-ideal and rarely exhibit perfectly linear behaviour. In low-turbulence or transitional flow regimes, the particle may experience sudden, localized deviations in the surrounding velocity field. Such phenomena include micro-vortical structures with extremely short timescales, particle–wall impact interactions, impulsive turbulent cores, as well as mechanical or vibrational perturbations from the environment. These events are inherently impulsive, inducing nearly discrete jumps in particle velocity that cannot be faithfully represented by smooth harmonic functions such as v₀cos(ωt).

To rigorously account for these abrupt and localized velocity fluctuations, the baseline system is extended through the inclusion of Dirac delta functions, δ(t−t_k), which represent discrete impulsive events occurring at specific instants t_k. The extended system can be expressed in vectorial form as:

$$\begin{array}{l}{\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{u}_p\\ v_p\end{array}\right]=\left[\begin{array}{c}A\left(U_g-u_p\right)-B\left(v_g^{\text{'}}-v_p\right)\\ A\left(v_g^{\text{'}}-v_p\right)-B\left(U_g-u_p\right)-g\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)\end{array}\right]\\ +{\displaystyle \sum _{k=1}^N}\left[\begin{array}{c}{I}_{x,k}\\ I_{y,k}\end{array}\right]\mathsf{\delta }(\mathrm{t}-t_k)\end{array}$$

(32)

where:

-

I_k=[_Ix,k,I_y,k]^T - denotes the impulsive contributions in the longitudinal and transverse directions.

Figure 2 provides a conceptual view of the impulsive forcing term introduced in the governing equations, illustrating how discrete events produce instantaneous changes in the particle velocity components superimposed on the continuous motion.

Integration of the system in the immediate vicinity of t_k yields the corresponding velocity jump condition:

$$u_p\left(t_k^{+}\right)-u_p\left(t_k^{-}\right)=I_k,$$

(33)

demonstrating that particle velocity undergoes a discrete discontinuity, whereas acceleration exhibits a singular component.

The incorporation of impulsive terms offers a scientifically robust framework for capturing short-duration acceleration phenomena, including deposition, rebound, or detachment events at the wall, while preserving the linear character of the underlying system. Light particles (τp≪1) respond almost instantaneously to impulsive forcing, whereas heavy particles (τp≫1) display attenuated responses dominated by inertial effects. This extended formulation therefore provides a rigorous mathematical tool for the analysis of realistic boundary layer flows, seamlessly integrating smooth harmonic fluctuations with discrete, impulsive events that typify turbulent or impact-driven processes.

3.2. Matrix Formulation of the Particle Motion System

The governing equations describing the motion of a particle in a horizontal boundary layer can be expressed in compact matrix–vector form, which consolidates the longitudinal and transverse velocity components into a unified dynamical representation. This is not merely a matter of notation: the two components are intrinsically coupled through viscous relaxation and lift effects and therefore should not be treated as independent scalar equations. The matrix formulation further enables continuous forcing and impulsive events to be incorporated within the same linear framework, thereby preserving analytical clarity and physical interpretability.

Introducing the state vector:

$$X\left[t\right]=\left[\begin{array}{c}{u}_p\left(t\right)\\ v_p \left(t\right)\end{array}\right],$$

the particle–motion system can be written as:

$${\displaystyle \frac{dX\left[t\right]}{dt}}=\mathrm{M}\mathrm{X}\left(\mathrm{t}\right)+F\left(t,y_p\right)+\sum _{k=1}^N{I}_k\mathsf{\delta }\left(\mathrm{t}-t_k\right)$$

(34)

where:

$$\begin{array}{l}M= \left[\begin{array}{cc}-A& B\\ B& -A\end{array}\right],\\ F\left(t\right)=\left[\begin{array}{c}A{U}_g\left(y_p\right)-Bv\text{'}\left(t\right)\\ -B{U}_g\left(y_p\right)+Av\text{'}\left(t\right)-g\left(1-{\displaystyle \frac{1}{{\epsilon }_p}}\right)\end{array}\right],\\ I_k=\left[\begin{array}{c}{I}_{x,k}\\ I_{y,k}\end{array}\right]\end{array}$$

The continuous forcing term F(t) encompasses the longitudinal carrier flow, the imposed transverse pulsations and the gravitational contribution, whereas the impulsive term I_k δ(t−t_k) represents short-duration events such as particle–wall interactions or micro-vortical impacts. Integrating the system across an impulsive instant t=t_k yields the jump condition

$$X\left(t_k^{+}\right)-X\left(t_k^{-}\right)=I_k,$$

(35)

which provides a mathematically consistent description of instantaneous velocity changes without altering the linear operator governing the system.

The impulsive contribution preserves linearity because it is introduced as an additive forcing term rather than as a modification of M. Consequently, the intrinsic dynamical properties—stability, characteristic time scales and natural frequencies—are dictated by the spectrum of M, while the impulses affect only the particular solution. In many situations involving steady or purely periodic excitation, the impulsive part may be safely omitted without changing the qualitative conclusions regarding stability or phase behaviour. Conversely, when transient phenomena or short-lived interactions are of interest, the impulsive representation provides the necessary analytical flexibility while remaining within a strictly linear formulation.

The matrix M encapsulates the essential physical mechanisms: the parameter A corresponds to viscous relaxation toward the carrier-phase velocity, and the parameter B introduces transverse coupling associated with lift effects. Light particles with small relaxation time (τ_p≪1) follow the carrier flow almost instantaneously, whereas heavy particles (τ_p≫1) exhibit inertia-dominated dynamics with pronounced phase lag.

Overall, the matrix formulation offers a concise, analytically transparent and physically interpretable description of the particle dynamics and forms a natural basis for the subsequent examination of transient, forced and impulsive regimes in boundary-layer flows, which is developed in the next section.

3.3. Analytical Solutions, Phase–Space Representation and Asymptotic Behaviour

The analytical strategy adopted in this section is motivated by the need to clearly distinguish, within a unified and physically consistent framework, the transient relaxation, the forced oscillatory response, and the long-term asymptotic behaviour of a single solid particle moving in a boundary layer subjected to transverse velocity pulsations. Since the governing equations have already been formulated as a linear time-invariant matrix system, the analysis that follows does not introduce additional modelling assumptions but instead applies standard operator-based tools of linear system theory directly to the established physical model.

The coupled longitudinal and transverse particle dynamics are written compactly in state-space form as:

$${\displaystyle \frac{d}{dt}}X\left[t\right]=\mathrm{M}\mathrm{X}\left(\mathrm{t}\right)+F\left(t\right),X\left[t\right]=\left[\begin{array}{c}{u}_p\left(t\right)\\ v_p \left(t\right)\end{array}\right]$$

where the system matrix M encapsulates the internal dynamics associated with viscous relaxation and lift-induced coupling between velocity components, while the forcing vector F(t) accounts for the mean carrier flow, the imposed harmonic transverse pulsations, and gravitational effects through the density ratio ε_p.

To separate transient effects associated with the initial conditions from the response induced by external forcing, the unilateral Laplace transform is applied to the full matrix system. The Laplace-transformed state vector $\tilde{X}\left(s\right)$and forcing vector $\tilde{F}\left(s\right)$ are defined in the standard manner, while the time derivative transforms according to:

$${\displaystyle \frac{dX}{dt}}\to s\tilde{X}\left(s\right)-X\left(0\right)$$

(36)

This yields the algebraic equation:

$$s\tilde{X}\left(s\right)-X\left(0\right)=M\tilde{X}\left(s\right)+\tilde{F}\left(s\right),$$

(37)

which can be rearranged to give the general solution in the Laplace domain:

$$\tilde{X}\left(s\right)={\left(sI-M\right)}^{-1}X\left(0\right)+{\left(sI-M\right)}^{-1}\tilde{F}\left(s\right)$$

(38)

Here, s is the complex Laplace variable, III is the identity matrix, and (sI−M⁾⁻¹ denotes the resolvent operator of the linear dynamical system. This operator plays the role of a transfer operator that maps both the initial condition and the external forcing to the system response and therefore contains the complete dynamical information of the particle–fluid interaction.

The spectral properties of the resolvent operator are governed by the eigenvalues of the system matrix M, which satisfy:

$$\mathit{det}(M-\lambda I)={(\lambda +A)}^2+B^2=0 {\lambda }_{\mathrm{1,2}}=-A\pm iB$$

(39)

These eigenvalues correspond to exponentially decaying oscillatory modes with decay rate A and intrinsic angular frequency B. In the phase plane (u_p,v_p), solutions trace spiral trajectories approaching a stationary point, with narrow rapid spirals for light particles (τ_p≪1) and elongated inertial trajectories for heavy particles (τ_p≫1).

The forcing vector can be naturally decomposed into a stationary and an oscillatory component,

$$F\left(t\right)=F_{stat}+F_{osc}\left(t\right)$$

(40)

where F_stat contains the time-averaged contributions of the base flow and gravity, and F_osc(t) arises exclusively from the imposed transverse velocity pulsations. In the long-time limit, the transient contribution associated with the initial condition vanishes for A>0, while the oscillatory forcing produces a bounded periodic response. The asymptotic stationary state is therefore determined by the balance between the internal dynamics and the mean external forcing and is given by:

$$X_{\infty }=M^{-1}{F}_{stat}$$

(41)

Physically, X∞ represents the mean drift state of the particle in the boundary layer, around which periodic oscillations are superimposed when harmonic forcing is present.

Because the governing system is linear and time-invariant, a harmonic excitation of the form v_g′(t)=v₀cos(ωt), after decay of transients, to a steady periodic particle response at the same frequency. The oscillatory component of the solution can therefore be expressed in complex-amplitude form as

$$X\left(t\right)=\mathfrak{R}\left\{\widehat{X}\left(\omega \right)e^{i\omega t}\right\}$$

(42)

enables an explicit evaluation of amplitude and phase relative to the forcing frequency ω. Resonance emerges when ω≈B, demonstrating that the imposed oscillations can significantly modify the transverse motion of the particle under certain parameter conditions.

The complete time-domain solution can be written as a superposition of homogeneous and forced contributions:

$$\mathrm{X}\left(t\right)=e^{Mt}X\left(0\right)+{\int }_0^t{e}^{M\left(t-\tau \right)}F\left(\tau \right)d\left(\tau \right),$$

(43)

where the first term governs the transient relaxation and the integral term defines the forced response. This representation makes explicit how the observed particle motion arises from the interaction between the intrinsic modal structure of the system and the applied external forcing.

The logical structure of the operator-based analysis and the connection between the Laplace-domain formulation, spectral properties, time-domain response, and asymptotic behaviour are summarized conceptually in Figure 3.

The figure serves as an explanatory aid, highlighting the central role of the resolvent operator (sI−M⁾⁻¹ in governing transient decay, resonant response, and long-term particle dynamics, without representing numerical simulations or previously published visualizations.

3.4. Dynamic Formulation, Modal Spectrum, and Lyapunov Stability

The dynamic formulation of the particle-motion system defines a time-dissipative linear model in which the temporal evolution of the state vector X(t) is governed by the spectral properties of the operator M and the external forcing F(t). This representation allows a comprehensive description of transient and steady regimes and enables the assessment of dynamic stability resulting from the balance between relaxation, inertia, and imposed oscillatory effects. The system therefore constitutes a genuine dynamic model rather than a static configuration, with a full temporal evolution and a well-defined internal modal structure.

The spectral representation of the system matrix is obtained through:

$$M=V\Lambda V^{-1}, \Lambda =diag\left({\lambda }_1,{\lambda }_2\right), {\lambda }_{\mathrm{1,2}}=-A\pm iB$$

(44)

where:

-

V=[v₁,v₂] is the matrix of right eigenvectors of M.

The system thus exhibits damped oscillations whose decay rate is controlled by A, while the modal frequency is given by B. For a periodic excitation, e.g., v_g′(t)=v₀cos(ωt), a complex representation yields:

$$\left(i\omega I-M\right)\tilde{X}=\tilde{F}, \mathrm{H}\left(\mathsf{\omega }\right)={\left(i\omega I-M\right)}^{-1}$$

(45)

where:

-

H(ω) denotes the dynamic transfer operator determining amplitude, phase, and resonant amplification (ω≃B).

Dynamic stability can be assessed through a quadratic Lyapunov function:

$$V\left(x\right)={\displaystyle \frac{1}{2}}{‖X‖}^2, {\displaystyle \frac{dV}{dt}}=X^{T}MX<0 (A>0)$$

(46)

which guarantees exponential decay of the homogeneous component and the existence of an attracting stationary state X∞=M⁻¹F_stat. Accordingly, the system belongs to the class of exponentially stable linear dissipative dynamic systems, where the response amplitude remains bounded and fully controllable through the parameters A,B,τ_p and the forcing frequency ω. The spectral properties thus provide a rigorous framework for analysing resonance, forced response, and long-term behaviour under realistic boundary-layer conditions.

4. Numerical Analysis and Physical Interpretation

4.1. Baseline Analytical and Numerical Framework

The analytical formulation developed in Section 3 provides a rigorous mathematical basis for analysing the coupled transverse and longitudinal dynamics of finite-size solid particles in a horizontal laminar boundary layer subjected to harmonic cross-flow velocity pulsations. By expressing the governing equations in compact matrix form and applying the Laplace transform, the particle motion is decomposed into transient, oscillatory, and asymptotic components, allowing a clear physical interpretation of the mechanisms governing particle transport under continuous forcing.

In order to illustrate the analytical solutions in a controlled and transparent manner, a baseline numerical configuration is introduced. Water at standard laboratory conditions is adopted as the carrier fluid, and the physical properties of the fluid together with the characteristic parameters of the boundary layer are summarised in Table 1. These parameters define the shear scale and flow conditions responsible for the generation of lift and viscous relaxation effects acting on the particle.

The properties of the dispersed phase are given in Table 2, where two representative cases are considered: a light particle characterised by a short relaxation time and a heavy particle exhibiting inertia-dominated behaviour.

The particles differ only in density and relaxation time, while identical initial conditions are prescribed for both cases. This ensures that any observed differences in particle motion arise solely from inertial effects rather than from differences in the initial kinematic state.

The numerical integration of the governing equations is performed using the parameters listed in Table 3, including the time step, total simulation duration, and numerical integration scheme. These values are selected to accurately resolve both the transient decay and the established oscillatory regime induced by the harmonic excitation.

At each time step, the longitudinal carrier-phase velocity is evaluated at the instantaneous particle position as U_g(y_p(t)), ensuring consistent coupling between transverse migration and longitudinal acceleration.

Within this first numerical stage, the impulsive forcing term retained in the general analytical formulation is set to zero, such that the results correspond to continuous harmonic forcing. The specific numerical datasets used to construct the graphical results presented below are provided as Supplementary Tables S1–S4, ensuring full reproducibility of the reported simulations.

4.2. Transverse Particle Dynamics

The transverse motion constitutes the primary mechanism responsible for segregation between light and heavy particles in a horizontal boundary layer subjected to cross-flow velocity pulsations. The analytical expression for the transverse particle velocity v_p(t), derived in Section 3, is evaluated numerically using the baseline parameter set defined in Table 1, Table 2, Table 3 and Table 4. The resulting time histories of the transverse velocity are shown in Figure 4 (Supplementary Table S1).

Two qualitatively distinct dynamical regimes are clearly observed. For the light particle, the transient response decays rapidly and the velocity oscillates about a positive mean value, indicating a sustained upward drift toward regions of higher carrier-fluid velocity. In contrast, the heavy particle exhibits an inertia-dominated response characterised by a negative mean transverse velocity, corresponding to gradual migration toward slower near-wall sublayers.

In both cases, the oscillatory behaviour reflects the imposed harmonic excitation, while the decay rate is governed by the viscous relaxation coefficient and the particle relaxation time.

The results presented in the figures correspond to physically calibrated responses; the applied rescaling accounts for finite-size effects and uncertainties in local shear estimation and does not introduce additional physical mechanisms.

The transverse displacement y(t), obtained by time integration of v_p(t), is presented in Figure 5 (Supplementary Table S2). The light particle undergoes continuous upward migration, whereas the heavy particle experiences a pronounced downward drift dominated by the effective gravitational contribution. Although oscillatory modulations are superimposed on the trajectories, the long-term trend is independent of the initial phase and is instead controlled by particle inertia, relaxation time, and the relationship between the system’s natural frequency and the excitation frequency.

These results demonstrate that even in the absence of turbulence or spatial shear gradients, harmonic transverse velocity pulsations alone are sufficient to induce vertical particle segregation in laminar boundary-layer flows.

4.3. Longitudinal Particle Dynamics

The longitudinal motion of the particles arises as a direct consequence of their transverse migration, since different regions of the boundary layer are associated with different values of the carrier-phase velocity. The analytical solution for the longitudinal particle velocity u_p(t), derived in Section 3, is evaluated numerically using the same baseline configuration. The corresponding time histories are shown in Figure 6 (Supplementary Table S3).

For the light particle, the longitudinal velocity gradually increases relative to its initial value, reflecting its migration into faster regions of the boundary layer. Conversely, the heavy particle exhibits a progressive decrease in longitudinal velocity as a result of its downward drift toward slower sublayers. Periodic oscillations are observed during the initial stage for both particles due to the transverse harmonic excitation; however, their amplitude and decay rate depend strongly on the particle relaxation properties.

The longitudinal displacement x(t), obtained by integrating u_p(t), is shown in Figure 7 (Supplementary Table S4). The light particle advances rapidly and moves significantly ahead of the heavy particle, while the latter exhibits a clear lag associated with residence in regions of lower transport velocity. This divergence intensifies with time, despite identical initial longitudinal velocities for both particles.

From a physical standpoint, the longitudinal dynamics cannot be interpreted independently of the transverse behaviour. The transverse migration determines access to faster or slower regions of the boundary layer and thus governs the longitudinal transport of the particles. This sections therefore complete the baseline numerical illustration of the analytical model under continuous harmonic forcing, providing a reference framework for the extended numerical examples considered in subsequent stages.

4.4. Impulsive Particle Dynamics

The impulsive extension of the particle–motion model introduced in Section 3.1 is illustrated numerically in order to quantify the response of light and heavy particles to short–duration perturbations. In the present simulations, the impulsive component is incorporated exactly as defined in Section 3.1, while all remaining physical and numerical parameters coincide with the baseline configuration used in the previous sections. The impulsive scenario represents intermittent disturbances commonly encountered in boundary-layer flows, such as micro-vortical structures, particle–wall interactions, or short turbulent bursts, which cannot be captured by smooth harmonic forcing alone. A finite sequence of impulsive velocity jumps is prescribed at selected instants, with identical impulse amplitudes applied to both light and heavy particles; the impulse times and the corresponding longitudinal and transverse velocity jumps are reported in Supplementary Table S5.

Figure 8 illustrates the phase–space response of light and heavy particles subjected to the impulsive forcing scenario. The left part of the figure highlights the short-time behaviour dominated by discrete velocity jumps induced by impulsive events, whereas the right part shows the long-time evolution revealing the spiral attractor of the continuous matrix dynamics. Immediately after each impulse, the particle velocity exhibits a discontinuous change, followed by a smooth relaxation governed by viscous damping and lift-induced coupling.

A pronounced difference between particle types is observed. Light particles rapidly relax toward a compact spiral attractor due to strong viscous relaxation, while heavy particles display wider phase-space excursions and a significantly slower return toward the attractor, reflecting the persistence of inertial effects. Despite these transient differences, the global attractor structure remains unchanged, indicating that impulsive forcing primarily produces temporary displacements rather than altering the intrinsic stability properties of the system.

These numerical results confirm the analytical predictions of Section 3.1 and demonstrate that the impulsive differential extension provides a physically meaningful and analytically consistent description of particle dynamics in boundary layers subjected to intermittent disturbances.

4.5. Matrix Dynamics and Transient Relaxation

In this section, the matrix formulation of the particle–motion system introduced in Section 3.2 is examined through numerical simulations in order to characterise the intrinsic continuous dynamics of the coupled longitudinal and transverse particle velocities. The analysis focuses on the transient relaxation behaviour in the absence of impulsive disturbances, so that the observed particle motion is governed solely by the linear time-invariant system matrix and its spectral properties.

The governing dynamics are defined by a two-dimensional matrix system whose eigenvalues are given by:

$${\lambda }_{\mathrm{1,2}}=-A\pm iB$$

This eigenstructure corresponds to exponentially damped oscillatory motion and guarantees the existence of a stable spiral attractor in the phase space (u_p,v_p) for A>0. Numerical integration of the matrix system is performed using identical initial conditions for light and heavy particles, while the model and numerical parameters employed in the simulations are summarised in Supplementary Table S7.

Figure 9 presents the phase–space trajectories obtained from the matrix formulation under continuous harmonic excitation. The left part of the figure illustrates the short-time response immediately following the initial condition, where transient relaxation effects dominate, whereas the right part shows the long-time evolution, revealing the full structure of the spiral attractor associated with the matrix dynamics. In both cases, the trajectories remain bounded and converge asymptotically toward the stationary solution of the system.

A clear distinction between the dynamical behaviour of light and heavy particles is observed. Light particles exhibit rapid contraction of their phase–space trajectories and quickly approach a compact spiral attractor, reflecting strong viscous relaxation and weak inertial memory. In contrast, heavy particles display wider phase–space excursions and a markedly slower convergence toward the attractor, indicating inertia-dominated dynamics and extended transient behaviour.

Despite these differences in relaxation time scales and trajectory geometry, the qualitative structure of the phase–space attractor remains the same for both particle types. This demonstrates that the stability of the motion and the existence of a spiral attractor are intrinsic properties of the matrix formulation itself and are not induced by impulsive or stochastic perturbations. The numerical results therefore provide direct confirmation of the analytical predictions derived in Section 3.2 and support the Lyapunov stability analysis discussed in Section 3.4.

The time-resolved particle velocity data used to construct the phase–space trajectories shown in Figure 9 are provided in Supplementary Table S7, ensuring full reproducibility of the numerical results.

4.6. Resonance Analysis Under Harmonic Forcing

The effect of harmonic transverse velocity pulsations on particle dynamics is examined by analysing the stationary forced response of the linear system derived in Section 3.3 and Section 3.4. In this subsection, attention is restricted to the long-time oscillatory regime, assuming that transient contributions associated with the homogeneous solution have fully decayed due to viscous dissipation.

As established in Section 3.3, under harmonic excitation the state vector of the particle motion admits a representation of the form(eq.42, section 3.3).

$$X\left(t\right)=\mathfrak{R}\left\{\widehat{X}\left(\omega \right)e^{i\omega t}\right\},$$

where:

- X(t)=[u_p(t), v_p(t)]^t and $\widehat{X}$ (ω) denotes the complex amplitude vector associated with the forcing frequency ω.

In the present analysis, attention is restricted to the transverse oscillatory response. The scalar quantity H(ω) denotes the amplitude of the transverse component of the steady harmonic particle response. H(ω)=amplitude of v_p(t) in the steady periodic regime.

Consequently, the frequency dependence of H(ω) is entirely governed by the system coefficients and by the spectral properties of the system matrix discussed in Section 3.4, and no additional modelling assumptions are introduced at this stage.

The resonance behaviour is controlled by the intrinsic oscillatory dynamics associated with the complex conjugate eigenvalues of the system matrix - λ_1,2=−A±iB which define a characteristic frequency scale of the order:

$$\omega \approx \sqrt{A^2{+B}^2}$$

(47)

When the forcing frequency approaches this intrinsic scale, the magnitude H(ω) increases significantly, indicating resonance-like amplification of transverse particle oscillations. At low forcing frequencies, the particle motion remains nearly in phase with the imposed oscillations, while near resonance a pronounced phase delay develops, reflecting the inertia-controlled response of the particle.

The influence of particle inertia is encapsulated through the relaxation time τ_p, which enters the coefficients of the governing equations and controls the frequency dependence of the steady oscillatory response.

Particles with weak inertial effects exhibit pronounced amplification over a broad frequency range, reflecting their ability to closely follow the imposed transverse fluctuations of the carrier flow. In contrast, particles with strong inertial resistance display a markedly attenuated response across the entire frequency spectrum, as inertial damping suppresses their capacity to respond to rapid transverse forcing.

From a physical perspective, harmonic transverse pulsations therefore act as a selective segregation mechanism within the boundary layer. Light particles experience enhanced transverse excursions when the forcing frequency approaches the intrinsic system frequency, promoting migration toward regions of higher carrier-phase velocity. Heavy particles, on the other hand, remain largely confined to near-wall regions due to inertial damping and reduced coupling to the oscillatory forcing. This behaviour provides a direct numerical manifestation of the modal properties and Lyapunov stability characteristics established in Section 3.4.

The frequency-dependent steady-state amplitudes H(ω) used to construct the resonance curves shown in Figure 10 are obtained numerically from time-domain simulations and reported in Supplementary Table S8.

4.7. Multi-Harmonic Transverse Forcing

In practical boundary-layer flows, transverse velocity fluctuations are rarely characterised by a single harmonic component. Instead, the excitation typically consists of a superposition of oscillatory contributions with different amplitudes, frequencies and phases. To examine the particle response under such conditions, the transverse carrier-phase pulsation is generalised to a multi-harmonic form

$$v_g^{\text{'}}\left(t\right)=\sum _{n=1}^N{v}_n\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{n}t+{\varphi }_n)$$

(48)

where the index nnn labels the individual harmonic components included in the excitation.

Each harmonic contribution is fully characterised by its amplitude v_n, angular frequency ω_n, and phase ϕ_n. The amplitude v_n is defined as:

$$v_n=max\left|v_{g,n}^{\text{'}}\left(t\right)\right|$$

(49)

where:

$$v_{g,n}^{\text{'}}\left(t\right)=v_n\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{n}t+{\varphi }_n)$$

(50)

denotes the n-th harmonic component of the transverse carrier-phase velocity. The angular frequency is given by:

$${\omega }_n=2\pi f_n$$

(51)

with f_n denoting the frequency of the n-th harmonic component, while ω_n has the dimensionality rad s⁻¹. The phase shift is expressed as:

$${\varphi }_n={\omega }_n{t}_{n,0}$$

(52)

where t_n,0 represents the phase offset of the n-th harmonic component relative to the reference time t=0.

Owing to the linearity of the governing matrix system established in Section 3.2, the particle response to multi-harmonic forcing can be interpreted as a superposition of the responses to the individual harmonic components. After the decay of transient effects, the transverse particle velocity exhibits a steady oscillatory motion composed of multiple periodic contributions, each associated with one of the imposed forcing frequencies. The resulting time signal is therefore characterised by amplitude modulation, envelope formation and beating phenomena, which are typical of weakly unsteady or transitional boundary-layer environments.

This behaviour is illustrated in the phase–space representation shown in Figure 11, where the superposition of multiple forcing frequencies leads to deformed trajectories in the (u_p,v_p) plane. Numerical simulations based on the representative multi-harmonic parameter sets listed in Supplementary Table S9 demonstrate a pronounced dependence of the particle response on inertia. Light particles effectively follow each harmonic component and exhibit strong amplitude modulation and beating, whereas heavy particles display a markedly attenuated response, with higher-frequency components being strongly damped by inertial effects.

From a physical perspective, multi-harmonic transverse forcing enhances the selectivity of particle segregation within the boundary layer. The simultaneous presence of several excitation frequencies increases the likelihood of partial resonance with the intrinsic oscillatory scale of the system, particularly for light particles. Heavy particles remain comparatively insensitive to such excitation, reinforcing the separation between inertial regimes already observed under purely harmonic forcing. Consequently, multi-harmonic excitation provides a more realistic and physically representative framework for analysing particle redistribution mechanisms in boundary-layer flows.

4.8. Asymptotic Long-Term Behaviour

The analytical framework developed in Section 3.3 and Section 3.4 predicts that, at sufficiently long time scales, the particle dynamics approach an asymptotic regime governed by the spectral properties of the linear system matrix and the persistent external forcing. In this regime, transient contributions associated with the homogeneous solution decay exponentially, and the particle motion is dominated by the forced response induced by transverse velocity pulsations.

For the transverse motion, the exponential decay of the natural modes associated with the eigenvalues of the system matrix (see Section 3.4) ensures that all initial-condition effects vanish as t→∞. The remaining motion is therefore fully determined by the forcing characteristics and the particle relaxation properties. Under harmonic or multi-harmonic excitation, the transverse velocity converges to a bounded oscillatory state, whose amplitude and phase are controlled by particle inertia and by the distribution of forcing frequencies.

The longitudinal dynamics exhibit a qualitatively different asymptotic behaviour. While oscillatory components persist due to coupling with the transverse motion, the longitudinal displacement displays a dominant drift component arising from sustained particle migration across the boundary layer. As a result, even after the transverse oscillations reach a steady periodic regime, the longitudinal position continues to evolve in time, leading to cumulative separation between particles with different inertial properties.

A clear contrast emerges between light and heavy particles. Light particles maintain relatively large transverse oscillation amplitudes and oscillate around an upward-shifted mean position, reflecting their strong coupling to the imposed transverse pulsations and their ability to closely follow rapid fluctuations of the carrier flow. In contrast, heavy particles exhibit a markedly attenuated transverse response due to inertial damping and converge toward a downward-shifted mean position governed primarily by the effective gravitational contribution.

This asymptotic behaviour is illustrated in Figure 12 (Supplementary Table S10), where phase–space trajectories extracted from a restricted steady-state window reveal non-closed, asymmetric geometries for both particle types. Although the detailed trajectory shapes differ, both cases converge toward stable attractor structures determined by the dissipative nature of the system. The persistence of bounded motion confirms the Lyapunov stability of the particle dynamics established in Section 3.4, while the separation of mean positions highlights the long-term segregation mechanism induced by particle inertia.

From a physical perspective, the asymptotic regime represents the ultimate outcome of the combined action of viscous relaxation, lift-induced coupling, gravity, and external transverse forcing. Once transient effects have decayed, the system retains a memory of particle inertia through distinct mean positions and oscillation amplitudes, leading to sustained spatial segregation even in the absence of turbulence. This behaviour underscores the importance of long-time analysis when assessing particle transport and redistribution in boundary-layer flows subjected to oscillatory forcing.

5. Discussion

The present study develops an analytical framework for describing the coupled transverse and longitudinal dynamics of inertial particles in a horizontal boundary-layer flow subjected to transverse velocity pulsations. By formulating the governing equations in a linear matrix form, the analysis establishes a direct connection between transient relaxation, resonance phenomena, and asymptotic long-term behaviour through the spectral properties of the system.

Unlike many previous investigations that rely predominantly on numerical simulations or turbulence-based transport models, the present results demonstrate that deterministic oscillatory forcing alone can induce sustained particle migration and segregation even in laminar or weakly perturbed flows. The analytical solutions reveal how viscous relaxation, lift-induced coupling, and gravity interact to shape particle trajectories, without invoking stochastic turbulence effects.

The resonance behaviour identified in the transverse particle response provides a clear physical interpretation of inertia-dependent sensitivity to oscillatory forcing. Particles with weak inertial resistance exhibit enhanced oscillatory amplitudes over a broad frequency range, whereas inertia-dominated particles remain strongly damped. This distinction is consistent with qualitative observations reported in earlier studies, while the present formulation offers a rigorous mathematical explanation based on modal dynamics and Lyapunov stability.

The asymptotic analysis further shows that, after the decay of transient contributions, particle motion converges toward stable attractor structures whose geometry and mean position depend on particle inertia. Although the transverse oscillations remain bounded, the associated longitudinal drift leads to cumulative separation between particles with different inertial properties. This mechanism highlights a purely deterministic pathway for long-term particle segregation in boundary-layer flows.

From a broader perspective, the results indicate that oscillatory forcing mechanisms—arising from flow pulsations, vibrations, or external actuation—may significantly influence particle transport and redistribution even in the absence of turbulence. The compact analytical nature of the proposed framework makes it particularly suitable for parametric studies and preliminary system design, complementing numerical simulations and experimental investigations.

6. Conclusions

A linear matrix-based analytical framework is developed to describe inertial particle motion in a horizontal boundary-layer flow subjected to transverse velocity pulsations. The resulting closed-form solutions capture transient relaxation, resonance phenomena, and asymptotic long-term behaviour in a unified manner.

The analysis demonstrates that oscillatory transverse forcing can induce pronounced inertia-dependent particle dynamics even in laminar flows. Particles with weak inertial resistance exhibit strong coupling to the imposed pulsations, leading to amplified transverse oscillations and upward-shifted mean positions, whereas inertia-dominated particles experience strong damping and converge toward downward-shifted states governed primarily by gravity.

Resonance phenomena emerge naturally from the spectral properties of the governing system and are shown to enhance transverse particle motion over specific frequency ranges. At long time scales, the exponential decay of natural modes ensures convergence toward stable attractor structures, confirming the Lyapunov stability of the particle dynamics while enabling sustained longitudinal separation between particles with different inertial properties.

Overall, the results provide a clear physical and mathematical interpretation of particle segregation mechanisms driven by oscillatory forcing in boundary layers. The proposed analytical framework offers a compact and predictive tool for studying inertial particle transport and may serve as a foundation for future investigations of controlled particle manipulation and separation in oscillatory flow environments.