# Symbolic regression modelling of the linear growth rates of stationary cross-flow instabilities

Lucas Pascal

► To cite this version:

Lucas Pascal. Symbolic regression modelling of the linear growth rates of stationary cross-flow instabilities. FLOW, applications of fluid mechanics, 2025, 5, 10.1017/flo.2025.10018. hal-05153606

HAL Id: hal-05153606

https://hal.science/hal-05153606v1

Submitted on 9 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Symbolic regression modelling of the linear growth rates of stationary cross-flow instabilities

Lucas Pascal

ONERA/DMPE, Université de Toulouse, Toulouse, France

Email: lucas.pascal@onera.fr

**Received:** 26 November 2024; **Revised:** 15 April 2025; **Accepted:** 5 May 2025

**Keywords:** Boundary layers; boundary layer stability; crossflow instability; symbolic regression; transition to turbulence

## Abstract

Accurately predicting the growth rates of stationary cross-flow instabilities is crucial for understanding transition mechanisms in swept-wing configurations, which can significantly impact aerodynamic performance. This paper introduces a model for the prediction of cross-flow instability growth rates, focusing on both accuracy and ease of implementation. The proposed model consists of straightforward expressions involving key boundary layer quantities. Validation against established methods demonstrates that the new model achieves comparable or superior accuracy in predicting growth rates. Additionally, tests conducted on a three-dimensional (3-D) prolate spheroid show strong alignment with transition lines computed by means of exact linear stability. Overall, this model provides a practical and efficient alternative for accurately predicting cross-flow transitions in complex 3-D geometries, contributing to improved aerodynamic design and analysis.

## Impact Statement

The transition from laminar to turbulent flow has a significant impact on the skin friction of an airplane wing. Accurate modelling of this phenomenon improves aircraft performance predictions and is essential for the design of laminar wings, leading to energy and cost savings. On a standard commercial aircraft, the wings are swept to enhance performance and efficiency in the transonic regime. As a result, the boundary layer flow becomes three-dimensional, and the cross-flow velocity component develops an inflection point, which is a destabilising factor that leads to the growth of cross-flow waves. The growth rate of stationary cross-flow waves is calculated by solving an eigenvalue problem derived from the linearised Navier-Stokes equations. However, for a three-dimensional wing, this computation is computationally demanding. This paper presents a simplified model for the growth rate of stationary cross-flow waves, significantly reducing computation time when predicting the transition location on three-dimensional geometries using linear stability analysis.

## 1. Introduction

In the field of aerodynamics, accurately predicting the onset of transition from laminar to turbulent flow on aerodynamic surfaces remains a critical challenge, particularly in the context of swept-wing configurations where stationary cross-flow instabilities play a dominant role under free-flight conditions (Arnal et al. (2008)). Transition prediction is of paramount importance because the onset of turbulence significantly influences drag and heat transfer, directly impacting aircraft performance, fuel efficiency and overall aerodynamic design. Among the various types of boundary layer instabilities, cross-flow instability, driven by favourable pressure gradients, is particularly influential on swept wings, affecting both the stability characteristics and the eventual transition location.

© The Author(s), 2025. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Traditionally, predicting cross-flow-induced transition has relied on linear stability theory (LST) calculations, which can accurately determine the linear amplification ratio $e^N$ in the boundary layer (Arnal et al. (2008)). While highly reliable, these methods are computationally intensive and require intricate knowledge of boundary layer profiles at each point on the surface. This is particularly challenging for practical applications, where computational efficiency and ease of implementation are crucial. As a result, empirical and semi-empirical models have been developed to approximate cross-flow instability growth rates without the full complexity of LST, aiming to balance accuracy with computational practicality. In the case of Tollmien–Schlichting (TS) waves, Arnal et al. (1984), Drela and Giles (1987) and Zafar et al. (2021) proposed directly modelling the N-factor curve. Other models provide the growth rate of TS instabilities; see Arnal (1989), Crouch et al. (2001) and Perraud et al. (2009). Due to the simplicity and reduced computational cost of some of these methods, they can be directly implemented in a RANS (Reynolds-Averaged Navier-Stokes) solver; see Coder and Maughmer (2014), Bégou et al. (2017) and Pascal et al. (2020). Dagenhart (1981) developed a model to estimate the growth rate of stationary cross-flow instabilities based on stability charts for ten velocity profiles over Pfenninger’s 970 airfoil, using three key boundary layer quantities. Another model, derived by Perraud et al. (2009), links the growth rates of stationary cross-flow instabilities with the generalised inflection point characteristics of cross-flow velocity profiles. Although effective on the NLF(2)-415 airfoil (Dagenhart and Saric (1999)), Dagenhart’s model relies on interpolations from stored tables, which complicates implementation. Finally, Arnal et al. (1984) and Langtry et al. (2015) derived two criteria predicting the transition location directly, based respectively on the transverse displacement thickness and on the so-called cross-flow strength derived from helicity.

While earlier studies derived growth-rate models manually, more recent approaches rely on machine learning applied to automatically generated datasets. For example, Zafar et al. (2021) applied a recurrent neural network to a database comprising 33,000 boundary layer flows from 53 airfoils to estimate the N-factor of TS instabilities. Additionally, Rouviere et al. (2023) presented a neural network model based on a dataset with 750 elements to determine the extra amplification of TS instabilities due to surface defects. Despite their effectiveness, neural networks are often difficult to interpret, share and implement across different codes. Conversely, symbolic regression offers interpretable, shareable models (Cranmer (2023)). This paper leverages symbolic regression to create models for calculating the growth rate of stationary cross-flow instabilities, aiming to provide simplicity and accuracy for practical applications. A large database of stability charts, generated from three different swept airfoils over a range of flow conditions, underpins this model. Validation is performed by comparison with existing methods on both swept-wing and fully three-dimensional geometries.

In section 2, LST is briefly introduced, and existing models are presented. The database generated for the development of the new model is described in section 3, and the actual derivation of the model is outlined in section 4. Section 5 presents validations of the model on both two-dimensional and three-dimensional geometries. Finally, concluding remarks are provided in section 6.

## 2. Linear stability theory and transition prediction for steady cross-flow waves

Linear stability theory (see Reed et al. (1996) for an extensive review) is a powerful tool for understanding the transition process. The local linear stability equations are derived from the linearised Navier–Stokes equations under the parallel flow assumption, where the perturbation is expressed as a wave ansatz $\vec{q}'(\vec{x},t) = \vec{q}(n) \exp(i(\alpha s + \beta z - \omega t))$. Here, $s$ is the curvilinear abscissa following the velocity at the boundary layer edge, $n$ is the wall distance and $z$ is the coordinate in the direction perpendicular to the boundary layer edge streamline (see Figure 1).

To study stationary cross-flow instabilities, the wave radial frequency $\omega$ is set to zero and the transverse wavenumber $\beta$ is treated as a fixed real parameter. The complex streamwise wavenumber $\alpha$ is computed as the eigenvalue of the resulting dispersion equation. Of particular interest for transition prediction is the growth rate $\sigma = -\mathcal{J}(\alpha)$, which is used in the $e^N$ method. In this method, the total amplification is evaluated by integrating the growth rate along the direction of the group velocity. In practice, the group velocity direction is very close to the direction of the streamline at the boundary
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Definition of the boundary layer coordinate system (s, n, z) and the associated velocities. The thick solid line represents the streamline at the boundary layer edge.</figcaption></figure>

