In this work, an open-source solver Gerris [27, 28], based on octree and multi-material interface capturing technique, is used to solve the incompressible unsteady Navier-Stokes Equations 1, 2. The solution process includes spatial discretization and time discretization. ∂ ρ u ∂ t + ∇ ⋅ ρ u u = − ∇ p + ∇ ⋅ D + σ κ δ s n , (1) ∂ ρ ∂ t + ∇ ⋅ ρ u = 0 , (2) where u = ( u , v , w ) is the fluid velocity, ρ = ρ ( x , t ) is the fluid density, μ = μ ( x , t ) is the dynamic viscosity, p is the pressure, D is the viscous stress tensor, defined by Equation 3. δ s is the Dirac distribution function (which means the surface tension is concentrated on the liquid-gas surface), σ is the surface tension coefficient, κ is the mean interface curvature, and n is the unit normal vector along the interface. D = 2 μ S − 2 3 μ ∇ ⋅ u I , S = 1 2 ∇ u + ∇ u T . (3)

For two-phase flow, the volume fraction α = α ( x , t ) of the first phase is introduced, then the density and viscosity are defined by Equation 4: ρ = α ρ 1 + 1 − α ρ 2 , μ = α μ 1 + 1 − α μ 2 , (4) where ρ 1 and ρ 2 denote the densities of the primary and secondary phases, respectively. μ 1 and μ 2 are the dynamic viscosities of the primary and secondary phases, respectively. The convection equation of volume fraction α ( x , t ) can be expressed as Equation 5: ∂ α ∂ t + ∇ ⋅ α u = 0 . (5) where the phase fraction α is defined as Equation 6: α x , t = 0 p h a s e 1 , 0 < α < 1 i n t e r f a c e , 1 p h a s e 2 (6)

The gas–liquid interface is tracked and reconstructed using the Volume of Fluid (VOF) interface capture method. The volume fraction equation is solved for each time step to determine the distribution of volume fractions in the flow field, and then reconstruct the interface information and obtain the flow flux of the grid cells near the interface. The piecewise-linear geometrical Volume-of-Fluid scheme is used for reconstruction [29]. This method achieves second-order accuracy, does not require special treatment of phenomena such as interface breakup and fusion, and effectively captures sharp interface with high precision.

The convection term and diffusion term are discretized by Bell–Colella–Glaz second-order upwind scheme [30] and Crank–Nichoslon scheme [31], respectively.

Moreover, space is discretized by hierarchical quadtree partitioning (octree in 3D) in Gerris. Since jet atomization is a multi-scale physical process, employing a uniform mesh across the entire calculation domain would incur prohibitively high computational costs. Therefore, the Adaptive Mesh Refinement (AMR) method is adopted in the present study to maintain accuracy while minimizing computational cost. By applying interface and vorticity as criteria, the grid is dynamically refined near the interface at each time step, thus ensuring precise tracking of interface deformation.

The accuracy of the numerical method used in this study has been validated in our previous work on swirling liquid sheet [26]. The breakup process of an infinite liquid jet subjected to a small initial disturbance was simulated. The numerical results are in good agreement with the Rayleigh linearized theoretical solution, confirming the model’s ability to capture the instability and breakup dynamics of jets.

The physical configuration of the base case consists of a liquid diesel fuel stream injected through a circular orifice with diameter D = 100   μ m into a gaseous crossflow environment. The properties of diesel fuel and the gas are summarized in Table 1. The crossflow velocity matches the liquid jet velocity at U g = U l = 30 m/s. As shown in Figure 1, the dimensions of the computational domain is 40 D × 20 D × 10 D . The gaseous flow direction aligns with the positive x axis, while the liquid circular jet propagates along − y . A velocity boundary condition is applied to both the gas inlet (left side) and the liquid inlet (upper side), with an outlet boundary condition applied to the right side and symmetric boundary conditions applied to the remaining three boundaries.

