{"id":"2c5b3a76-b7a0-4e57-ac37-5487265ca52c","arxiv_id":"2501.18397","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":6,"one_line_summary":"A weakly compressible SPH method with adaptive dissipation, limited particle transport, and a constant-y-plus wall treatment achieves convergent k-epsilon RANS simulations of wall-bounded turbulent channel flows.","lead":"This paper adapts a particle-based fluid simulation method (WCSPH) to solve the k-epsilon RANS model for wall-bounded turbulent flows, adding fixes for near-wall particles and a boundary-offset trick for consistent mesh refinement. It reports, for the first time in this family of methods, convergent velocity and turbulence profiles across several benchmark channels, matching experiments and standard CFD.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Convergence claim rests on unverified consistency of the ARD switch (Eq. 32) and near-wall blending (Eq. 40); if these hand-calibrated terms do not vanish with resolution, the reported convergence is tuning-dependent.","rationale":"Agree with the reader: the weakest point is the hand-calibrated numerical terms. The paper has real strengths: code released, straight-channel convergence over four resolutions, comparisons against FDM/FVM/DNS and experiments, and the boundary-offset technique is a thoughtful way to keep y+ fixed. However, the central claim is global, including separated flow. The straight-channel convergence test does not exercise ARD, and the HCD case is where ARD is load-bearing; Figs. 24-26 show large qualitative differences with and without ARD. The weight ws changes from 0.1 to 0.5 when ARD is enabled, which is an explicit recalibration; without a sensitivity sweep or a scaling argument, the reader cannot know whether the converged k profiles are a property of the SPH-RANS discretization or of the chosen constants. The paper also self-reports occasional non-unique solutions from initial relaxation, which further complicates the convergence claim. None of this proves the method is wrong; it means the headline assertion is conditional on a consistency check that is currently missing. A budget-based scaling test would settle it. If the numerical terms scale to zero (or to the wall-model limit), the CONDITIONAL verdict can be upgraded; if not, the convergence claim should be downgraded.","tokens_in":19229,"tokens_out":8395,"duration_ms":80526,"concrete_test":"Re-run the HCD case at Nf=40, 80 and 120 with ARD active and ws=0.5. In post-processing, accumulate time-averaged per-particle contributions to the k-equation budget: (i) the ARD-induced diffusion (the part of Eq. (21) using mu_R when Eq. (32) selects it), and (ii) the weighted P-layer replacement in Eq. (40) subtracted from the SPH velocity-gradient production. Normalize each by the physical production G_k and dissipation epsilon. If the normalized numerical contribution does not decrease consistently as Nf increases from 40 to 120 (or does not asymptote to a fixed wall-model value), then the convergence shown in Figs. 29-30 is not convergence to the RANS equations the paper claims to solve.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central claim is that SPH-RANS now achieves resolution convergence for velocity and turbulent kinetic energy in wall-bounded flows. That claim requires the added numerical mechanisms to be consistent discretizations of the RANS/wall-model equations, or at least to have errors that vanish with particle spacing. This is not established for the two most important modifications. First, ARD (Eq. 32, mu_c = max(mu_R, e_mu_ij)) replaces the physical eddy-viscosity contribution whenever the Riemann-based numerical viscosity mu_R exceeds it. The straight-channel convergence study (Figs. 10-12) deliberately deactivates ARD, so the strongest convergence evidence does not test the scheme that is essential for the separated HCD case. For HCD, the ARD indicator (Eq. 42) is shown to decrease with refinement, but no quantitative budget demonstrates that the ARD contribution to the momentum or k/epsilon equations becomes negligible relative to physical production and dissipation as dp tends to zero. Second, the weighted near-wall gradient compensation (Eq. 40) mixes the wall-model gradient into the production term with ws=0.1 without ARD and ws=0.5 with ARD. The weight is explicitly recalibrated when ARD is switched on, and no sensitivity study or error estimate is given. Because the P-layer is present at every resolution, this blending term is active in the continuum limit; whether it converges to the intended wall-function limit is untested. Without such evidence, the observed agreement across resolutions could reflect the chosen ws and the ARD switch rather than convergence of the SPH discretization.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes a weakly compressible SPH method for solving the k-epsilon RANS equations in wall-bounded turbulent flows. The central methodological contributions are: an adaptive Riemann-eddy dissipation (ARD) switch, a limited transport velocity formulation (LTVF), a particle-based Lagrangian wall model with a weighted near-wall gradient compensation, and a constant-y_p strategy combined with a level-set boundary-offset technique for resolution-independent wall-adjacent placement. The method is validated on a straight channel, mildly- and strongly-curved channels, and a half converging-diverging (HCD) channel, with comparisons against FDM, FVM, DNS, and experimental data. The headline claim is that this is the first SPH-RANS method to achieve good resolution convergence for both mean velocity and turbulent kinetic energy in wall-bounded flows.","tokens_in":19744,"tokens_out":5647,"duration_ms":53267,"significance":"If the convergence claim is substantiated, this is a meaningful advance for particle-based RANS simulation of engineering wall-bounded flows. The paper identifies a real inconsistency between the Lagrangian particle dynamics and the RANS/wall-model representation, and it proposes several systematic remedies rather than a purely ad hoc fix. The external validation base is genuinely useful: the straight channel is compared with DNS, the curved channels with experiments and FVM, and the HCD channel with two independent FVM codes. The release of the implementation in the open-source SPHinXsys repository is also a concrete strength. However, the paper's central claim currently rests on numerical mechanisms whose consistency with the RANS equations at vanishing particle spacing is not demonstrated: the ARD switch and the weighted near-wall gradient compensation contain hand-set parameters and are not subjected to a quantitative consistency or sensitivity analysis. The convergence evidence is also weaker than stated for parts of the validation matrix. The contribution is therefore promising but not yet fully supported.","major_comments":[{"comment":"The resolution-convergence study that supports the paper's headline claim is performed with the ARD switch deliberately deactivated, as stated in the text: 'Since there is no strong flow separation in this case, the ARD technique is not activated for the convergence test.' Yet the separated HCD case is only convergent when ARD is active. The paper therefore demonstrates convergence for the no-ARD variant and stable, resolution-improving behavior for the ARD variant, but it does not demonstrate that the ARD variant is itself a consistent discretization. Please add a quantitative budget showing the ARD contribution to the momentum and k/epsilon equations (Eqs. (5), (7), (8) with Eq. (32)) decreasing relative to the physical production and dissipation terms as dp tends to zero, or run a straight-channel convergence study with ARD active.","section":"§5.1, Figs. 10-12"},{"comment":"The modulation function mu_c = max(mu_R, e_mu_ij) replaces the physical eddy viscosity by the Riemann-based numerical viscosity whenever the latter is larger. This is an ad hoc switch, not a scale-separated correction: mu_R contains rho0 c0 h and the limiter beta_ij, and the text gives no order-of-convergence estimate for the error introduced when mu_R exceeds e_mu_ij. Because the same switch is essential for HCD stability, the convergence claim requires either a proof that the switch does not alter the RANS balance in the continuum limit or a convergence test with the switch active and its contribution explicitly monitored.","section":"§3.1, Eq. (32)"},{"comment":"The weighted near-wall gradient compensation uses ws=0.1 without ARD and ws=0.5 with ARD. The weight is a free parameter that is explicitly recalibrated when ARD is switched on, and no sensitivity study or error estimate is provided. Since P-layer particles exist at every resolution, this term is active in the continuum limit. The paper must show that the compensated gradient asymptotically approaches the wall-model gradient (Eq. (38)) as dp decreases and that the results are not sensitive to ws in a range independent of ARD. Without this, the observed convergence may reflect the chosen weight rather than the consistency of the discretization.","section":"§4.3, Eq. (40)"},{"comment":"For the mildly-curved channel the two resolutions use different yp values (yp=1.59e-3 for Nf=20 and yp=7e-4 for Nf=40), unlike the constant-yp strategy used for the straight and HCD channels. A convergence test with a changing wall-model parameter is not a resolution-convergence test in the same sense. Please either repeat the mildly-curved study at fixed yp or explicitly label the comparison as a combined resolution/wall-model sensitivity study.","section":"§5.2, Fig. 18"},{"comment":"The straight-channel SPH simulation is initialized from the converged FDM solution to accelerate flow development. This weakens the convergence evidence: with the same equations, parameters, and a nearly converged initial condition, the SPH integration might preserve the FDM profile without independently generating it. Please quantify the evolution from the initial state, for example by starting from a uniform or coarse initial condition and showing that the SPH solution reaches the same converged profile, or by reporting the decay of residuals from the FDM initialization, to rule out initialization bias.","section":"§5.1, first paragraph"},{"comment":"The statement that 'occasionally, with different initial particle relaxation distribution, a different vortex shape would appear' indicates that the reported HCD solution is not unique for a fixed resolution and parameter set. Since the HCD case is one of the central convergence demonstrations, please quantify the run-to-run variability (number of relaxation seeds, range of velocity and turbulent-kinetic-energy errors) and show that the convergence trend is larger than the seed-to-seed scatter.","section":"§5.3, final paragraph"}],"minor_comments":[{"comment":"'an level-set-based' should be 'a level-set-based'.","section":"Abstract"},{"comment":"'Duel-criteria time stepping' should be 'Dual-criteria time stepping'.","section":"§2.5 heading"},{"comment":"The dashed line at y/H=0.05 in the legends is not explained in the text; please define it in the caption or in the text.","section":"§5.1, Figs. 10-12"},{"comment":"The notation ||R^0_nabla_phi||_2 should be explicitly defined; Eq. (26) defines R^0_nabla_phi as a vector, but the limiter in Eq. (33) uses it as a scalar magnitude without stating the norm.","section":"§3.2, Eq. (33)"},{"comment":"The relation R=2h-dp and the statement that yp=dp/2 is applied to all particles in the Pext layer would benefit from a one-sentence derivation, since the figure alone does not make clear why a single yp is applied to a layer of finite thickness.","section":"§4.2"},{"comment":"The time-averaging procedure for the SPH results is not specified: the averaging window, number of samples, and stationarity checks should be reported for all quantitative profiles, since the comparisons in Figs. 10-12, 18, 20-21, and 29-30 rely on time-averaged data.","section":"§5.1-§5.3"}],"recommendation":"major_revision","confidential_remarks":null},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is the first SPH-RANS paper I've seen that shows resolution-converged velocity and k for straight, curved, and separated wall-bounded flows against DNS, experiments, and FVM. That is a real advance for the particle-method community. The ingredients—ARD, LTVF with a limiter, the particle-based wall model, and the constant-y_p boundary-offset strategy—are new in this combination, and the four test cases are well chosen. The code is open source in SPHinXsys, which makes the work reproducible and should be credited.\n\nWhat the paper does well: the straight-channel convergence study (Nf = 10 to 80) is carefully set up, and the boundary-offset technique genuinely fixes the y_p inconsistency that makes convergence tests of wall-model RANS methods ambiguous. The HCD case shows that without ARD the method becomes unstable in separated flow, and the ARD indicator does decrease with refinement. The authors also report the occasional non-unique solutions honestly, which is more than many CFD papers do.\n\nWhere I'd push back: the central claim, first convergent SPH-RANS, is not fully established. ARD (Eq. 32) replaces the physical eddy viscosity with numerical viscosity when the Riemann dissipation exceeds it, and the weighted near-wall gradient compensation (Eq. 40) blends the wall-model gradient into production with ws = 0.1 or 0.5. These are hand-calibrated, and no error estimate or grid-convergence analysis is given for them. The strongest convergence evidence (straight channel) deliberately deactivates ARD, so it does not exercise the mechanism that is essential for the separated case. In the mildly-curved channel, y_p is changed between resolutions (1.59e-3 vs 7e-4), which weakens the convergence claim for that case. Also, the straight-channel run is initialized from the FDM solution, so it tests the discretization's ability to preserve a developed profile more than its ability to develop one from scratch. Time-averaging is mentioned but no averaging window or statistical error bars are given, so it is hard to judge whether the reported differences are noise.\n\nNone of this contradicts the paper's own data, and the stress-test concern about ARD and ws not vanishing is fair but not disqualifying: the visual evidence that the ARD indicator decreases and the three-resolution HCD results are convergent. Still, a referee should ask for a quantitative budget of the ARD contribution and a sensitivity study for ws before the convergence claim is taken as solid.\n\nWho it is for: SPH and particle-method researchers working on turbulence and FSI. It deserves serious peer review; I would send it out with a request for those missing consistency checks rather than desk reject.","headline":"A serious SPH-RANS method paper with genuine benchmark progress; the convergence claim is plausible but rests on hand-calibrated switches that are not shown to vanish, so it needs revision rather than rejection.","tokens_in":20127,"tokens_out":1851,"would_cite":true,"duration_ms":16593,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that a particle-based SPH solver for the k–epsilon RANS equations achieves, for the first time, resolution-converged mean velocity and turbulent kinetic energy in wall-bounded turbulent flows.","keywords":["Smoothed particle hydrodynamics","RANS turbulence model","k-epsilon model","Wall-bounded turbulent flow","Wall model","Lagrangian particle method","Convergence study","Turbulent kinetic energy over-prediction"],"falsifier":"Repeat the straight-channel convergence sweep beyond Nf = 80 with ARD active and w_s fixed at 0.5, then compute the difference of k and U from the FDM and DNS references; if the difference stops shrinking or the wall-nearest k bias reappears, the claimed convergence owes to the tuned blending rather than to the resolved physical balances.","tokens_in":19025,"feed_emoji":"🌊","tokens_out":6542,"duration_ms":61270,"temperature":0.7,"pith_summary":"The paper sets out to make Smoothed Particle Hydrodynamics (SPH) usable for wall-bounded turbulent flow by solving the two-equation k–epsilon RANS equations in a Lagrangian particle framework. Its central assertion is that, with a package of corrections, the method converges with resolution: both mean velocity and turbulent kinetic energy settle toward reference values as particles are refined, in straight, curved, and mildly separated channels. If true, this removes a known obstacle, as particle-based RANS has typically over-predicted turbulent kinetic energy and lacked benchmark validation in wall-bounded cases. The authors identify the root causes as over-damping in the core flow, spurious production in the plug-flow region, and the kinematically sharp shear discontinuity at the wall, and they propose targeted fixes for each.","feed_headline":"SPH turbulence solver converges for wall-bounded flow","feed_subtitle":"Adaptive dissipation plus a fixed wall-adjacent distance lets particle RANS match DNS and experiments.","key_machinery":"The load-bearing objects are four interacting corrections. ARD uses $\\mu_c = \\max(\\mu_R, \\tilde{\\mu}_{ij})$ to pick, per pairwise interaction, whichever dissipation is larger, the Riemann numerical viscosity or the eddy viscosity, so stability is supplied where eddy viscosity is low and removed where it would over-damp. LTVF multiplies the transport-velocity position correction by a limiter $\\beta_{tvf} = \\min(m h^2 \\|R^0_{\\nabla\\phi}\\|_2, 1)$, so particles are re-centered only when the consistency residue is actually large. The weighted near-wall gradient compensation replaces the SPH velocity gradient with a blend of the wall-function gradient and the SPH gradient for neighbors in the wall layer, with weight $w_s$ set to 0.1, or 0.5 when ARD is active. The constant $y_p$ strategy sets the wall-adjacent layer height as a model parameter, and the level-set boundary-offset technique shifts the dummy-particle wall inward by $\\delta = y_p - d_p^r/2$ so that $y^+$ stays in the wall-function's effective range as resolution rises.","core_discovery":"The authors claim that the inconsistency between the Lagrangian character of SPH and the RANS wall-model setup is the reason SPH-RANS has not previously converged for wall-bounded turbulence, and that it can be fixed. In the main stream, an Adaptive Riemann-eddy Dissipation (ARD) switches dissipation between the Riemann solver's numerical viscosity and the eddy viscosity, preventing both particle clustering in high-shear layers and over-damping in the core; a Limited Transport Velocity Formulation (LTVF) suppresses spurious turbulent kinetic energy production caused by position corrections in plug-flow regions. Near the wall, a particle-based wall model evaluates shear stress and velocity gradients from tangential velocities, a weighted compensation scheme restores the under-estimated velocity gradient in the layer just above the wall, and a constant y_p strategy with a level-set boundary offset keeps the non-dimensional wall distance fixed as resolution increases. The result, the paper claims, is the first SPH-RANS convergence study in which both velocity and k converge and agree with DNS, finite-volume, and experimental references at engineering-acceptable resolutions.","pith_inferences":["If the claimed convergence survives un-tuned resolution sweeps with the same w_s and no recalibration, the ARD switch effectively behaves like an automated sub-particle dissipation model, which would make ARD worth testing as an implicit LES-type closure for non-RANS SPH.","The constant-y_p device forgoes the boundary offset once FSI is involved; the paper's own limitation note suggests a wall model independent of y^+ would be needed for rigorous FSI convergence, and that is the natural next test.","The same diagnosed inconsistency, namely kinematic shear discontinuity plus particle migration near walls, should appear in other Lagrangian solvers such as MPS, so the LTVF limiter and weighted gradient compensation may transfer there.","The weight w_s being recalibrated from 0.1 to 0.5 when ARD is active hints that the compensation is coupled to the dissipation balance; a resolution- or flow-dependent w_s would be the decisive check of the mechanism."],"forward_implications":["SPH-RANS can be validated on standard wall-bounded benchmarks such as channels, curved ducts, and converging–diverging ducts, reaching the same mean-flow and k answers as mesh-based RANS at comparable coarse resolutions.","The k over-prediction known from earlier particle RANS implementations is attributed to TVF position corrections in plug-flow regions and to particle-vortex production in high-shear regions, and the LTVF and ARD fixes address those mechanisms directly.","Because wall dummy particles are retained, the framework can be extended to turbulent fluid–structure interaction without redesigning the boundary treatment.","The convergence machinery, namely constant y_p, boundary offset, and y^+-consistent refinement, provides a template for convergence studies of other RANS models in particle methods.","The level-set boundary offset makes the approach applicable to complex geometries, not just straight channels."],"supporting_citations":[{"why":"Supplies the DNS reference data for the straight-channel convergence test that the paper's central convergence claim is measured against.","marker":"[53]"},{"why":"Provides the standard k–epsilon turbulence closure and the stepwise wall function used for near-wall boundary conditions.","marker":"[28]"},{"why":"Supplies the kinetic-based wall shear stress formulation and the operative y+ range that motivates the constant-y_p strategy.","marker":"[29]"},{"why":"Forms the base WCSPH discretization via a low-dissipation Riemann solver that the ARD scheme modifies.","marker":"[30]"},{"why":"Provides the consistency-driven transport-velocity formulation whose zero-order consistency residue the LTVF limiter acts on.","marker":"[38]"},{"why":"Supplies the linear limiter that the Limited TVF uses to suppress over-correction in plug-flow regions.","marker":"[42]"},{"why":"Represents the prior SPH k–epsilon approach that required very high resolution and whose unresolved issues the paper builds on.","marker":"[10]"},{"why":"Documents the k over-prediction problem in particle RANS that LTVF and the weighted near-wall compensation are designed to fix.","marker":"[32]"},{"why":"Provides the level-set-based body-fitted particle generation tool that the boundary-offset technique uses for complex geometries.","marker":"[45]"}],"fun_headline_variants":["SPH-RANS converges for wall-bounded flow via adaptive dissipation","First SPH-RANS convergence for velocity and k in channels","Wall turbulence SPH solver now matches DNS after fixes","SPH solver overcomes RANS inconsistency for wall flows","Adaptive SPH turbulence method converges on wall-bounded cases"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The hand-set numerical dissipation and the hand-calibrated blending weight w_s are assumed not to corrupt the RANS mean-flow and turbulence balances as the resolution increases.","fun_headline_variants_meta":{"raw":{"variants":["SPH-RANS converges for wall-bounded flow via adaptive dissipation","First SPH-RANS convergence for velocity and k in channels","Wall turbulence SPH solver now matches DNS after fixes","SPH solver overcomes RANS inconsistency for wall flows","Adaptive SPH turbulence method converges on wall-bounded cases"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000503,"raw_usage":{"total_tokens":2528,"prompt_tokens":1086,"completion_tokens":1442,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":702,"completion_tokens_details":{"reasoning_tokens":1358}},"tokens_in":702,"tokens_out":1442,"duration_ms":12961,"temperature":1.0,"reasoning_tokens":1358,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-09T23:37:43.053522+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Repeat the straight-channel convergence sweep beyond Nf = 80 with ARD active and w_s fixed at 0.5, then compute the difference of k and U from the FDM and DNS references; if the difference stops shrinking or the wall-nearest k bias reappears, the claimed convergence owes to the tuned blending rather than to the resolved physical balances.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the DNS reference data for the straight-channel convergence test that the paper's central convergence claim is measured against."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the kinetic-based wall shear stress formulation and the operative y+ range that motivates the constant-y_p strategy."},{"cited_title":"Zhang, X","cited_arxiv_id":null,"evidence_quote":"Forms the base WCSPH discretization via a low-dissipation Riemann solver that the ARD scheme modifies."},{"cited_title":"The efficient implementation of transport velocity formulation","cited_arxiv_id":"2411.13992","evidence_quote":"Supplies the linear limiter that the Limited TVF uses to suppress over-correction in plug-flow regions."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Represents the prior SPH k–epsilon approach that required very high resolution and whose unresolved issues the paper builds on."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Documents the k over-prediction problem in particle RANS that LTVF and the weighted near-wall compensation are designed to fix."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the level-set-based body-fitted particle generation tool that the boundary-offset technique uses for complex geometries."}],"review_version":1}