layer edge (Arnal (1994)), and as a result, the integration is performed along the s-direction. For cross-flow instabilities, the most unstable mode occurs at a wave angle $\psi = \arctan(\beta/\mathcal{R}(\alpha))$ between 85° and 89° (Arnal (1994)), which roughly corresponds to the z direction.

In this paper, the N factors $N_β$ and $N_{max}$ are used for transition prediction. They are defined as

$$ \begin{cases} N_{\beta}(s) = \max_{\beta} \int_{0}^{s} \sigma(\xi, \beta) \, d\xi, \\ N_{\max}(s) = \int_{0}^{s} \sigma_{\max}(\xi) \, d\xi \quad \text{where} \quad \sigma_{\max}(\xi) = \max_{\beta} \sigma(\xi, \beta). \end{cases} \qquad (1) $$

Hence, $\sigma$ must be computed for a wide range of wavenumber $\beta$ and at different locations along the boundary layer edge streamline. For TS waves, rapid transition is observed once the nonlinear saturation stage is reached. Although this does not hold true for cross-flow instabilities, linear stability analysis remains valuable in practice for computing the transition location (Arnal et al. (2008)). In practice, the transition location is identified as the point where a given threshold $N^T$ is reached. While the use of $N_β$ is more common, Srokowski and Orszag (1977) observed that both methods yield a comparable spread of $N^T$ among the experimental cases published by Boltz et al. (1960) (this topic is further discussed in Appendix A.1).

## 2.1 Existing models

In this section, five existing models for stationary cross-flow are explained. Each of these models is based on distinct boundary layer characteristics.

### 2.1.1 Dagenhart (1981)

In the coordinate system (s, n, z), the cross-flow velocity $w$ (along z) is zero at the wall, reaches a maximal value $w_{max}$ at $n_{wmax}$, and decreases to zero at the edge of the boundary layer (see Figures 1 and 2). Pfenninger (1977) defines $\delta_{10}$ as the location between $n_{wmax}$ and the boundary layer edge where the cross-flow velocity is $0.1w_{max}$ (see Figure 2). The following Reynolds number is then defined: $Re_{\delta_{10}} = \delta_{10}w_{max}/\nu_e$, where $\nu_e$ is the kinematic viscosity at the boundary layer edge.

These parameters are central to the model developed by Dagenhart (1981). Dagenhart computed stability diagrams $\sigma(Re_{\delta_{10}}, \beta \times \delta_{10}) \times \delta_{10}$ for ten velocity profiles characterised by their values<fnref n="1" />

<footnote n="1">¹The ratio of $n_{wmax}$ over $\delta_{10}$ is called the cross-flow boundary layer shape factor.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Cross-flow velocity profile and definitions of $n_{wmax}$, $n_{gip}$, $\delta_{10}$, $w_{max}$ and $u_{gip}$.</figcaption></figure>

($n_{wmax}/\delta_{10})_i$ and ($w_{max}/U_e)_i$ ($i \in [1, 10]$), where $U_e$ is the velocity at the boundary layer edge. For the ith stability chart, the associated critical Reynolds number ($Re_{\delta_{10},cr}$)$_i$ is computed as a function of ($n_{wmax}/\delta_{10})_i$. The expression for the critical Reynolds number $Re_{\delta_{10},cr}$ is given by

$$Re_{\delta_{10},cr} = 0.4329 \left( \frac{n_{wmax}}{\delta_{10}} \right)^{-3.436} + 38.96. \quad (2)$$

Equation (2) (Xu et al. (2019)) approximates the data presented in Dagenhart (1981, Figure 12).

To compute the growth rate $\sigma$ of a cross-flow instability with wavenumber $\beta$ for a boundary layer profile characterised by $\delta_{10}$, $n_{wmax}$ and $Re_{\delta_{10}}$, the following procedure is applied by Dagenhart:

(i) Identify the stability chart of index $i$ corresponding to the velocity profile whose ($Re_{\delta_{10},cr}$)$_i$ is closest to $Re_{\delta_{10},cr}$.

(ii) Define the adjusted Reynolds number

$$\tilde{Re}_{\delta_{10}} = Re_{\delta_{10}} + (Re_{\delta_{10},cr})_i - Re_{\delta_{10},cr}.$$

(iii) Read $\hat{\sigma}^*(\tilde{Re}_{\delta_{10}}, \beta\delta_{10})$ from the selected stability chart.

(iv) Finally, calculate the growth rate $\sigma$ using the relation

$$\sigma \delta_{10} = \hat{\sigma}^* \frac{\left( \frac{w_{max}}{U_e} \right)}{\left( \frac{w_{max}}{U_e} \right)_i}.$$

### 2.1.2 Perraud et al. (2009)

Given the inflectional origin of cross-flow instability (see Figures 1 and 2), Casalis and Arnal (1996) derived a model for unsteady cross-flow instability based on the characteristics of the generalised inflection point of the velocity profile $u_\psi$ in a plane rotated by an angle $\psi$ relative to the streamline direction. In this model, $n_{gip}$ is the wall distance at which $\frac{d}{dn}\left(\rho\frac{du_\psi}{dn}\right) = 0$, $U_{gip} = u_\psi(n = n_{gip})/U_e$ and $P_{gip} = \left[n/U_e\frac{du_\psi}{dn}\right]_{n=n_{gip}}$ (see Figure 2 for $\psi = 90^\circ$).

A similar approach was later extended by Perraud et al. (2009) for stationary cross-flow instabilities. The model is expressed as

$$\sigma(\psi) = \sigma_{\infty} \left( 1 - \frac{R_c}{Re_{\delta_{1,i},\psi}} \right). \quad (3)$$
In this model, $R_c$ depends on $P_{gip}$ while $\sigma_\infty$ is a function of $U_{gip}$. The Reynolds number $Re_{\delta_{1,i},\psi}$ is defined using the incompressible displacement thickness of the velocity profile $u_\psi$.

The corresponding N-factors can then be computed as