To simplify the model and concentrate on the influence of pulsation at the nozzle exit, both gasoues and liquid flow are considered as uniform flows. If turbulence is taken into account, there might exist a coupling or competitive relationship between the disturbance frequency and the turbulent pulsation, which would complicate the causal relationship in the subsequent analysis. In the research conducted by Shinjo and Umemura [32], the inlet velocity was also set as a uniform distribution, and other nozzle interference factors such as turbulence were ignored to establish a clear causal relationship in the analysis. Therefore, for the same reason, an idealized nozzle configuration is adopted in this study.

Figure 2 shows the mesh resolution employed in the present study. By appling the interface-basded and vorticity-based refinement criteria, the grid is dynamically refined both at the gas-liquid interface and in regions of high vorticity. In the present study, the minimum grid size used in the calculation is 2 μ m [25]. investigated the streamwise forcing effects on the liquid circular jet using a minimum grid size of 2 μ m , closely matching the present mesh resolution. Yang and Turan [33] demonstrated that a grid resolution with a minimum scale of 2 μ m adequately resolves ligament structures, droplets, and the jet breakup process through comprehensive quantitative and qualitative validation. Given the identical grid resolution adopted in our study, the reliability of the present results is thereby confirmed. In present study, the simulated jet structures and mean jet diameter obtained with the 2 μ m minimum grid size demonstrate strong agreement with the existing results in the literature. Additionally, simulations using a coarser grid with a minimum scale of 4 μ m were also conducted, and the results were compared in Figure 3. From the transient gas-liquid interface profiles and jet trajectories, it is evident that although the jet trajectories achieved with these two grids align closely, the coarse grid fails to accurately resolve surface waves and column breakup effectively. Thus, the 2 μ m minimum grid scale adopted in this work better resolves the breakup process.

Building upon the grid independence verification, the authors further examine how domain boundaries affect flow characteristics. To investigate the effect of domain size, we conducted additional simulations with an extended computational domain ( 50 D × 20 D × 10 D ) and compare these results with the baseline case, and the distance between the gas inlet and the center of the circular jet was increased from 5 D to 15 D . Figure 4 shows the comparison between the two simulation cases. Based on the gas-liquid profiles in Figures 4A,B, the liquid core and large-scale structures closely match. In Figure 4B, where the liquid trajectories are compared, the two curves coincide with each other very well. These results demonstrate that the important characteristics of the liquid jet remain unaffected by changing the domain size. Consequently, to avoid excessive computational expense in subsequent simulations, we adopt the smaller domain configuration where the crossflow inlet is positioned 5 D from the jet centerline.

In the present study, the key dimensionless parameters are defined and given as follows: the liquid Weber number W e l = ρ l D U l 2 σ = 2544 ; the gas Weber number W e g = ρ g D U g 2 σ = 103.5 ; the liquid-gas momentum flux ratio q = ρ l U l 2 ρ g U g 2 = 24.58 ; and the Ohnesorge number O h = μ l ( ρ l D σ ) 1 / 2 = 0.0569 . According to the primary breakup regime classification proposed by Wu et al. [2], the operating condition in present study lies within the vicinity of the multi-mode and shear breakup transition. It should be noted that previous literature has proposed various classification maps for primary breakup regimes [3]. Wu et al. [2] utilized W e g − q parameter space for classification, while Mazallon et al. [4] employs W e g − O h as criteria. However, there are currently no universally accepted standard classification criteria, and further research is required to establish consensus.

First, we simulated the undisturbed liquid jet in crossflow, as illustrated in Figure 5. The characteristic breakup process of JICF [2] can be clearly identified. The left column of Figure 5 shows the liquid-gas interface through a contour corresponding to a liquid volume fraction of 0.5, where color mapping indicates axial velocity. Figure 5 reveals that the axial velocity of the liquid column progressively diminishes with increasing distance from the nozzle, while surface breakup phenomena emerge due to crossflow-induced shear effects. Further, the leeward vortex of the liquid jet redirects droplet velocities towards the + y direction. The right column of Figure 5, representing the gas-liquid interface at z = 0 , demonstrates both column breakup and a large number of ligaments and droplets.