$$ \left\{ \begin{array}{l} N_{\psi}(s) = \max_{\psi} \int_{0}^{s} \sigma(\xi, \psi) \, d\xi, \\ N_{\max}(s) = \int_{0}^{s} \sigma_{\max}(\xi) \, d\xi \quad \text{where here } \sigma_{\max}(\xi) = \max_{\psi} \sigma(\xi, \psi). \end{array} \right. \qquad (4) $$

Here, $N_\psi$ is closely related to $N_\beta$ because 'a fixed value of $\beta$ is associated with a practically constant value of $\psi$' (Arnal et al. (2008)).

### 2.1.3 Xu et al. (2019)

From stability computations on Falkner-Skan-Cooke boundary layer profiles, Xu et al. (2019) derived the following model for the cross-flow instability growth rate:

$$ \sigma_{\beta}\delta_{10} = 2.128 \left(\frac{w_{max}}{U_e}\right)^{1.07} \frac{n_{wmax}}{\delta_{10}} \left(1 + \left|\frac{n_{wmax}}{\delta_{10}} - 0.35\right|^{1.5}\right). \qquad (5) $$

Here, $\sigma_\beta$ is defined such that

$$ N_{\beta}(s) = \int_{0}^{s} \sigma_{\beta}(\xi) \, d\xi. \qquad (6) $$

### 2.1.4 Arnal et al. (1984) and Langtry et al. (2015)

Unlike the previous models, the models by Arnal et al. (1984) and Langtry et al. (2015) do not predict the growth rate of cross-flow instabilities but instead directly determine the transition location:

* Langtry et al. (2015) developed a RANS transition model that accounts for cross-flow transition. This model relies on the so-called non-dimensional cross-flow strength defined as $H_{cf} = nHe/U$, where He is the helicity given by $He = \|\vec{U} \cdot \vec{\Omega}\|/U$. Here, $\vec{U}$ is the local velocity (and U its magnitude) and $\vec{\Omega}$ is the vorticity.

* The C1 (Arnal et al. (1984)) criterion defines the transition location based on a threshold value of $Re_{\delta_2}$ where $\delta_2$ is transverse displacement thickness.

## 3. Building a database of stationary cross-flow characteristics

In this section, the workflow used to compute a database of boundary layer profiles on swept wings and their associated stability characteristics with respect to stationary crossflow (CF) waves is presented. The diagram in Figure 3 illustrates how the process operates.

This workflow is quite similar to the optimisation process described by Sudhi et al. (2023) for transonic airfoil design. The inputs of the workflow are the airfoil geometry and aerodynamic parameters: Reynolds number ($Re_\infty$), Mach number ($M_\infty$), (negative) angle of attack $\alpha$ and sweep angle $\varphi$. As outputs, the workflow generates the boundary layer profiles on the upper surface and the corresponding stability diagrams.

### 3.1 Chosen parameters

The ONERA-D, NLF(2)-415 and NACA 64$_2$ A015 airfoils were selected to build the database. These three airfoils are known for their relevance in studying cross-flow transition (see respectively the references Schmitt and Manie (1979), Dagenhart and Saric (1999) and Boltz et al. (1960)). For each airfoil, six sets of aerodynamic conditions were chosen, as detailed in tables 1, 2 and 3. Different angles of attack
<table><caption>**Table 1.** Aerodynamic conditions for the ONERA-D airfoil</caption><thead><tr><th></th><th colspan="5">ONERA-D</th></tr><tr><th></th><th>#0</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th></tr></thead><tbody><tr><th>Re<sub>∞</sub> × 10<sup>-6</sup></th><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><th>α[°]</th><td>-2.0</td><td>-2.0</td><td>-2.0</td><td>-6.0</td><td>-6.0</td><td>-6.0</td></tr><tr><th>φ[°]</th><td>40.0</td><td>60.0</td><td>80.0</td><td>40.0</td><td>60.0</td><td>80.0</td></tr></tbody></table>

<table><caption>**Table 2.** Aerodynamic conditions for the NLF(2)-415 airfoil</caption><thead><tr><th></th><th colspan="5">NLF(2)-415</th></tr><tr><th></th><th>#0</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th></tr></thead><tbody><tr><th>Re<sub>∞</sub> × 10<sup>-6</sup></th><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td></tr><tr><th>α[°]</th><td>-2.0</td><td>-2.0</td><td>-2.0</td><td>-4.0</td><td>-4.0</td><td>-4.0</td></tr><tr><th>φ[°]</th><td>40.0</td><td>60.0</td><td>80.0</td><td>40.0</td><td>60.0</td><td>80.0</td></tr></tbody></table>

<table><caption>**Table 3.** Aerodynamic conditions for the NACA 64₂ A015 airfoil</caption><thead><tr><th></th><th colspan="5">NACA 64₂ A015</th></tr><tr><th></th><th>#0</th><th>#1</th><th>#2</th><th>#3</th><th>#4</th><th>#5</th></tr></thead><tbody><tr><th>Re<sub>∞</sub> × 10<sup>-6</sup></th><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td><td>5.0</td></tr><tr><th>α[°]</th><td>-1.5</td><td>-1.5</td><td>-1.5</td><td>-1.5</td><td>-1.5</td><td>-1.5</td></tr><tr><th>φ[°]</th><td>30.0</td><td>40.0</td><td>50.0</td><td>60.0</td><td>70.0</td><td>80.0</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Database generation chain.</figcaption></figure>

(α) and sweep (φ) were selected to broaden the range of the cross-flow profiles. These angles were chosen such that stationary cross-flow instabilities develop on the upper surface of the airfoil. No more than six sets of aerodynamic conditions per airfoil were selected to limit the size of the database. The chord Reynolds number value is not crucial for deriving the model as the boundary layer profiles are rescaled when computing the stability diagram (see section 3.4). Therefore, $Re_{\infty} = 5 \times 10^6$ was chosen for all cases. The Mach numbers are selected within the incompressible regime, specifically $M_{\infty} \in [0.05, 0.22]$. This choice is not critical, as stationary cross-flow instabilities exhibit nearly identical behaviour in both incompressible and transonic flows (Arnal (1994)).

## 3.2 Inviscid pressure distribution computation

The pressure distribution is computed using the ISES software (Drela and Giles (1987)) and the simple sweep theory (Meier (2010)). Unlike the approach used by Sudhi et al. (2023), ISES is run with the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The NLF(2)-415 airfoil, the grey region indicates the locations where stability diagrams are computed.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The NACA 64$_2$ A015 airfoil, the grey region indicates the locations where stability diagrams are computed.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The ONERAD airfoil, the grey region indicates the locations where stability diagrams are computed.</figcaption></figure>

inviscid model, where only the two-dimensional (2-D) Euler equations are solved using the finite volume method without coupling to the integral boundary layer equations. The simple sweep theory is then applied to extend the 2-D computation to the 2.5-D flow, allowing the pressure distribution around an infinite swept wing to be calculated.

### 3.3 Solving of the boundary layer equations

Following the infinite swept-wing assumptions (i.e. invariance is assumed along the span), the velocity at the boundary layer edge is derived from the pressure distribution. The dimensionless results computed by ISES are dimensionalised using the following upstream temperature and pressure: $T_∞ = 300$K and $P_∞ = 101325$Pa. This velocity distribution is then used as input for the boundary layer equation solver 3C3D (Houdeville (1992)) (in contrast, Sudhi et al. (2023) used the COCO solver). The 3C3D solver employs the method of characteristics and finite difference discretisation to solve the 3-D boundary layer equations. Although the solver can handle laminar, transitional and turbulent flows, only laminar flow computations are considered in this study. The equations are solved using a marching method, starting from the leading edge. However, the solver cannot account for upstream characteristics, meaning it halts when separation occurs.

### 3.4 Local linear stability computations

Finally, the 1,260 boundary layer profiles computed by 3C3D are passed to ONERA's in-house linear stability solver MAMOUT (Brazier (2015)) (the linear stability LILO solver is used in the work of Sudhi et al. (2023)). The local linear stability equations are discretised by means of a compact finite differences method.

For each case presented in tables 1, 2 and 3, the stability diagrams $\sigma^*(Re_{\delta_{10}}, \beta^*)$ (with $\sigma^* = \sigma \times \delta_{10}$ and $\beta^* = \beta \times \delta_{10}$) are computed at every available station along the wing. For 23 boundary layer profiles, no unstable stationary CF mode is found and the automated stability chart computation fails and these profiles are discarded. Figures 4, 5 and 6 show the locations on each airfoil where the stability diagrams are computed.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Stability diagram $\sigma^*(Re_{\delta_{10}}, \beta^*)$, case no. 1 of the NLF(2)-0415 airfoil at $x/c \approx 1.6\%$.</figcaption></figure>

The different values of $Re_{\delta_{10}}$ used to build the stability diagram are obtained by rescaling the boundary layer profiles. This is achieved by multiplying the two velocity components, $u$ and $w$, of the boundary layer profiles by $Re_{\delta_{10}}/\bar{Re}_{\delta_{10}}$, where $\bar{Re}_{\delta_{10}}$ is the value of $Re_{\delta_{10}}$ without rescaling of the boundary layer profiles.

As an example, the stability diagram obtained for case no. 1 of the NLF(2)-0415 airfoil at $x/c \approx 1.6\%$ is plotted in Figure 7.

**4. Derivation of model by means of symbolic regression**

In this section, the `pysr` library (Cranmer (2023)) is used to derive models for $\sigma_{max}$ and $\sigma$ from the database built in section 3. The aim is to provide accurate models for computing $N_{\beta}$ and $N_{max}$ that are easier to implement and use than those developed by Dagenhart (1981) and Perraud et al. (2009).

Only two open-source, easy-to-install Python libraries were found for symbolic regression: `gplearn` and `pysr`.<fnref n="2" /> Initial tests showed that `pysr` significantly outperforms `gplearn` for the cases considered here, both in terms of accuracy and computational efficiency. `pysr` is based on a multi-population, multi-evolution evolutionary algorithm and is highly customisable. The standard binary operators addition, subtraction, multiplication and division – are selected, along with the unary operators: logarithm (log), hyperbolic tangent (tanh) and inverse ($x \to 1/x$). Additional regressor settings prevent variables from appearing in exponents and limit the nesting of the operators tanh and log. During symbolic learning, 70 % of the data are randomly selected for training, and the remaining 30 % are used for validation. Consequently, the validation set corresponds to interpolation. There is no *a priori* reason to expect the resulting symbolic expressions to yield accurate results in extrapolation.

The regressor inputs are based on cross-flow boundary layer variables used in the models presented in section 2.1: $\delta_{10}$, $n_{wmax}$, $w_{max}$, $Re_{\delta_{10}}$ (models of Dagenhart (1981) and Xu et al. (2019), see sections 2.1.1 and 2.1.3), $H_{cf,max}$ (maximum value of $H_{cf}$ in the boundary layer, defined by Langtry et al. (2015), see section 2.1.4), $n_{gip}$, $U_{gip}$, $P_{gip}$ (model of Perraud et al. (2009), see section 2.1.2; only $\psi = 90^\circ$ is considered here) and $\delta_2$ (model of Arnal et al. (1984), see section 2.1.4). Additional boundary layer quantities are also included: $\delta_1$ (displacement thickness), $\theta_{11}$ (momentum thickness) and $U_e$ (velocity at the boundary layer edge). The composite variable $C_f Re_{\theta_{11}}$ is also used, where $C_f$ is the skin friction coefficient defined as $C_f = \tau_w/(1/2\rho_e U_e^2)$, with $\rho_e$ the density at the boundary layer edge, $\tau_w$ the wall shear stress and $Re_{\theta_{11}} = \theta_{11}U_e/v_e$.

<footnote n="2">²Table 1 of (Cranmer (2023)) compares 11 libraries for symbolic regression.</footnote>
These quantities are combined to form 21 dimensionless variables: $H_{cf,max}$, $P_{gip}$, $C_f Re_{\theta_{11}}$, $\delta_2/\delta_1$, $\theta_{11}/\delta_1$, $\theta_{11}/\delta_2$, $n_{wmax}/\delta_1$, $n_{wmax}/\delta_2$, $n_{wmax}/\theta_{11}$, $\delta_{10}/\delta_1$, $\delta_{10}/\delta_2$, $\delta_{10}/\theta_{11}$, $\delta_{10}/n_{wmax}$, $n_{gip}/\delta_1$, $n_{gip}/\delta_2$, $n_{gip}/\theta_{11}$, $n_{gip}/n_{wmax}$, $n_{gip}/\delta_{10}$, $w_{max}/U_e$, $u_{gip}/U_e$ and $u_{gip}/w_{max}$. Redundancy is expected among these cross-flow variables, but the parsimony constraint imposed by the symbolic regression algorithm ensures that only the most relevant parameters are retained. Moreover, while variables such as $\theta_{11}/\delta_2$ could be computed by combining $\theta_{11}/\delta_1$ and $\delta_2/\delta_1$, doing so would incur a higher cost under the parsimony constraint in pysr than directly including $\theta_{11}/\delta_2$.

As the pysr algorithm involves randomness, each regression is run 20 times and the best-performing model is retained.

## 4.1 Expression for $\sigma_{max}$

The variable $Re_{\delta_{10}}$ is added to the input dataset, and $\sigma_{max}^* = \sigma_{max} \times \delta_{10}$ is set as the output. The latter is obtained by extracting the maximum value of $\sigma^*$ from the stability diagram at each $Re_{\delta_{10}}$ (see Figure 7).

The resulting input database consists of 22 variables and 25,604 entries.

The following expression is returned by pysr:

$$ \sigma_{max} = \frac{H_{cf,max}}{\frac{n_{gip}}{n_{wmax}} + 1.2110786 + \frac{488.197}{1.935369Re_{\delta_{10}} - 96.19734196179}} \quad (7) $$

As expected from Figure 7, pysr returns an increasing function with respect to $Re_{\delta_{10}}$ which asymptotically approaches a limit. This expression yields a critical Reynolds number $Re_{\delta_{10},cr} = 49.7$. This constant value aligns well with Figure 12 of Dagenhart (1981) given the values of $\delta_{10}/n_{wmax} \in [0.30, 0.53]$ found in the database computed in section 3. The parameter $\sigma_{max}$ is found to depend solely on boundary layer quantities related to the cross-flow velocity profile.

The aim of this paper is to develop and validate a model for the growth rate of stationary cross-flow instabilities, rather than to conduct a sensitivity analysis. Therefore, all digits provided by pysr are retained.

Since the derived model consists of a relatively simple expression for the growth rate of stationary cross-flow instabilities, it could be combined with the amplification factor transport method of Coder and Maughmer (2014) for predicting transition in RANS CFD (Computational Fluid Dynamics) software. This integration could be implemented in a CFD solver capable of computing quantities along wall-normal lines (Cliquet et al. (2008), Pascal (2023)).

On the one hand, the proposed model for $N_{max}$ requires:

(i) Computing the boundary layer parameters $H_{cf,max}$, $n_{gip}$, $n_{wmax}$ and $Re_{\delta_{10}}$ at each location; and

(ii) Evaluating Equation (7) at each location and computing the integral (1).

On the other hand, exact linear stability analysis requires:

(i) Solving an eigenvalue problem to compute $\sigma(\xi, \beta)$ at each location for multiple $\beta$ values (typically $\ge 12$); and

(ii) Evaluating Equation (1).

Therefore, the proposed model would significantly reduce the computation time, as demonstrated in section 5.2.

## 4.2 Expression for $\sigma$

Compared with $\sigma_{max}$, modelling $\sigma$ is more complex. An initial attempt was made to build a database using each point from every stability diagram computed in section 3. However, the formulas generated
by pysr were unsatisfactory in terms of accuracy. Consequently, it was necessary to provide the model with a priori knowledge. The stability diagram in Figure 7 is characteristic of inflectional instabilities. For such cases, Casalis and Arnal (1996) modelled the growth rate using two half-parabolas. In this study, however, the wavelength is chosen as the parameter for the parabolas rather than the Reynolds number, and the proposed model is

$$ \sigma(Re_{\delta_{10}}, \beta) = \sigma_{\max} \left( 1 - \left[ \frac{\beta - \beta_M}{\beta_k - \beta_M} \right]^2 \right)^2, \quad \beta_k = \beta_0 \text{ if } \beta < \beta_M \text{ else } \beta_1. \quad (8) $$

While $\sigma_{\max}$ is already known (see Eq. (7)), additional regressions are required to model $\beta_0$, $\beta_1$ and $\beta_M$. Following the approach of Casalis and Arnal (1996), these three wavenumbers were initially sought in the form

$$ \beta_X \times \delta_{10} = a_X + b_X / Re_{\delta_{10}}, \quad (9) $$

where the subscript X denotes either 0, 1 or M. However, it was found that better agreement with the exact LST is achieved by setting $a_0 = b_M = 0$.

For each of the 1,237 available stability diagrams, an initial set of curve fits is performed to determine $\beta_M$, $\beta_0$ and $\beta_1$ at each $Re_{\delta_{10}}$. A second set of curve fits is then conducted to compute the parameters $b_0$, $a_1$, $b_1$ and $a_M$ for each stability diagram. In 86 cases, the regression coefficient $R^2$ was below 0.95, and these diagrams were discarded. Finally, the parameters $b_0$, $a_1$, $b_1$ and $a_M$ were regressed using pysr on a database containing 1,151 entries. The following formulas were obtained:

$$ \left\{ \begin{aligned} a_M &= 1.0431601 + (3.2593474 - 0.007606023 \frac{n_{wmax}}{\delta_2}) \tanh\left(\left[\frac{n_{gip}}{\delta_2} \frac{u_{gip}}{U_e}\right]^{4.5976334}\right), \\ b_0 &= (3.8361076703356 \frac{\delta_{10}}{\theta_{11}} + 270.70172 \tanh(a_M)^{735.47253}) \frac{\log\left(\frac{\delta_{10}}{n_{wmax}}\right)}{\tanh\left(\frac{n_{gip}}{\delta_{10}}\right)}, \\ a_1 &= \frac{0.0049740043b_0}{\tanh(C_f Re_{\theta_{11}})} - 2.93882 \frac{n_{gip}}{n_{wmax}} + 8.263331, \\ b_1 &= 8.446098 \left[ \left( \frac{\delta_{10}}{n_{wmax}} \right)^{1.7631726} + \left( \frac{n_{wmax}}{\delta_{10}} \frac{\delta_1}{\theta_{11}} \right)^{4.8353176} \right] \\ & \quad - 79.171326 (a_1 + C_f Re_{\theta_{11}}) - 53.115273 + 656.5801 \frac{n_{wmax}}{\delta_{10}}. \end{aligned} \right. $$

While the expressions for $b_0$, $a_1$, $b_1$ and $a_M$ primarily depend on boundary layer cross-flow quantities, the streamwise-related variables $\delta_1$, $\theta_{11}$ and $C_f Re_{\theta_{11}}$ also appear in the formulas. The cross-flow shape factor $n_{wmax}/\delta_{10}$, as defined by Dagenhart (1981), appears four times. In contrast to Eq. (7), no dependency on $H_{cf,max}$ is observed.

# 5 Validation

## 5.1 The 2.5-D infinite swept wing: database airfoils

In this section, equations (7) and (8), derived earlier, are validated for infinite swept-wing configurations. The validation cases are constructed using the same tool chain as in section 3 (illustrated in Figure 3). While the same three airfoils are employed, the aerodynamic parameters are selected from existing experimental results published in Schmitt and Manie (1979), Dagenhart and Saric (1999) and Boltz et al. (1960), and are listed in table 4 (with Mach numbers in the range [0.05, 0.07]).

In tables 5, 6 and 7, the values of $N_{max}$ and $N_{\beta}$ computed by integrating the growth rate using equations (7) and (8) are compared with exact linear stability results from MAMOUT. These comparisons are made at locations where either $N_{max}$ or $N_{\beta}$ exceeds 3.0. For each airfoil, the maximum and mean absolute and relative errors are reported. Additional comparisons are made with the methods of Dagenhart (1981) and Perraud et al. (2009).<fnref n="3" /> For $N_{\beta}$, comparisons were also made with the model of Xu et al. (2019);

<footnote n="3">³For this method, only the computation of $N_{max}$ was implemented; therefore, $N_{\beta}$ is not reported.</footnote>
<table><caption>**Table 4.** Aerodynamic parameters for validation computations</caption><thead><tr><th colspan="3">ONERA-D</th><th colspan="3">NLF(2)-415</th><th colspan="3">NACA 64<sub>2</sub> A015</th></tr><tr><th>Re<sub>∞</sub></th><th>α[°]</th><th>φ[°]</th><th>Re<sub>∞</sub></th><th>α[°]</th><th>φ[°]</th><th>Re<sub>∞</sub></th><th>α[°]</th><th>φ[°]</th></tr></thead><tbody><tr><td>1.0 × 10<sup>6</sup></td><td>-6.0</td><td>40</td><td>1.92 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td>9.0 × 10<sup>6</sup></td><td>-1.0</td><td>30</td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>-6.0</td><td>50</td><td>2.19 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td>11.0 × 10<sup>6</sup></td><td>-1.0</td><td>30</td></tr><tr><td>1.0 × 10<sup>6</sup></td><td>-6.0</td><td>60</td><td>2.37 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td>6.0 × 10<sup>6</sup></td><td>-2.0</td><td>30</td></tr><tr><td>1.5 × 10<sup>6</sup></td><td>-6.0</td><td>30</td><td>2.73 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td></td><td></td><td></td></tr><tr><td>1.5 × 10<sup>6</sup></td><td>-6.0</td><td>40</td><td>3.27 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td></td><td></td><td></td></tr><tr><td>1.5 × 10<sup>6</sup></td><td>-6.0</td><td>50</td><td>3.73 × 10<sup>6</sup></td><td>-4.0</td><td>45</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 5.** Mean and maximum errors of N-factors, ONERA-D airfoil</caption><thead><tr><th rowspan="2"></th><th colspan="2">N<sub>max</sub></th><th colspan="2">N<sub>β</sub></th></tr><tr><th>Mean error</th><th>Max error</th><th>Mean error</th><th>Max error</th></tr></thead><tbody><tr><td>Eqs. (7) and (8)</td><td>0.05/1.39 %</td><td>0.33/10.5 %</td><td>0.17/4.7 %</td><td>0.22/5.84 %</td></tr><tr><td>Dagenhart (1981)</td><td>0.43/11.1 %</td><td>0.72/15.0 %</td><td>0.33/8.5 %</td><td>0.65/15.1 %</td></tr><tr><td>Perraud et al. (2009)</td><td>0.73/18.7 %</td><td>1.16/23.7 %</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 6.** Mean and maximum errors of N-factors, NLF(2)-415 airfoil</caption><thead><tr><th rowspan="2"></th><th colspan="2">N<sub>max</sub></th><th colspan="2">N<sub>β</sub></th></tr><tr><th>Mean error</th><th>Max error</th><th>Mean error</th><th>Max error</th></tr></thead><tbody><tr><td>Eq. (7) and (8)</td><td>0.04/0.74 %</td><td>0.13/3.11 %</td><td>0.24/4.5 %</td><td>0.42/5.6 %</td></tr><tr><td>Dagenhart (1981)</td><td>0.21/4.3 %</td><td>0.43/9.9 %</td><td>0.09/2.0 %</td><td>0.22/5.2 %</td></tr><tr><td>Perraud et al. (2009)</td><td>0.77/12.4 %</td><td>1.68/15.1 %</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 7.** Mean and maximum errors of N-factors, NACA 64<sub>2</sub> A015 airfoil</caption><thead><tr><th rowspan="2"></th><th colspan="2">N<sub>max</sub></th><th colspan="2">N<sub>β</sub></th></tr><tr><th>Mean error</th><th>Max error</th><th>Mean error</th><th>Max error</th></tr></thead><tbody><tr><td>Eq. (7) and (8)</td><td>0.28/3.9 %</td><td>0.58/16.1 %</td><td>0.23/3.7 %</td><td>0.50/11.6 %</td></tr><tr><td>Dagenhart (1981)</td><td>0.23/2.8 %</td><td>0.56/12.7 %</td><td>0.24/4.7 %</td><td>0.77/19.7 %</td></tr><tr><td>Perraud et al. (2009)</td><td>2.03/20.2 %</td><td>8.91/87.1 %</td><td>N/A</td><td>N/A</td></tr></tbody></table>

however, this model consistently produced values at least twice those given by exact LST, rendering its inclusion in tables 5, 6 and 7 irrelevant.

The proposed models, based on equations (7) and (8), yield satisfactory results for all three airfoils, for both $N_{max}$ and $N_{β}$. The model by Dagenhart (1981) delivers comparable accuracy. The method of Perraud et al. (2009) achieves only moderate accuracy for $N_{max}$. Among all models, the best and worst agreements are observed for the NLF(2)-415 airfoil and NACA 64<sub>2</sub> A015 airfoil, respectively. The excellent performance of the Dagenhart (1981) model on the NLF(2)-415 airfoil was previously demonstrated in (Dagenhart and Saric (1999)).

## 5.2 The 2.5-D infinite swept wing: NACA0012, DTP-A and DTP-B airfoils

While the derived models yield satisfactory results in section 5.1, their generalisability and robustness for infinite swept-wing computations are further assessed by considering airfoils that do not belong to the training database. Three cases relevant to the study of cross-flow instabilities (Arnal et al. (2007), Tokugawa et al. (2005)) are considered in this section:
<table><caption>**Table 8.** Mean and maximum errors on N-factors, NACA 0012, DTP-A and DTP-B airfoils</caption><thead><tr><th rowspan="2"></th><th colspan="2">N<sub>max</sub></th><th colspan="2">N<sub>β</sub></th></tr><tr><th>Mean error</th><th>Max error</th><th>Mean error</th><th>Max error</th></tr></thead><tbody><tr><td>Eqs. (7) and (8)</td><td>0.08/1.22 %</td><td>0.23/7.02 %</td><td>0.20/3.1 %</td><td>0.50/7.4 %</td></tr><tr><td>Dagenhart (1981)</td><td>0.87/8.8 %</td><td>1.91/11.2 %</td><td>0.94/8.6 %</td><td>2.79/16.0 %</td></tr><tr><td>Perraud et al. (2009)</td><td>0.96/11.3 %</td><td>2.55/23.3 %</td><td>N/A</td><td>N/A</td></tr></tbody></table>

(i) NACA 0012 airfoil: $\alpha = -12°$, $\varphi = 40°$ and $Q_\infty = 20m/s$ (upstream velocity). The upstream temperature and pressure are set to $T_\infty = 300K$ and $P_\infty = 101325Pa$, respectively.

(ii) DTP-A airfoil: $\alpha = 0°$, $\varphi = 40°$, $Q_\infty = 36m/s$, $Re_\infty = 2.8 \times 10^6$ and $P_\infty = 101325Pa$.

(iii) DTP-B airfoil: $\alpha = 6°$ (the lower surface is considered), $\varphi = 40°$, $Q_\infty = 70m/s$, $Re_\infty = 3.26 \times 10^6$ and $P_\infty = 101325Pa$.

The validation process follows the same methodology as in section 5.1: the maximum and mean absolute and relative errors on $N_{max}$ and $N_β$ are presented in table 8. Errors are computed by comparing the results of exact linear stability computations performed with MAMOUT against those obtained using equations (7) and (8), as well as the models of Dagenhart (1981) and Perraud et al. (2009). Only regions supporting cross-flow instabilities are retained by considering locations where either $N_{max}$ or $N_β$ exceed 3.0.

The models defined by equations. (7) and (8) yield satisfactory results. The magnitudes of the mean and maximum errors are comparable to those obtained in section 5.1 (see in particular table 5). The derived models outperform those of Dagenhart (1981) and Perraud et al. (2009) across all three cases. In particular, the model of Perraud et al. (2009) is once again found to provide only moderate accuracy.

All computations were performed on a single core of a workstation for the three cases considered in this section. Exact linear stability analysis with the MAMOUT solver required 58 minutes to compute $\sigma$ at each location and for 20 values of $\beta$. Less than a second was then needed to evaluate both integrals $N_{max}$ and $N_β$ using Eq. (1). The simplified model derived in this paper reduced the computation time to 26 s for $N_β$ and 16 s for $N_{max}$.

## 5.3 The 3-D prolate spheroid

Transition on the 3-D prolate spheroid with aspect ratio $a/b = 1.2$ was measured by Kreplin et al. (1985). Flow conditions $Re = 6.54 \times 10^6$, $\alpha = 15°$ were selected to ensure a transition scenario dominated by cross-flow instability (Stock (2006)). Despite the geometric and flow complexity in this configuration, the potential solution can be computed analytically following the approach of Cebeci et al. (1978). The analytically derived velocity distribution is provided as input to 3C3D to compute the boundary layer velocity profiles over the geometry. The surface mesh is generated using step sizes of $1°$ and $0.5°$ along the angular coordinate $\chi$ (with $x = -a \cos(\chi)$) and the polar angle $\phi$ in the $(y, z)$ plane, respectively. To avoid the velocity singularity at the leading edge, the mesh starts at $\chi = \chi_{min} = 0.049$, corresponding to $x/a = -0.9988$.

The growth rate from linear stability analysis is computed at each location using ONERA's in-house stability solver MAMOUT, as well as the models developed in this paper and that of Dagenhart (1981). Growth-rate values are interpolated along 360 streamlines (relative to the potential flow), enabling the computation of N-factors along each. The streamline starting locations are uniformly distributed along the azimuth (i.e. with a step size $\Delta\phi = 0.5°$) at $\chi = \chi_{min}$.

The threshold value $N_β^T = 5.5$ given by Stock (2006) is applied. According to Stock (2006, Figure 15 b), stationary cross-flow instability is dominant within the region $\phi \in [-135°, -80°]$. This region is therefore the focus of the following discussion.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Computed transition locations with eq. (7) (□), with the model of Dagenhart (1981) (△) and with MAMOUT (solid line). Streamlines are plotted as dotted lines.</figcaption></figure>

### 5.3.1 Computation of $N_{max}$

In this section, equation (7) is tested by comparing transition locations obtained using $N_{max}^T = 8.2$. This value is selected based on the ratio $N_{max}^T/N_β^T$ reported by Srokowski and Orszag (1977) (see Appendix A.1). The resulting transition lines are shown in Figure 8 and compared with predictions from the model of Dagenhart (1981) and with exact linear stability results from MAMOUT. Both the model developed in this study and that of Dagenhart (1981) show strong agreement with exact stability results. However, equation (7) predicts a transition line that lies slightly upstream of the exact linear stability prediction, whereas the model of Dagenhart (1981) places the transition line slightly downstream.

### 5.3.2 Computation of $N_β$

In this section, transition lines based on $N_β$ are computed and compared. Figure 9 shows that the transition lines obtained using equation (8) and the method of Dagenhart (1981) align closely with the results of exact linear stability computation. As previously observed, the transition line predicted by equation (8) lies slightly upstream of the MAMOUT prediction, while the transition locations given by the method of Dagenhart (1981) are found both upstream and downstream. Very good agreement is observed in the region of interest with the computations of Krimmelbein (2021), whereas the transition line computed by Stock (2006) appears slightly further downstream. No explanation can be offered for these discrepancies.

## 6. Conclusions

In this study, a new model was derived to estimate the growth rate of stationary cross-flow instabilities, combining ease of use with good accuracy. The model was developed using a comprehensive database of 1,237 boundary layer profiles computed for three representative airfoils (ONERA-D, NLF(2)-415 and NACA 64₂ A015) under a range of aerodynamic conditions. Linear stability characteristics of these profiles were analysed using the ONERA MAMOUT solver, and symbolic regression techniques were applied to derive expressions for the growth rate of cross-flow instabilities. As expected the resulting formulas were found to depend primarily on cross-flow boundary layer quantities.

Validation of the model on infinite swept-wing configurations showed satisfactory results. It performed consistently well across the ONERA-D and NLF(2)-415 airfoils, although slightly higher errors were observed for the NACA 64₂ A015 airfoil. The new model outperforms the method of Perraud *et al*.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Computed transition locations with eq. (8) (□), with the model of Dagenhart (1981) (Δ) and with MAMOUT (solid line). Streamlines are plotted as dotted lines. The transition lines computed by Krimmelbein (2021) and Stock (2006) are plotted as dash-dotted and dashed lines, respectively.</figcaption></figure>

(2009). The method of Dagenhart (1981) yields slightly better accuracy on the NACA 64$_2$ A015 airfoil for $N_β$ and on the NLF(2)-415 airfoil for $N_{max}$. The new model also produces satisfactory results for configurations involving airfoils that are not part of the training database (NACA 0012, DTP-A, DTP-B).

To assess the predictive capability of the derived models on 3-D geometries, the prolate spheroid configuration was considered. Under flow conditions dominated by cross-flow instabilities, comparisons of transition lines based on $N_{max}$ and $N_β$ show that the models developed in this paper closely align with the results from exact linear stability computations. Specifically, the transition lines predicted using the growth-rate equations generally show excellent agreement with the MAMOUT results, with the model-based transition lines lying slightly upstream.

The proposed model offers a simpler and more accurate alternative to the model developed by Perraud et al. (2009), delivering improved precision without added complexity. Compared with the method by Dagenhart (1981), the accuracy is comparable. However, the implementation of the present model is significantly more straightforward, as it eliminates the need for storage and interpolation of precomputed tables, making it highly practical for broader application in transition prediction.

The overall performance indicates that this model offers a robust and efficient tool for predicting cross-flow transition in boundary layers. Future work could focus on implementing the model in RANS solvers, with potential applications in aerodynamic design and optimisation.

**Acknowledgements.** The author acknowledges K. Dotse and V. Mouysset for providing the program used to compute streamlines on the prolate spheroid.

During the preparation of this work the author used ChatGPT and DeepSeek in order to improve language and readability. Moreover ChatGPT was used to generate a skeleton for the abstract. After using this AI tool, the author reviewed and edited the content as needed and takes full responsibility for the content of the publication.

**Data availability statement.** The data and code used in this work are currently proprietary and confidential.

**Funding statement.** The author gratefully acknowledges funding by the ONERA project MATRIX directed by F. Méry.

**Competing interests.** The author declares no conflict of interest.

**Ethical standards.** The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
<figure><img src="image_11.png" /><figcaption>Figure 10. Computed transition lines defined by $N_β = N_β^T$ (solid line) and by $N_{max} = r \times N_β^T$ (dash-dotted line). Streamlines are plotted as dotted lines.</figcaption></figure>

# Appendix

## A.1 Relationship between $N_β$ and $N_{max}$

Srokowski and Orszag (1977) determined both $N_β^T$ and $N_{max}^T$ corresponding to the transition locations published by Boltz et al. (1960) (for the NACA 64_2 A015 airfoil). They found that $N_β^T \in [6.7, 7.2]$, corresponding to a spread of 7.2% around the arithmetic mean while $N_{max}^T \in [9.7, 11]$ which corresponds to a spread of 6.1%.

A similar investigation using experimental results from Dagenhart and Saric (1999) yielded $N_β^T \in [5.6, 6.8]$ and $N_{max}^T \in [6.9, 8.3]$, resulting in relative spreads of 20.0% and 18.5%, respectively.

While both integration strategies are considered valid for transition prediction, $N_β^T$ is typically reported. Here, we attempt to deduce $N_{max}^T$ from $N_β^T$. For each airfoil and set of flow conditions computed in section 5.1, a subset was created with $N_β > 3$, resulting in a dataset of 850 entries with 23 features (the 22 features from section 4.1 plus $N_β$). Using *pysr* as in section 4, we derived the following expression for $N_{max}/N_β$:

$$r = \frac{N_{max}}{N_{\beta}} = (p_{gip} + 1.1400409) \left( \frac{\delta_2}{\theta_{11}} \left( \frac{n_{gip}}{\delta_{10}} - 0.23753823 \right) - 0.09712514 \right) + 1.1903278. \quad (10)$$

Equation (10) is validated on the 3-D prolate spheroid case from section 5.3. In Figure 10, transition lines computed using MAMOUT growth rates are plotted by comparing $N_β$ with $N_β^T$ and $N_{max}$ with $r \times N_β^T$ (only where $N_{max} > N_β^T$). The two transition lines display excellent agreement, demonstrating that using $N_β^T$ along with $r$ allows accurate prediction of the transition location based on $N_{max}$ computation.

# References

* Arnal, D. (1989). Transition prediction in transonic flow. In Symposium Transsonicum III, pp. 253–262, Springer.

* Arnal, D., Casalis, G., & Houdeville, R. (2008). Practical transition prediction methods: subsonic and transonic flows. In VKI Lectures Series Advances in Laminar-Turbulent Transition Modelling.

* Arnal, D., Habiballah, M., & Coustols, E. (1984). Laminar instability theory and transition criteria in two and three-dimensional flow. *La Recherche Aérospatiale*, 2, 125–143.

* Arnal, D., Piot, E., Archambaud, J.-P., Casalis, G., Content, C., Dandois, J., & Colamartino, S. (2007). TRACMIR: TRAnition control by Micron-sized roughness elements.
Arnal, D. (1994). Boundary layer transition: predictions based on linear theory, In In AGARD 793, Special Course on Progress in Transition Modelling, volume 1. AGARD.

Bégou, G., Deniau, H., Vermeersch, O., & Casalis, G. (2017). Database approach for laminar-turbulent transition prediction: Navier-Stokes compatible reformulation. *AIAA Journal*, 55(11), 3648–3660. https://doi.org/10.2514/1.J056018

Boltz, F. W., Kenyon, G. C., & Allen, C. Q. (1960). Effects of sweep angle on the boundary-layer stability characteristics of an untapered wing at low speeds NASA, Technical Report TN-D-338.

Brazier, J.-P. (2015). Modules d'Analyse MOdale Unidimensionnelle avec les polynômes de Tchebychev. Instructions pour l'emploi des programmes Spectre et Moditer Version 6 Technical Report RT 1/24067 DMAE, ONERA.

Casalis, G., & Arnal, D. (1996). Database method-development and validation of the simplified method for pure crossflow instability at low speed ELFIN II European Project, Tech. Rep. 145.

Cebeci, T., Hirsh, R. S., & Kaups, K. (1978). Calculation of three dimensional boundary layers on bodies of revolution at incidence Douglas Aircraft Company, Technical Report MDC J7643.

Cliquet, J., Houdeville, R., & Arnal, D. (2008). Application of laminar-turbulent transition criteria in Navier-Stokes computations. *AIAA Journal*, 46(5), 1182–1190. https://doi.org/10.2514/1.30215

Coder, J. G., & Maughmer, M. D. (2014). Computational fluid dynamics compatible transition modeling using an amplification factor transport equation. *AIAA Journal*, 52(11), 2506–2512. https://doi.org/10.2514/1.J052905

Cranmer, M. (2023). Interpretable machine learning for science with PySR and SymbolicRegression.jl, arXiv preprint arXiv: 2305.01582. https://doi.org/10.48550/ARXIV.2305.01582

Crouch, J., Crouch, I., & Ng, L. (2001). Estimating the laminar/turbulent transition location in three-dimensional boundary layers for CFD applications. In 15th AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.2001-2989

Dagenhart, J. R. (1981). Amplified crossflow disturbances in the laminar boundary layer on swept wings with suction. Technical Report Technical paper, 1902.

Dagenhart, J., & Saric, W. S. (1999). Crossflow stability and transition experiments in swept-wing flow NASA, Technical Report NASA/TP-1999-209344.

Drela, M., & Giles, M. B. (1987). Viscous-inviscid analysis of transonic and low Reynolds number airfoils. *AIAA Journal*, 25(10), 1347–1355. https://doi.org/10.2514/3.9789

Houdeville, R. (1992). Three-dimensional boundary layer calculation by a characteristic method. In Fifth Symposium on Numerical and Physical Aspects of Aerodynamic Flows, Long Beach.

Kreplin, H. P., Vollmers, H., & Meier, H. U. (1985). Wall Shear Stress Measurements on an Inclined Prolate Spheroid in the DFVLR 3m x 3m Low Speed Wind Tunnel, Göttingen DFVLR, Technical Report 222-84.

Krimmelbein, N. (2021). Transition prediction method for three-dimensional flows using linear stability theory', URL https://elib.dlr.de/205660/.

Langtry, R., Sengupta, K., Yeh, D., & Dorgan, A. (2015). Extending the Gamma-rethetat correlation based transition model for crossflow effects. *45th AIAA Fluid Dynamics Conference*, 1–12. https://doi.org/10.2514/6.2015-2474

Meier, H.-U. (2010). *German development of the swept wing 1935–1945*. Reston, VA: American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/4.867552

Pascal, L., Delattre, G., Deniau, H., & Cliquet, J. (2020). Stability-based transition model using transport equations. *AIAA Journal*, 58(7), 2933–2942. https://doi.org/10.2514/1.J058906

Pascal, L. (2023). Implementation of stability-based transition model by means of transport equations. In *AIAA SciTech forum*. American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.2023-1622,

Perraud, J., Arnal, D., Casalis, G., Archambaud, J.-P., & Donelli, R. (2009). Automatic transition predictions using simplified methods. *AIAA Journal*, 47(11), 2676–2684. https://doi.org/10.2514/1.42990

Pfenninger, W. (1977). Laminar flow control laminarization. In *AGARD spec. Course on concepts for drag reduction*. vol. 654, pp. 3.1–3.62. AGARD.

Reed, H. L., Saric, W. S., & Arnal, D. (1996). Linear stability theory applied to boundary layers. *Annual Review of Fluid Mechanics*, 28(1), 389–428. https://doi.org/10.1146/annurev.fl.28.010196.002133

Rouvière, A., Pascal, L., Méry, F., Simon, E., & Gratton, S. (2023). Neural prediction model for transition onset of a boundary layer in presence of two-dimensional surface defects. *Flow*, 3, E20. https://doi.org/10.1017/flo.2023.17

Schmitt, V., & Manie, F. (1979). Ecoulements subsoniques et transsoniques sur une aile a fleche variable. La Recherche Aérospatiale. No. 4.

Srokowski, A., & Orszag, S. (1977). Mass flow requirements for LFC wing design. *Aircraft Design and Technology Meeting*. https://doi.org/10.2514/6.1977-1222

Stock, H. W. (2006). eN transition prediction in three-dimensional boundary layers on inclined prolate spheroids. *AIAA Journal*, 44(1), 108–118. https://doi.org/10.2514/1.16026

Sudhi, A., Radespiel, R., & Badrya, C. (2023). Design exploration of transonic airfoils for natural and hybrid laminar flow control applications. *Journal of Aircraft*, 60(3), 716–732. https://doi.org/10.2514/1.C036968
Tokugawa, N., Takagi, S., Ueda, Y., & Ido, A. (2005). Influence of the external disturbances on natural boundary-layer transition in rectangular wing flows, *24*(6), 629-639.

Xu, J., Han, X., Qiao, L., Bai, J., & Zhang, Y. (2019). Fully local amplification factor transport equation for stationary crossflow instabilities. *AIAA Journal*, *57*(7), 2682-2693. https://doi.org/10.2514/1.J057502

Zafar, M. I., Choudhari, M. M., Paredes, P., & Xiao, H. (2021). Recurrent neural network for end-to-end modeling of laminar-turbulent transition. *Data-Centric Engineering*, *2*, e17. https://doi.org/10.1017/dce.2021.11

<footnote>Cite this article: Pascal L. (2025). Symbolic regression modelling of the linear growth rates of stationary cross-flow instabilities. *Flow*, *5*, E25. https://doi.org/10.1017/flo.2025.10018</footnote>