The temporal evolution of a jet in crossflow is illustrated in Figure 6, which reveals clear evidence of both surface and column breakup phenomena. Following injection from the nozzle, the jet’s leading edge deforms into a mushroom-like structure due to crossflow impingement. Under sustained aerodynamic loading from the gaseous crossflow, the liquid jet core undergoes progressive flattening and curvature development, exhibiting pronounced bending alignment with the crossflow direction. Figure 7 provides a systematic analysis of cross-sectional deformation dynamics. As the liquid jet propagates downstream, its transverse dimension (along the z -direction) demonstrates progressive expansion, while the cross-sectional geometry gradually transitions from circular to crescent morphology. This geometric transformation enhances windward surface exposure to transverse flow, thereby amplifying aerodynamic drag and accelerating liquid column bending. The interaction between high-speed gas flow shear stresses at lateral edges and vortex-induced coiling mechanisms leads to progressive thinning of the liquid film at the periphery. This thinning process culminates in subsequent shedding of ligaments and droplets. Further downstream, the crescent-shaped liquid core gradually reverts to an irregular geometry through dynamic restructuring. A distinctive transverse edge develops on the deformed liquid core, from which secondary ligament and droplet shedding occurs under sustained gas flow shear–a process bearing mechanistic resemblance to secondary droplet breakup phenomena.

It should be emphasized that under the coupled transverse-axial gas flow effects, bag breakup occurs near the liquid jet tip. To clarify this mechanism, Figure 8 displays the temporal evolution of bag breakup phenomena, with critical stages marked by red arrows. The elevated stagnation pressure on the windward side of the flattened liquid column induces bag formation at its downstream center, where the bag progressively enlarges in the transverse gas flow direction. As shown in Figure 8, when approaching maximum expansion, the liquid bag detaches from the column apex through a detachment process analogous to droplet bag breakup. However, according to breakup regime classifications established in existing literature, the dimensionless numbers in the present study deviate from those characterizing conventional bag breakup. Our simulation results demonstrate that local breakup patterns depend not only on the initial global field but also exhibit significant coupling with local flow conditions. Consequently, exclusive reliance on dimensionless numbers for regime classification is inadequate; this finding is crucial for improving atomization modeling accuracy.

Under the current JICF conditions, we quantify the dominant wavelengths of three competing instability modes–Rayleigh-Taylor (RT), Kelvin-Helmholtz (KH), and capillary instabilities–through theoretical analysis. This wavelength comparison enables definitive identification of the prevailing instability mechanism governing the present JICF breakup process.

The evolution of the liquid jet on the windward side is presented in Figure 10. Consistent with the analysis in the previous section, increasing distance from the nozzle induces progressive destabilization of the jet surface through Rayleigh-Taylor (RT) instability, resulting in the formation of surface wave structures on the windward side. As shown by black circles in Figures 10A–I, airflow-induced stretching causes both sides of the “ Λ ”-shaped waves to gradually break into ligaments from the troughs. These ligaments subsequently contract and ultimately disintegrate into small droplets, as marked by black circles in Figures 10J–L. The bag breakup process in surface wave troughs is highlighted with red circles in Figure 10. Transverse airflow generates a high-velocity zone near the lateral edges of surface wave troughs, forcing liquid migration toward the crossflow direction ( + x ) and forming liquid bags. These bags become progressively thinner and eventually develop a central perforation. Subsequent expansion of these perforations leads to the formation of closed ligaments around the holes, as demonstrated by red circles in Figures 10D–F. Ultimately, these closed ligaments rupture into striped ligaments, shown by red circles in Figures 10G–I.

Ligaments formed through surface wave crest shedding exhibit larger dimensions, consequently producing correspondingly larger droplets. Additionally, as indicated by blue arrows in Figure 10A, intense shear effects at the liquid jet’s spanwise edge - caused by airflow motion along the surface ( z -direction) - induce shear breakup modes. This breakup mechanism generates comparatively smaller droplets.

The aforementioned breakup modes occur periodically near the nozzle. Both bag and shear breakup modes produce relatively small ligaments and droplets, while surface wave crest shedding generates larger ones. Small droplets, characterized by low inertia, experience significant influence from the surrounding flow field, resulting in broad dispersion. Conversely, larger droplets remain predominantly concentrated near the liquid core.

We next analyze and compare the temporal evolution of the disturbed liquid jet under the operating conditions outlined in Table 2.

Figure 11 illustrates the temporal evolution of the liquid jet’s global structure in Case 2. The liquid column breakup demonstrates a periodic flapping behavior, exhibiting a significantly reduced breakup length relative to undisturbed conditions. Velocity pulsation at the nozzle induces periodic contracting and stretching of the liquid core, resulting in wave node formation. Blue arrows in Figure 11 track the temporal development of these wave nodes, where spanwise dilatation ( z -direction) of the inter-node liquid core occurs due to crossflow-induced aerodynamic forcing. The subsequently flattened liquid film undergoes thinning until ultimate breaks. This breakup pattern recurs with each pulsation cycle.

Figure 12 shows the temporal evolution of the jet’s windward surface viewed along the x -direction. Rayleigh-Taylor (RT) instability-induced fluctuations persist in the planar liquid core between wave nodes. Though, the formation of characteristic undisturbed jet features - specifically, the Λ -shaped surface waves - is prevented by the flow field perturbationthe induced by the pulsed jet. Consequently, no residual spine-like structure remains at this excitation frequency, with jet column disintegration occurring between proximal wave nodes near the nozzle, yielding reduced column breakup length.

Figure 13 presents the temporal evolution of the liquid jet’s global structure in Case 3 - where the excitation frequency precisely matches the intrinsic Rayleigh-Taylor (RT) instability frequency - with black lines depicting the gas-liquid interface in the central plane ( z = 0 ) . Contrasting with Figure 12, periodic flapping breakup of the liquid core is completely absent at this excitation frequency. Nozzle velocity pulsation induces periodic surface fluctuations. Through klystron effect, knot-like liquid structures emerge with sustained wave node amplitude over time. This markedly distinct behavior observed here differs from streamwise-forced circular jet dynamics, where the surface wave on the liquid core, induced by a 300 kHz frequency disturbance, gradually diminishes [25]. The discrepancy originates from distinct dominant breakup mechanisms governing directly injected circular jets versus crossflow liquid jets under identical conditions. Our prior investigations [25] established capillary instability as the primary mechanism governing surface wave dynamics in circular liquid jets. In contrast, the current crossflow jet analysis under undisturbed conditions (Section Atomization of Undisturbed Liquid Jet in Crossflow) reveals Rayleigh-Taylor (RT) instability as the dominant surface destabilization mechanism. This departure in governing physics is also quantitatively evidenced by the characteristic wavelength scaling: the RT-dominant wavelength λ R − T ≈ D (jet diameter) emerges as significantly smaller than the capillary instability wavelength λ C a . Such wavelength disparity, coupled with experimental observations of λ R − T -scale wave patterns, fundamentally explains the contrasting interfacial evolution between circular and crossflow jet systems under equivalent perturbation frequencies. Specifically, the amplitude of the surface waves on the transverse jet induced by the disturbance does not exhibit significant attenuation, whereas the amplitude of the surface waves on the circular jet exhibits gradual decay.

Figure 14 illustrates the temporal evolution of the liquid jet’s global structure in Case 4. Distinct surface wave structures remain observable in the immediately adjacent jet segment to the nozzle under this specific excitation frequency. Comparative analysis with Figure 13 indicates that surface waves excited at this characteristic frequency exhibit shorter relative wavelengths. Furthermore, their amplitudes display progressive decay along the jet development, particularly within the near-nozzle region. This behavior aligns closely with previous circular jet observations under equivalent pulsation frequencies [25]. In contrast to the amplitude attenuation observed near the nozzle, surface waves exhibit amplitude recovery in downstream regions ( x / D > 2.5 ) , demonstrating RT instability’s predominance over capillary effects. These transitional zones are explicitly identified in Figure 15 through solid/dashed red enclosures.

Figure 16 documents the temporal evolution of the liquid jet’s global structure in Case 5. While nozzle exit velocity pulsations generate localized perturbations within the proximal jet region ( x / D < 0.5 ) , capillary-driven damping rapidly dissipates these short-wavelength disturbances. Consequently, the jet’s flow morphology beyond ( x / D > 0.5 ) closely matches the unperturbed scenario. Crucially, the temporal evolution maintains correspondence with crossflow jet behavior under undisturbed conditions. Surface waves exhibiting the characteristic “ Λ ”-shape are preserved under this perturbation regime, as identified by black trajectories in Figure 16. This demonstrates crossflow liquid jet breakup’s insensitivity to high-frequency pulsations, consistent with circular jet observations reported by Zhou et al. [25].

The vortical structures inherent to transverse jets (JICF) - encompassing counter-rotating vortex pairs, upright vortices, shear-layer vortices, horseshoe vortices, and wake vortices [37] - serve as critical determinants of liquid jet disintegration, particularly under turbulent jet regimes. These turbulent vortex configurations exert direct control over interfacial instability initiation mechanisms [34]. Given the heightened susceptibility of jet vorticity fields to external perturbations, nozzle-introduced velocity fluctuations considerably modify transverse jet vortex formation dynamics, thereby governing both initial interface destabilization and subsequent breakup processes. This causality necessitates systematic investigation of nozzle perturbation-induced vortex evolution. Among these vortical features, annular and wake vortices play a dominant role in instability initiation, as evidenced by their persistent manifestation in transient flow regimes [38]. Subsequent analysis focuses on perturbation-induced evolution of these two vortex classes.

Figure 17A depicts the instantaneous vortical topology visualized through the Q -criterion, where the velocity field is normalized by the initial jet velocity. The three-dimensional vortex features surrounding the liquid column reveals two dominant structural components: 1) Counter-rotating vortex pairs (CVPs) aligned parallel to the liquid column axis, directly linked to surface breakup mechanisms; and 2) Windward-side vortices governing the onset of interfacial instability. Additionally, as shown in Figure 17A, the spatiotemporal evolution of wake, vertical, and droplet-associated vortices critically modulates both mixing efficiency and atomization dynamics. Notably, analogous vortical configurations observed in supersonic transverse jets [19] reveal mechanistic similarities in atomization processes across flow regimes.

Figure 17 documents the perturbation-modulated vortical structures surrounding the jet under multiple excitation conditions. Nozzle-imposed disturbances induce significant modifications to the vortex formation dynamics. The windward annular vortex structure demonstrates non-monotonic scaling with increasing actuation frequency, exhibiting initial expansion followed by subsequent contraction. At elevated frequencies, the jet’s vortical configuration approximates undisturbed baseline conditions, demonstrating constrained influence of high-frequency perturbations on breakup mechanisms.

Figure 18 compares windward surface morphologies under varying perturbation frequencies. The unperturbed jet exhibits extended breakup length, where symmetry-preserving waves on the windward interface generate an elongated spine-like liquid ligament along the symmetry plane, deflecting toward the crossflow to enhance liquid core extension. Nozzle-imposed perturbations induce non-monotonic breakup length evolution: initial reduction followed by subsequent elongation with increasing frequency. When surface wavelengths approach RT instability’s dominant scale ( f = 0 − 300 k H z ) , intensified perturbation effects progressively suppress both surface undulations and spinal ligament formation, reducing breakup length through periodic wave packet shedding. However, arrow-marked features in Figure 18D reveal trapezoidal wave configurations emerging at f = 600kHz, triggering renewed breakup length enhancement. At ultrahigh frequencies ( f = 1500 kHz), nozzle-induced wave effects diminish, permitting RT instability-dominated breakup with reemerging “ Λ ”-shaped waves. This regime restores symmetry-plane liquid ligament deflection and breakup length characteristics analogous to unperturbed conditions.

Subsequent analysis focuses on central plane ( z = 0 ) morphological evolution to examine disturbance frequency effects on jet trajectory (Figure 19). When the excitation frequency falls below the capillary instability’s dominant wavelength threshold ( f = 66.6 k H z ) , three distinctive alterations emerge in the perturbed jet compared to its unperturbed counterpart: 1) reduction in breakup length, 2) shortening of column breakup timescale, and 3) periodic lateral oscillations of jet trajectory about its mean path. These collectively enhance crossflow-directed ( x -direction) dispersion of ligaments and droplets. Near RT instability’s critical frequency ( f = 300 k H z ) , accelerated surface wave growth induces marked crossflow deflection (Figure 20), resulting in earlier breakup of the jet column. At elevated frequencies ( f = 600 ,   1500 k H z ) , breakup length restoration occurs alongside trajectory convergence toward unperturbed characteristics. This non-monotonic breakup length response - initial decrease followed by recovery - reaches minimum values below RT instability’s characteristic frequency. The underlying mechanism involves disturbance-induced wavelength modulation: Sub-RT-scale wavelengths transform “ Λ ”-shaped waves into linear configurations while suppressing spinal ligaments, yielding abrupt breakup length reduction. Super-RT wavelengths generate trapezoidal wave morphologies that delay breakup. Ultimately, high-frequency perturbations ( f > 1500 k H z ) experience rapid near-nozzle attenuation, restoring both “ Λ ”-shaped wave patterns and breakup characteristics analogous to unperturbed jets.

The trajectory and penetration depth of liquid jets in crossflow (JICF) constitute critical performance metrics, governing fuel distribution patterns within combustion chambers and consequently determining evaporation rates and oxidizer mixing efficiency [3]. This necessitates rigorous quantitative characterization of both jet trajectory and breakup location. Existing predictive formulations, categorized by operating condition regimes such as standard temperature and pressure (STP), high-pressure and standard temperature (STHP), standard pressure and high-temperature (HTSP), and high temperature and pressure (HTP), along with dimensionless parameters including the momentum flux ratio, Weber number, viscosity ratio, and Reynolds number [3], exhibit inherent limitations. Current models split into formulations for unity nozzle discharge coefficients versus sub-unity cases, yet all require experimental calibration for trajectory prediction accuracy. Direct numerical simulations precisely eliminate flow disturbances within nozzle geometries and experimental measurement noise, enabling more reliable revelation of the underlying physical mechanisms governing jet breakup processes. While their idealized boundary conditions require complementary validation through experimental datasets, the high-fidelity numerical results provide critical theoretical foundations for establishing predictive models of jet trajectory evolution.

The trajectory depth of non-disturbed (unperturbed) liquid jets has been comprehensively characterized through well-established empirical correlations in the technical literature. A survey of published models reveals numerous correlations proposed to link trajectory depth with governing parameters including momentum ratio. However, as demonstrated in our cross-comparative analysis (Figure 21), significant discrepancies exist among existing trajectory depth prediction models. This pronounced divergence strongly suggests that trajectory depth cannot be uniquely determined by momentum ratio alone but depends on additional dimensionless parameters whose functional relationships remain poorly understood. Through systematic evaluation against our DNS dataset, the multi-parameter framework proposed by Eslamian et al. [39] emerges as the optimal predictor, exhibiting excellent agreement with numerical simulations in trajectory depth quantification. This model, rigorously selected from existing formulations through comparative analysis, achieves near-perfect alignment with DNS-derived penetration profiles, confirming its superior predictive accuracy. The Eslamian correlation, expressed as Equation 11: y D = 0.191 x D 0.43 q 0.30 R e ch 0.12 R e l 0.14 , (11) where R e ch designates the gas-phase Reynolds number based on channel hydraulic diameter, and R e l corresponds to the liquid-phase Reynolds number normalized by jet diameter. Distinct from conventional formulations, this framework incorporates two critical dimensionless parameters: R e ch and R e l . These findings establish that the momentum flux ratio constitutes only one component of jet trajectory determination, necessitating inclusion of additional scaling parameters - a critical requirement emphasized by Chai et al. [38].

The existing predictive models for jet trajectory depth, as discussed above, are principally applicable to either unperturbed liquid jets or those subjected to high-frequency excitation. However, in low-frequency excitation regimes, our detailed analysis reveals that inlet perturbations at the nozzle exert significant influence on both the modal characteristics of surface wave formation and their subsequent development dynamics. Under such perturbed excitation conditions, the resultant jet trajectories and breakup lengths exhibit pronounced unsteady characteristics and nonlinear behavior dominated by multimodal interactions. Consequently, establishing reliable predictive frameworks for jet trajectory and breakup length in low-frequency excitation scenarios requires systematic investigations encompassing broader parametric ranges of excitation characteristics–a critical research gap that constitutes the focus of our planned investigations.