A heavy body translating in a boundary layer: ‘crash’, ‘fly away’ and ‘bouncing’ responses

Abstract

The study concerns a slender, heavy body moving with streamwise velocity in a boundary layer. Modelling assumptions on body size reduce the governing equations for the body motion to a pair of nonlinear integro-differential equations (IDEs) which displays a wide range of distinguished behaviours, including eventual collision with the wall (‘crash’), escape to infinity (‘fly away’) and repeatedly travelling far from the wall and back again without ever colliding or escaping (‘bouncing’). The paper gives a survey of the variety of behaviour, as well as asymptotic analysis and insight into each category of fluid/body interaction and the conditions under which crash, fly away and bouncing occur.

1. Introduction

The problem being considered is on the dynamic interaction between a thin, heavy rigid body moving through a boundary layer and the surrounding fluid, with the challenge being to predict the body motion. This issue is important in the study of aircraft icing, a phenomenon where ice particles in the boundary layer of an aircraft wing at high altitude may adhere to the surface and inhibit the aerodynamics, in worst cases resulting in serious accidents (Gent, Dart & Cansdale 2000; Schmidt, Young & Benard 2010; Purvis & Smith 2016; Norde 2017). Thus there is significant industrial interest in predicting the motion and behaviour of such particles, which can vary widely in size, shape, speed and angle of incidence, to improve anti-icing systems.

In the aircraft-icing context, the ice particles of concern are in many cases shards or ice crystals and they tend to be isolated, entering the air boundary layer one by one. However, the large range of variables involved requires the problem to be studied for many different bodies in a variety of settings. The studies by Smith & Wilson (2011), Smith et al. (2019)and Smith & Ellis (2010) focus on irrotational flow, corresponding to a body far from the wall, while more recent studies by Smith & Johnson (2016), Smith (2017), Palmer & Smith (2019), Smith & Servini (2019), Palmer & Smith (2021) and Smith & Palmer (2019), as well as the current, have considered bodies in a boundary layer or channel flow for which there is substantial incoming vorticity. These predominantly analytical contributions are the most relevant here. Other very interesting contributions on fluid/body interactions are the experimental ones by Hall (1964), Einav & Lee (1973), Petrie et al. (1993), Schmidt & Young (2009), Wang & Levy (2006) and the direct simulations by Palmer & Smith (2020), Wang & Eldredge (2015), Wu et al. (2018), Eldredge (2008). See also the studies by Gavze & Shapiro (1997), Kishore & Gu (2010), Frank et al. (2003), Loth & Dorgan (2009), Poesio et al. (2006), Yu, Phan-Thien & Tanner (2007) in which particular attention is given to whether a particle will move towards or away from the wall, and those by Dehghan & Basirat Tabrizi (2014) and Portela, Cota & Oliemans (2002) which are concerned with particles in turbulent flows. In previous analytical studies set in a boundary layer, the body is sufficiently small as to be able to enter the near-wall layer in which viscous-inviscid interaction occurs, while the present work will focus on bodies large enough for such effects to be ignored. The body in this work is also taken to be undergoing streamwise translation, a situation allowing for a large range of angles of incidence as the body enters the boundary layer, which in the rest frame of the body means non-zero wall velocity.

The investigation will focus on nonlinear dynamic fluid/body interactions, in which a thin body is nearly but not quite aligned with the surrounding boundary-layer flow, and include vorticity. The Reynolds number and Froude number of interest here are large. The use of the physically relevant assumption that the body is much heavier than the surrounding fluid allows significant analytical progress to be made as the flow is quasi-steady on the time scale of the interaction with the body, reducing the problem to a set of nonlinear integro-differential equations (IDEs). This system has a diverse set of interactive solutions, resulting in many outcomes for the long-term behaviour of the body. The first half of the paper concerns ‘crash’ solutions, i.e. finite time collision (impact) with the wall. This is done both for incoming uniform shear flow and for arbitrary velocity profiles in a boundary layer, with the arbitrary case having already been discussed in brief in Jolley, Palmer & Smith (2021). The latter half of this investigation will focus on ‘lift-off’, i.e. the body departing or otherwise not colliding with the wall; such phenomena have been studied for inviscid settings by Smith & Wilson (2013) and Balta & Smith (2018). It will transpire that the presence of vorticity causes a huge range of intriguing behaviours which are qualitatively different from each other and from those found in previous studies. Particular attention will be given to behaviours termed ‘fly away’, i.e. departure of the body far from the wall, and ‘bouncing’, i.e. repeated departure followed by return (without collision), although other interesting solutions do exist and are discussed briefly.

The layout of the paper is as follows. Section 2 will set out the modelling of the fluid and body motion and derive the nonlinear system to be studied in the remaining sections. Section 3 and an accompanying appendix derive the asymptotic behaviour of the fluid/body interaction in a crash scenario for the cases of constant vorticity and arbitrary incoming velocity profile, respectively. Focusing on the constant vorticity case from then on, § 4 analyses ‘fly away’ while § 5 focusses on ‘bouncing’, for which the modelling is related to fly away, and § 6 provides a brief survey of other types of possible behaviour. Section 7 provides further discussion, including parameter values and the checking of assumptions, followed by conclusions.

2. The fluid–body interaction

The fluid–body interaction and the resulting governing equations are presented in full in the current section. This is to set out a clear basis for the new interactive solutions described in subsequent sections.

The scaling and governing equations of the problem are the same as those in Jolley et al. (2021); however, approximating the oncoming boundary layer flow as initially close to uniform shear will allow significant simplifications to the system in that paper which, in turn, allows analysis of several new solution categories. The body here is described by its mass, moment of inertia and underbody curve. It translates upstream relative to the oncoming flow with speed $u_c$, but we work in its rest frame (hence this corresponds to a positive wall velocity of magnitude $u_c$). The set-up is shown in figure 1. The body is assumed to be both ‘small’ and ‘thin’, by which specifically we mean its length $L$ is much less than 1 and its thickness $\varDelta$ is much less than $L$, where lengths here are non-dimensionalised relative to $\hat {l}$, the representative boundary layer length. The thickness $\varDelta$ is comparable with the thickness of the boundary layer. Here, we will use non-dimensional scaled coordinates throughout, such that $(\hat {x},\hat {y}) = \hat {l}(LX,\Delta Y)$, where $\hat {\ }$ represents a dimensional quantity.

Figure 1. Diagram (not to scale) showing the non-dimensional (and scaled) set-up: a thin body translates upstream in the boundary layer of a much larger body, which acts as a stationary wall – in the rest frame of the body, this results in a positive flow velocity $B$ at the wall. The height of the body's centre of mass (CoM) is $h(T)$ and the angle its chord line makes with the $X$-axis is $\theta (T)$; $T$ denotes scaled time. The incoming velocity profile is $u = u_0(Y) = AY+B$.

The body has two degrees of freedom in its motion: the height of its centre of mass $\Delta h(T)$ and the angle its chord line makes with the $X$-axis, $\Delta \theta (T)$. The aim of the problem in the present study is to find the body motion (i.e. $h(T)$ and $\theta (T)$) under interaction with the surrounding air flow. Hence it is necessary to compute the force and torque exerted on the body by the flow.

The flow here is two-dimensional and incompressible, and non-dimensionalisation with respect to $\hat {l}$, representative boundary layer velocity $\hat {U}$, fluid density $\hat {\rho }_F$ and kinematic viscosity $\hat {\nu }$ gives a Reynolds number $Re = \hat {U}\hat {l}/\hat {\nu }$ which is large for parameter ranges of interest, as discussed in the conclusion of the paper and in Appendix A. Relevant Froude numbers similarly are large. Hence viscous and gravitational forces on the body are to be neglected as a first-go model, an aspect which is discussed further in the appendix. The scales of the problem ($\varDelta = Re^{-1/2} \ll L \ll 1$) allow it to be modelled by the inviscid boundary layer equations (see (2.3)–(2.5) below), which, in turn, allows the pressure forces above the body to be negated also by identification with the free stream. Hence the only force in play is the pressure force underneath the body, leading to the governing equations,

(2.1)$$\begin{gather} Mh_{TT} = \int_0^1 p\, {\rm d}X, \end{gather}$$

(2.2)$$\begin{gather}I\theta_{TT} = \int_0^1 (X-\beta)p\, {\rm d}X, \end{gather}$$

where $T$ is a non-dimensional time coordinate scaled to be order one relative to the body motion, related to its dimensional counterpart by $\hat {t} = (\hat {l}/\hat {U}) \Delta (\hat {\rho }_B/\hat {\rho }_F)^{1/2} T$, and $p = \hat {p}/\hat {\rho }_F\hat {U}^2$ is the pressure field under the body. The non-dimensional mass and moment of inertia of the body are $M = \hat {M}/(\hat {\rho }_B\hat {l}^2L\varDelta )$ and $I = \hat {I}/(\hat {\rho }_B\hat {l}^4L^3\varDelta )$ (where $\hat {\rho }_B$ is the constant density of the body).

According to the time scale derived from (2.1)–(2.2), the body motion occurs slowly in comparison with time variations in the fluid provided that $\hat {\rho }_F/\hat {\rho }_B \ll (\varDelta / L)^2$. The Navier–Stokes equations thus reduce to the inviscid quasi-steady boundary layer equations,

(2.3)$$\begin{gather} uu_X+Vu_Y ={-}p_X, \end{gather}$$

(2.4)$$\begin{gather}0 ={-}p_Y , \end{gather}$$

(2.5)$$\begin{gather}u_X+V_Y = 0, \end{gather}$$

where we require $v(x,y,t) = (\varDelta /L)V(X,Y,T)$ owing to incompressibility. Viscous terms are negligible despite the body width being comparable with that of the boundary layer owing to its small horizontal length: while $u_{yy}/Re \sim 1$, the inertial terms and pressure gradient are of size $1/L$ and are thus significantly larger. At any given time, the curve describing the position of the underbody is given by

(2.6)\begin{equation} Y = F(X,T) = F_u(X)+h(T)+(X-\beta)\theta(T), \end{equation}

where $\beta$ is the $X$-coordinate of the centre of mass, which allows the boundary conditions to be formulated as

(2.7)$$\begin{gather} V(X,0,T) = 0, \end{gather}$$

(2.8)$$\begin{gather}V(X, F(X,T), T) = uF_X, \end{gather}$$

(2.9)$$\begin{gather}p(1,T) = 0, \end{gather}$$

with $Y = F(X,T)$ defining the unknown position of the underbody curve. The requirements (2.7), (2.8) correspond, in turn, to tangential flow at the wall and the kinematic boundary condition at the underbody surface. Condition (2.9) is the Kutta condition ensuring that velocity is finite at the trailing edge. Because (2.3)–(2.5) are hyperbolic, for the Kutta condition to be enforced, a relatively short Euler region must be present ahead of the leading edge in which pressure and velocity change dramatically in anticipation of the flow past the body (Jones & Smith 2003; Smith et al. 2003). The flow is quasi-steady in this region, so Bernoulli's theorem holds here as well as in the rest of the flow, and we can relate the post-Euler flow to the pre-Euler by

(2.10)\begin{equation} \tfrac{1}{2}u_0(Y(0^-))^2 = \tfrac{1}{2}u(0^+,Y(0^+),T)^2+p(0^+,T) = \tfrac{1}{2}u(X,Y(X),T)^2+p(X,T), \end{equation}

on a streamline given by $Y = Y(X)$. The incoming velocity profile $u = u_0(Y)$ is treated as known.

This system can now be significantly reduced by the assumption that the incoming flow is of the form $u_0(Y) = AY+B$ for some constants $A$ and $B$. This is a good approximation of a typical boundary layer profile in the near-wall region. (The case of a full boundary layer will be discussed later in the paper.) For two-dimensional, quasi-steady, inviscid flow, the vorticity equation reveals that the vorticity, equal at leading order to $-u_Y$ when scaled, is constant on streamlines, and because $u_Y = A$ on all incoming streamlines, we have $u_Y = A$ everywhere. Thus, the horizontal velocity is of the form $u = AY+b(X,T)$, and by (2.10), choosing the wall streamline on which $Y(X) =0$, the unknowns $b$ and $p$ are related by

(2.11)\begin{equation} \tfrac{1}{2}b(X,T)^2+p(X,T) = \tfrac{1}{2}B^2. \end{equation}

Assuming the flow at the trailing edge is fully forward, $b(1,T) = B$, by the Kutta condition (2.9). The conservation of volume flux, $q = \int _0^F u\, {\rm d}Y= ({A}/{2})F^2+bF,$ can be used to solve for $b$ and $p$ in terms of $F$,

(2.12)$$\begin{gather} b(X,T) = \frac{\dfrac{A}{2}(F(1,T)^2-F(X,T)^2)+BF(1,T)}{F(X,T)}, \end{gather}$$

(2.13)$$\begin{gather}p(X,T) = \frac{1}{2}B^2 - \frac{1}{2}\left(\frac{\dfrac{A}{2}(F(1,T)^2-F(X,T)^2)+BF(1,T)}{F(X,T)} \right)^2. \end{gather}$$

Substituting (2.13) into the body motion equations (2.1), (2.2), we arrive at a coupled pair of ODEs, which can be solved numerically for given initial conditions. Numerical solutions reveal a wide range of interesting behaviour, which can be split into the categories of ‘crash’ solutions, in which the body makes contact with (impacts on) the wall in finite time, and solutions in which the body height may tend to infinity or otherwise never approach the wall. In § 3, we focus on the case of a crash, and analyse the asymptotic behaviour of the system as contact with the wall is approached. Other classes of solution are discussed in §§ 4, 5 and 6.

3. ‘Crash’ solutions

The first category of solution we will consider is ‘crash’ solutions, which are solutions of the dynamical system (2.1)–(2.2), (2.13) for which the body collides with the wall in finite time, i.e. there exists a point $X_0$ and time $T_0$ such that $F(X_0,T_0) = 0$. These solutions are presented for an arbitrary incoming profile in Jolley et al. (2021) and agree with those presented here; see Appendix B for completeness. However the method presented in that paper cannot cope with solutions with regions of reversed flow and hence although the method described here is less general, it is more powerful (revealing several new solution categories) and its simplicity makes analysis and physical intuition more straightforward.

For instance, we can more easily study the conditions under which ‘crash’ and ‘fly away’ might initially be set into motion. In general, the initial motion of the body is sensitively shape-dependent, but the simple case of a flat body can provide intuition for the circumstances under which the body will initially lift-off. Consider small perturbations taken about the rest state of a flat body at zero angle of incidence, i.e. $F(X,T) = h_0 +\hat {\epsilon }(h_1(T)+(X-\beta )\theta _1(T)),$ for arbitrary small parameter $\hat {\epsilon }$. Then

(3.1)\begin{equation} p ={-}\hat{\epsilon} B\left(A+ \frac{B}{h_0}\right)(1-X)\theta_1, \end{equation}

and hence by (2.1) and (2.2),

(3.2)$$\begin{gather} Mh_{1TT} ={-}\frac{1}{2}B\left(A+ \frac{B}{h_0}\right)\theta_1, \end{gather}$$

(3.3)$$\begin{gather}I\theta_{1TT} ={-}\frac{1}{2}B\left(A+ \frac{B}{h_0}\right)\left(\frac{1}{3}-\beta\right)\theta_1. \end{gather}$$

Thus, $\theta _1$ is governed by a straightforward linear ordinary differential equation, while $h_1$ can be found easily from $\theta _1$. As $A, B, h_0 > 0$, the stability is determined by the sign of $\beta -1/3$. Solving for $h_1$ and identifying the fastest growing term yields the ‘lift-off’ criteria in the three cases,

(3.4)$$\begin{gather} \theta(0)+\theta_T(0)/\lambda < 0, \quad \text{for } \beta > 1/3, \end{gather}$$

(3.5)$$\begin{gather}\theta_T(0) < 0, \quad \text{for } \beta = 1/3, \end{gather}$$

(3.6)$$\begin{gather}h_T(0)-\frac{I\theta_T(0)}{M(1/3-\beta)} >0, \quad \text{for } \beta < 1/3, \end{gather}$$

where $\lambda = \sqrt {B(A+B/h_0)(\beta -1/3)/2I}$. These conditions and numerical results indicate that the crucial factor in inducing lift-off is a sufficiently negative or decreasing initial angle of the body, and similarly ‘crash’ solutions are strongly associated with initially positive or increasing angle.

3.1. First stage

After an initial descent triggered by the sensitive early-time interaction described above, the asymptotic behaviour of the body in approach to the crash can be determined as follows. Let the contact point between the body and the wall be $X=X_0$, and let the crash occur at time $T_0$. Then $F(X_0,T_0) = 0$ and $F_X(X_0,T_0) = 0$, because the gap width must be a minimum where the (smooth) body touches the wall. Considering $X-X_0$, $T-T_0$ to be small, the underbody position curve can be expanded locally,

(3.7)\begin{equation} F(X,T) = \tfrac{1}{2}F_u''(X_0)(X-X_0)^2+h_1(T)+(X_0-\beta)\theta_1(T), \end{equation}

where $h_1,\theta _1$ are the deviations of $h, \theta$ from their values $h_0 = h(T_0), \theta _0 = \theta (T_0)$ at the time of the crash. Assuming that $h-h_0, \theta -\theta _0$ are $O(T_0-T)^n$ say as the crash is approached, such that $F(X,T)$ is $O(T_0-T)^n, n>0$ in the $X$ range $X = X_0+O(T_0-T)^{n/2}$, we can solve for $n$ using the body motion equations to obtain $n = 4/5$. This value shows very close agreement with numerical results (shown in figure 2) for the behaviour of $h$ and $\theta$ as a crash is approached.

Figure 2. First-stage numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, and sinusoidal shape $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = Y+1$ are shown in blue. Initial conditions are $h(0) = 1, \theta (0) = 0.2, h_T(0) = -0.1, \theta _T(0) = 0.$ Red shows a curve varying as $(T_0-T)^{4/5}$ with the right-hand end-point fixed to match the corresponding $h$ or $\theta$ value there (and similarly in figures 3 and 4).

3.2. Second stage

As the body continues to approach contact, its velocity becomes large because $h_T = O(T_0-T)^{-1/5}$, which gives rise to a second time stage. In this second time stage, the pressure force on the body is dominated by that arising from a small region surrounding the contact point, dubbed the ‘inner’ region, for which $X-X_0 = O(T_0-T)^{2/5}$. Defining the order unity time coordinate $s = (T-T_0)/\tau$ for the relevant time scale $\tau$ which is to be found, we rescale variables according to their first-stage behaviour: $\tau ^{2/5}\xi = X-X_0, \tau ^{4/5}\eta = Y, \tau ^{-4/5}\tilde {u} = u, \tau ^{-2/5}\tilde {v} = V, \tau ^{-8/5}\tilde {p} \,{=}\, p, \tau ^{4/5}\phi (\xi,s) {=}\, F(X,T)$ $\text {and}\ \tau ^{4/5}(\tilde {h},\tilde {\theta }) = (h-h_0,\theta -\theta _0).$ This yields the inner region governing equations:

(3.8)$$\begin{gather} \tilde{u}\tilde{u}_\xi ={-}\tilde{p}_\xi, \end{gather}$$

(3.9)$$\begin{gather}0 ={-}\tilde{p}_\eta, \end{gather}$$

(3.10)$$\begin{gather}\tilde{u}_\xi+\tilde{v}_\eta = 0, \end{gather}$$

and boundary conditions for tangential flow, kinematic requirement and matching respectively

(3.11)$$\begin{gather} \tilde{v} = 0, \quad \text{at } \eta = 0, \end{gather}$$

(3.12)$$\begin{gather}\tilde{v} = \tilde{u}\phi_\xi, \quad \text{at } \eta = \phi(\xi,t), \end{gather}$$

(3.13)$$\begin{gather}\tilde{u},\tilde{v},\tilde{p} \to 0 \quad \text{as } \xi \to \pm \infty. \end{gather}$$

The flow here is effectively vorticity-free because $\tilde {u}_\eta = A\tau ^{8/5}$. One may expect that viscosity could come into play but, in fact, it remains small on this time scale (to be found later) provided that $\varDelta ^2 L^5 \ll \hat {\rho }_F/\hat {\rho }_B$. Hence, this momentum equation can be integrated to yield that $\tilde {u}^2/2+\tilde {p}$ is zero. As before, we can now use conservation of volume flux to obtain

(3.14a,b)\begin{equation} \tilde{u} = \frac{q(s)}{\phi(\xi,s)}, \quad \tilde{p} ={-}\frac{1}{2}\frac{q(s)^2}{\phi(\xi,s)^2}, \end{equation}

for the volume flux $q(s)$ (which is not scaled). By the expansion (3.7), the underbody positioning $\phi (\xi,s)$ is of parabolic shape to leading order in this region,

(3.15)\begin{equation} \phi = \tfrac{1}{2}\alpha\xi^2+\tilde{h}(s)+(X_0-\beta)\tilde{\theta}(s), \end{equation}

where we have defined $\alpha = F_u''(X_0)$, effectively the curvature. The body motion is governed entirely by the $O(\tau ^{-8/5})$ inner pressure and hence the body motion equations can now be integrated directly over the inner region to obtain

(3.16)$$\begin{gather} M\tilde{h}_{ss} ={-}\frac{{\rm \pi} q(s)^2}{2\sqrt{2\alpha}}(\tilde{h}+(X_0-\beta)\tilde{\theta})^{{-}3/2}, \end{gather}$$

(3.17)$$\begin{gather}I\tilde{\theta}_{ss} ={-}\frac{{\rm \pi} q(s)^2}{2\sqrt{2\alpha}}(X_0-\beta)(\tilde{h}+(X_0-\beta)\tilde{\theta})^{{-}3/2}. \end{gather}$$

It is now necessary to consider the outer region, where $X - X_0, Y$ are order unity and the body is effectively in contact with the wall. Then $(h,\theta ) = (h_0,\theta _0)+\tau ^{4/5}(\tilde {h},\tilde {\theta })$, but no further scaling is required as other quantities remain order one. Then the new governing equations are

(3.18)$$\begin{gather} u_X+V_Y=0, \end{gather}$$

(3.19)$$\begin{gather}\epsilon\tau^{{-}1}u_s+uu_X+Vu_Y ={-}p_X, \end{gather}$$

(3.20)$$\begin{gather}0={-}p_Y, \end{gather}$$

where $\epsilon = (L/\varDelta ) (\hat {\rho }_F/\widehat {\rho _B})^{-1/2}$. We have included the unsteady term here in anticipation of it coming into play shortly. The boundary conditions are those of tangential flow, kinematics and the Kutta constraint,

(3.21)$$\begin{gather} V=0,\quad \text{at } Y=0, \end{gather}$$

(3.22)$$\begin{gather}V = \epsilon\tau^{{-}1}F_s+uF_X,\quad \text{at }Y=F(X,s), \end{gather}$$

(3.23)$$\begin{gather}p = 0, \quad \text{at } X=1. \end{gather}$$

The velocity can be expanded $u = AY+b_0(X,s)+O(\tau ^{4/5})$ by (2.12), and then use of the kinematic boundary condition (3.22) algebraically eliminates the vorticity term from the $x$-momentum equation,

(3.24)\begin{equation} b_0b_{0X}+p_X = \epsilon\tau^{{-}1/5}u_s. \end{equation}

This can be integrated as before to yield a Bernoulli equation, but because the right-hand side is much larger in the inner region than in the outer, we must consider separately the cases of integrating across the inner region and that of integration only in the outer region. Choosing $\tau = \epsilon ^{5/7}$, we find that this term makes an order-unity contribution when passing over the inner region, which is related to the time derivative of velocity potential that appears in the unsteady Bernoulli equation (effectively there is a jump in velocity potential at the contact point). For this time scale, we then arrive at

(3.25)\begin{equation} \frac{1}{2}b(X,s)^2+p(X,s) = \begin{cases} \dfrac{1}{2}B^2, & X < X_0, \\ \dfrac{1}{2}B^2-\int_{-\infty}^{\infty} \left(\dfrac{q(s)}{\phi(\xi,s)}\right)_s \, {\rm d}\xi, & X > X_0,\end{cases} \end{equation}

as the ‘modified’ Bernoulli equation. Having determined $\tau$, it is clear that the flow is again quasi-steady, so as in the first stage, the volume flux conservation can be exploited to express $b_0$ in terms of the underbody curve $F_0(X) = F_u(X)+h_0+(X-\beta )\theta _0$ analogously to (2.12). Choosing $X=1$ in (3.25) to eliminate the pressure, we arrive at a first-order IDE for $q$,

(3.26)\begin{equation} \frac{1}{2}\left(\frac{q(s)}{F_0(1)}-\frac{A}{2}F_0(1)\right)^2 = \frac{1}{2}B^2-\int_{-\infty}^\infty \left(\frac{q(s)}{\phi(\xi,s)}\right)_s\, {\rm d}\xi. \end{equation}

This can be solved numerically by matching $q$ with the flow from the first stage to obtain an initial condition at large negative $s$, thereby solving for velocity and pressure in both the inner and outer regions. So at each time step, having computed $q$ via (3.26), the body motion equations can be numerically integrated to find $\tilde {h}$ and $\tilde {\theta }$. See figure 3 for results. This again agrees closely with results shown in Jolley et al. (2021) for $u_0(Y) = 2-e^{-Y}$.

Figure 3. Second-stage numerical solutions of $\tilde {h}$ and $\tilde {\theta }$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = Y+1$ are shown in blue. For initial conditions, we take $\tilde {h}, \tilde {\theta }$ large (10 and $-$10, respectively) and $\tilde {h}_s, \tilde {\theta }_s$ small to match with the first stage. Red shows curves varying as $(|s|)^{4/5}$ with the start point fixed to match the initial conditions and yellow shows a curve varying as $O(|s|)$ with the end-point fixed similarly.

The leading-order behaviour of this system on approach to a crash can be determined analytically in similar style to the first stage. Assume as $s\to 0$, where we define the time when the body hits the wall to be zero, that $\tilde {h}, \tilde {\theta } = O(s^N)$ and $k = O(s^M)$. The integral term in (3.26) cannot become large because the left-hand side is positive. Its leading-order term is $O(s^{M-N/2-1})$, which we expect to be large. Hence, it is possible to determine $M = N/2$ by equating this leading-order term of the integral to zero. By consideration of the body motion equations, we come to the conclusion that $N=1$ and $M=1/2$, which agrees well with the numerical results (shown in figure 3). The body now crashes into the wall with finite velocity, and so the solution is complete. This value of $N$ is in agreement with that found in Smith & Wilson (2011).

4. ‘Fly away’ solutions

Several other interesting solutions to the fluid–body interaction (2.3)–(2.10) exist for which the body never collides with the wall. This section will focus on ‘fly away’, where the body tends to depart far from the wall, that is, $h\to \infty$ as $T\to \infty$, starting with $h = O(1)$ at $T=0$. A previous work by Balta & Smith (2018) studied this feature in a similar case with zero vorticity, and solution schemes described in previous sections reproduce their findings that $h \to \infty, \theta \to -\infty$ parabolically for $T \to \infty$ for that zero vorticity case. However, for non-zero vorticity, fly away solutions are completely distinct and, in fact, $\theta$ oscillates and remains bounded. Some understanding of the cause of this distinction can be gained by inspecting the formulae (2.12) and (2.13) for velocity $b$ and pressure $p$. With zero vorticity, as the body tilts, the pressure gradient steepens, generally causing it to tilt further. However, with non-zero vorticity $A$, it is clear that as $\theta \to -\infty$, eventually there will be a region of reverse flow at the leading edge and correspondingly pressure reaches its maximum at the zero of $b$. Once the pressure reaches its maximum, the continued steepening of the pressure gradient causes a shift in distribution eventually altering the direction of rotation, causing the qualitative difference we see between these solutions and those of Balta & Smith (2018). Because the arbitrary velocity profile scheme discussed above does not apply to reverse flow, we return to the constant vorticity case for which our results (2.12) and (2.13) describing velocity and pressure in terms of underbody position still hold (flow is still fully forward in the pre-Euler region and at the trailing edge, so this analysis still applies).

To consider behaviour following on from such initial lift-off over long times, let time $T = \tau \tilde {T}$ with $\tau \gg 1$ and body height $h = \tau ^2H_2(\tilde {T})+\tau H_1(\tilde {T})+H_0(T)$, where the final term varies on the short time scale, while the first two vary on long time only. We also let $\theta = \theta _0(T)+O(\tau ^{-1})$ and seek to study the leading-order behaviour of the height $h$. By (2.12), the $x$-dependent part of the velocity is then

(4.1)\begin{equation} b = A\left[\theta_0 (1-X)-F_u(X)\right]+B+O(\tau^{{-}1}), \end{equation}

while from (2.13), the pressure is

(4.2)\begin{equation} p ={-}A\theta_0(1-X)(B-AF_u(X))-\tfrac{1}{2}A^2\theta_0^2(1-X)^2+ABF_u(X) -\tfrac{1}{2}AF_u(X)^2+O(\tau^{{-}1}). \end{equation}

The body motion equations are then

(4.3)$$\begin{gather} M\left(\frac{{\rm d}^2H_2}{{\rm d}\tilde{T}^2}+\frac{{\rm d}^2H_0}{{\rm d}T^2}\right) = a_0+a_1\theta_0+a_2\theta_0^2, \end{gather}$$

(4.4)$$\begin{gather}I\frac{{\rm d}^2\theta_0}{{\rm d}T^2} = c_0+c_1\theta_0+c_2\theta_0^2, \end{gather}$$

where

(4.5)\begin{equation} \left. \begin{aligned} a_0 & = \int_0^1 ABF_u(X) -\frac{1}{2}A^2F_u(X)^2\, {\rm d}X, \\ a_1 & ={-}\frac{1}{2}AB+A^2\int_0^1(1-X)F_u(X)\, {\rm d}X,\\ a_2 & ={-}\frac{1}{6}A^2, c_0 = \int_0^1 (ABF_u(X) -\frac{1}{2}A^2F_u(X)^2)(X-\beta) \, {\rm d}X, \\ c_1 & = AB\left(-\frac{1}{6}+\frac{\beta}{2}\right)+A^2\int_0^1 F_u(X)(1-X)(X-\beta)\, {\rm d}X \quad {\rm and}\\ c_2 & = \frac{1}{2}A^2\left(-\frac{1}{12}+\frac{\beta}{3}\right) \end{aligned} \right\} \end{equation}

are constants. Solving (4.4) shows that $\theta _0$ is a periodic function, oscillating on the fast time scale, and thus has a constant mean value $\langle \theta _0\rangle$ over a period in $T$. Then, after taking its mean, (4.3) can be straightforwardly integrated to obtain

(4.6)\begin{equation} H_2 = \frac{1}{2M}\left( a_0+a_1\langle\theta_0\rangle+a_2\langle\theta_0^2\rangle\right)\tilde{T}^2. \end{equation}

Provided the coefficient of $\tilde {T}^2$ here is positive, this is a fly away solution. Figure 4 shows that the solution from the full computation and the approximate parabolic curve from (4.6) overlap very closely, which indicates firm agreement; $\theta (T)$ is also shown to be an oscillating function of time as predicted analytically.

Figure 4. A fly away case. The numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = Y+1$ are shown in blue. Initial conditions are $h(0) = 1, \theta (0) = -0.2, h_T(0) = 0, \theta _T(0) = 0.$ Red shows the parabolic approximation given in (4.6) (it overlaps very closely with the blue curve).

4.1. The $\theta _0$ limit cycle

Of natural interest is the question of how (or if) one can tell if fly away will occur without running full simulations. In light of (4.6), it is clear that to do this requires an understanding of the $\theta _0$ limit cycle and its mean and variance. Equation (4.4) can be integrated to obtain

(4.7)\begin{equation} I\left(\frac{{\rm d}\theta_0}{{\rm d}T}\right)^2 = 2c_0\theta_0+c_1\theta_0^2+\frac{2}{3}c_2\theta_0^3 + k, \end{equation}

for some constant $k$, which is related to $\langle \theta _0\rangle$ by

(4.8)\begin{equation} \langle\theta_0\rangle = \frac{1}{T_0}\int_T^{T+T_0} \theta_0\, {\rm d}T = \frac{2\sqrt{I}}{T_0}\int_{\theta_{{min}}}^{\theta_{{max}}} \frac{u}{\sqrt{2c_0u+c_1u^2+\tfrac{2}{3}c_2u^3+k}}\, {\rm d}u, \end{equation}

where the minimum $\theta _{{min}}$ and the maximum $\theta _{{max}}$ of $\theta _0$ are functions of $k$ in the sense that they are the two real solutions to the equation $0= 2c_0\theta _0+c_1\theta _0^2+\frac {2}{3}c_2\theta _0^3 + k.$ Similarly, the period $T_0$ is given by

(4.9)\begin{equation} T_0 = \int_T^{T+T_0}\,{\rm d}T = 2\sqrt{I}\int_{\theta_{{min}}}^{\theta_{{max}}}\frac{1}{\sqrt{2c_0u+c_1u^2+\tfrac{2}{3}c_2u^3+k}}\, {\rm d}u. \end{equation}

It is not easy to determine the value of $k$ as it depends on the behaviour of the system prior to entry into the ‘large $h$’ regime. However, $k$ can be easily bounded to force the cubic in (4.7) to have two real roots (else a limit cycle cannot exist), which in turn can produce a range of viable $\langle \theta _0\rangle$ values. Because $\langle \theta _{0TT}\rangle = 0$, we have

(4.10)\begin{equation} \langle\theta_0^2\rangle ={-}\frac{1}{c_2}\left(c_0+c_1\langle\theta_0\rangle\right), \end{equation}

which can be used to bound the parabola coefficient. Thus, although we cannot determine the exact behaviour of the body for given initial conditions without running full simulations, we will know the range of possible values of the parabola coefficient which gives a range of possible behaviours, e.g. if the range is entirely positive, then fly away is guaranteed following lift-off, and if not, there may exist some other types of solution such as those discussed in §§ 5 and 6 below.

5. ‘Bouncing’ solutions

Numerical results reveal the existence of solutions for which $h$ ‘bounces’, i.e. reaches large values then comes back to $O(1)$ repeatedly. We will show cases where $H_2 \equiv 0$, and hence $a_0+a_1\langle \theta _0\rangle +a_2\langle \theta _0^2\rangle = 0.$ To examine these, we need to take the analysis to the next order, so let $h = \tau H_1(\tilde {T})+H_0(T)$ and $\theta = \theta _0(T)+\tau ^{-1}\theta _{-1}(T).$ Proceeding as before, the body motion equations at the next order are

(5.1)$$\begin{gather} M\frac{{\rm d}^2H_1}{{\rm d}\tilde{T}^2} = \frac{k_1}{H_1}+a_1\langle\theta_{{-}1}\rangle+2a_2\langle\theta_0\theta_{{-}1}\rangle, \end{gather}$$

(5.2)$$\begin{gather}I\frac{{\rm d}^2\theta_{{-}1}}{{\rm d}T^2} = \frac{k_2(T)}{H_1}+ c_1\theta_{{-}1}+2c_2\theta_0\theta_{{-}1}, \end{gather}$$

where

(5.3)$$\begin{gather} k_1 ={-}\int_0^1 \left\langle\vphantom{\frac{A}{2}}\left(A\left[(1-X)\theta_0-F_u(X)\right]+B\right)\left[(1-X)\theta_0-F_u(X)\right]\right.\nonumber\\ \left.\times\left(\frac{A}{2}\left[(1-X)\theta_0-F_u(X)\right]+B\right)\right\rangle\, {\rm d}X, \end{gather}$$

and

(5.4)$$\begin{gather} k_2(T)={-}\int_0^1 \left(A\left[(1-X)\theta_0-F_u(X)\right]+B\right)\left[(1-X)\theta_0-F_u(X)\right]\nonumber\\ \times\left(\frac{A}{2}\left[(1-X)\theta_0-F_u(X)\right]+B\right)(X-\beta)\, {\rm d}X, \end{gather}$$

define $k_1, k_2$ as the coefficients above.

5.1. Analysis for small $k_1$, $k_2$

We will assume that $k_1$ and $k_2$ are small: this is partly to make progress and partly because this has tended to be the case numerically. The possibility that $k_{1,2}$ are not small will be addressed at the end of the subsection. Then, for $H_1 = O(1),$ we have the system

(5.5)$$\begin{gather} M\frac{{\rm d}^2H_1}{{\rm d}\tilde{T}^2} = a_1\langle\theta_{{-}1}\rangle+2a_2\langle\theta_0\theta_{{-}1}\rangle, \end{gather}$$

(5.6)$$\begin{gather}I\frac{{\rm d}^2\theta_{{-}1}}{{\rm d}T^2} = c_1\theta_{{-}1}+2c_2\theta_0\theta_1. \end{gather}$$

Similarly to the previous solution, $\theta _{-1}$ is a periodic function and $\langle \theta _{-1}\rangle$ is a constant, which indicates that $H_1$ is parabolic to leading order. If the right-hand side of (5.5) is positive, then this describes the solution for all time, and we see a weaker but similar fly away solution. If it is negative, however, then $H_1 \to 0$ at some time $\tilde {T} = \tilde {T}_1$, and there is a second solution when $H_1$ is small. In this case, we rescale for small $H_1$ by letting $H_1 = k_1^{3/2}y$, $\theta _{-1} = k_1^{-1/2}\varTheta$ and $\tilde {T} = \tilde {T}_1+k_1t$ to find the new system

(5.7)$$\begin{gather} M\frac{{\rm d}^2y}{{\rm d}t^2} = \frac{1}{y}+k_1^{1/2}C(t), \end{gather}$$

(5.8)$$\begin{gather}I\frac{{\rm d}^2\varTheta}{{\rm d}T^2} = \frac{k_2(T)}{k_1}\frac{1}{y}+c_1\varTheta+2c_2\theta_0\varTheta, \end{gather}$$

where we have also defined $C$ as the right-hand side of (5.5) for convenience ($C = C(t)$ is not a constant in this boundary layer solution). The expansion $y = y_0+k_1^{1/2}y_1+\cdots$ yields

(5.9)$$\begin{gather} M\frac{{\rm d}^2y_0}{{\rm d}t^2} = \frac{1}{y_0}, \end{gather}$$

(5.10)$$\begin{gather}M\frac{{\rm d}^2y_1}{{\rm d}t^2} ={-}\frac{y_1}{y_0^2} + C(t). \end{gather}$$

The solution to (5.9) is

(5.11)\begin{equation} y_0 = y_{0\,min}\exp\left(\left[\mathrm{erfi}^{{-}1}\left(\sqrt{\frac{2}{{\rm \pi} M}}\frac{t}{y_{0\,min}}\right)\right]^2\right), \end{equation}

where $\mathrm {erfi}^{-1}$ is the inverse of the imaginary error function. We can understand its behaviour by integrating and considering its phase portrait,

(5.12)\begin{equation} \frac{M}{2}\left(\frac{{\rm d}y_0}{{\rm d}t}\right)^2 = \log |y_0| + \text{const.} \end{equation}

Here, starting from $t = -\infty$, $y_0$ decreases from $\infty$ close to linearly (like $t(\log t)^{1/2}$), passes through a minimum point at $(0,y_{0\,min})$, then leaves symmetrically to $+\infty$, so that $y_0$ and hence $H_1$ apparently ‘bounces’. Numerical results from full simulations (shown in figure 5) confirm this behaviour, and simulations of the reduced system (5.7) and (5.8) in figure 6 show that the effect of the bounce is ultimately to push $\theta _{-1}$ into a different limit cycle, so that in the next excursion, although $C$ is again constant, it takes a new value. This produces repeated excursions of varying height and length, as observed numerically.

Figure 5. A bouncing case. The numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = 3.2Y+1$ are shown in blue. Initial conditions are $h(0) = 1, \theta (0) = -0.2, h_T(0) = 0, \theta _T(0) = 0.$

Figure 6. Bouncing. The results of a simulation for which $\theta _0$ is calculated using (4.4), and $h$ and $\theta _{-1}$ are calculated using (5.5) and (5.6). Values for $k_{1,2}$ and $C$ as well as initial conditions for $\theta _0$ are taken as constants from full simulations. Because $C$ is taken as constant, we see $h$ repeats the same excursion – but in reality, the shift in the $\theta _{-1}$ limit cycle would produce distinct excursions as shown in the full results of figure 5.

The ‘bouncing’ behaviour here should be observed regardless of the sizes of $k_1, k_2$, provided they are similar in size; if $k_{1,2}/H_1$ is small, we see $\theta _{-1}$ oscillating and $H_1$ decreasing parabolically, until $H_1$ becomes small enough that $k_1/H$ becomes the dominant term in (5.5), causing the body to ‘bounce’ and re-enter the parabolic regime.

6. Other behaviour

The findings above point to further behaviours being possible in the current fluid/body interactions. These behaviours are associated with interactions in which the scaled body height $h$ oscillates as discussed below in § 6.1, or an alternative crash effect occurs which is described in § 6.2 or the scaled angle $\theta$ ceases to remain bounded as investigated in § 6.3 below.

6.1. Oscillating $h$ solutions

These are solutions for which $H_2 \equiv H_1 \equiv 0$. In this case, $h = H_0(T)$ and $\theta = \theta _0(T)$ are both periodic functions to leading order and $h$ never becomes large nor crashes. This is shown in figure 7. Because for these solutions $h$ and $\theta$ are of comparable size, analysis is difficult and less fruitful.

Figure 7. An ‘oscillating $h$’ solution. The numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = 4Y+1$ are shown in blue. Initial conditions are $h(0) = 1, \theta (0) = -0.2, h_T(0) = 0, \theta _T(0) = 0.$

6.2. Parabolic crash

If $h$ is large in the initial conditions, strong negative parabolas can be realised. The body then crashes with $\langle h\rangle$ approximately linear, and $h, \theta$ oscillating. This is shown in figure 8. On smaller scales such that the effect of $H_0$ is visible, the 4/5 behaviour of § 3 is still shown, but on large scales, it cannot be seen in comparison to the large linear term in $h$.

Figure 8. An ‘alternative crash’. The numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = 5Y+1$ are shown in blue. Initial conditions are $h(0) = 100, \theta (0) = -0.2, h_T(0) = 0, \theta _T(0) = 0.$

6.3. The $\theta$ escape solutions

Using (4.4) to examine the phase portrait for $\theta _0$, it is clear that there should also exist some ‘escape’ solutions, in which $\theta _0 \to \infty$. These can be realised if starting with large $h$. This also eventually causes a crash as the body turns on its side, but it might be physically questionable as large values of $\theta$ are not necessarily compatible with the small angle approximation that is used in the definition of the body curve. A numerical simulation of this case is shown in figure 9.

Figure 9. A ‘$\theta$-escape’ solution. The numerical solutions of $h$ and $\theta$ for a body with $M = 1, I=0.2$, $F_u(X) = -0.2\sin {\rm \pi}X$ and incoming profile $u_0 = Y+1$ are shown in blue. Initial conditions are $h(0) = 100, \theta (0) = 0.2, h_T(0) = 0, \theta _T(0) = 0.$

7. Discussion and conclusions

The work has focussed on dynamic fluid–body interaction near a wall with substantial incoming flow vorticity, as is appropriate for the boundary layer setting. Modelling assumptions that the body is thin and in the appropriate range of lengths, as well as heavy, allowed the reduction of the problem to a coupled pair of nonlinear IDEs. Numerical solutions confirmed and explained by accompanying analysis have revealed that the possible behaviours of the body are wide-ranging and heavily shape-dependent, and in particular, the solutions can be split into two categories. The first class of solutions is so-called crashes, for which the asymptotic dependency of the body height and angle on time has been determined in the limit of a crash, through two time stages. The final result that the body velocity is constant at the onset of a crash agrees with the results of the zero-vorticity study by Smith & Wilson (2011), as would be expected because the effect of vorticity on the leading order pressure contribution in this limit is negligible. The problem here differs from their work also in the assumption that the body is heavy which allows for modelling the flow as quasi-steady and, in turn, causes the need for two time stages. The crash problem has also been treated for arbitrary incoming velocity profile, with analysis and numerics indicating that any velocity profile, provided it is fully forward, yields the same local body motion asymptotically for a crash scenario, as we might expect given the underbody comes close to the wall in this case.

Several reasons for concentrating mostly on the constant vorticity case here are notable. Physically the fluid–body interaction for a body outside the boundary layer or in the outer reaches of a boundary layer are covered by Balta & Smith (2018). The interaction for constant vorticity describes the opposite regime where the body lies in the depths of the boundary layer, closer to the wall. Mathematically, the case of constant vorticity is simpler and clearer than the full case. We tackled it in this work for clash phenomena and then readily extended that to the case of a full boundary layer in an appendix. The constant-vorticity case also readily allows interesting new phenomena, namely fly away and bouncing as well as other phenomena, to be identified subsequently.

In contrast to the crash solutions, the present range of interesting non-crash behaviours involving oscillation of the body angle is unique to flows with vorticity and thus in stark contrast to Balta & Smith (2018). The presence of vorticity necessitates that as the body angle decreases, the pressure will eventually reach its maximum point and eventually cause the body to tilt the other way, causing an oscillating body angle. (Without vorticity, pressure cannot reach this point and so this oscillation then never happens.) Not only does this oscillation give rise to a wider range of behaviours, but it is perhaps more physically realistic because it does not result in nominally infinite values of the scaled body angle. A natural question for further work would be whether, as for crashes, the same qualitative behaviour is observed regardless of the velocity profile. It is already clear that the presence of vorticity changes the body behaviour substantially in these cases, and in reality, one would expect vorticity to eventually decrease as the body rises sufficiently far from the wall. To numerically solve for fly away cases in an arbitrary velocity profile would require adaptation of the modelling in Appendix B to accommodate reverse flow, i.e. potential changes in the sign of the square root in (B1). Smith & Palmer (2019) consider reverse flow and separation in a related setting.

On physical parameter values, in realistic ranges, the typical ice particles are of length $10^{-3}$ m to $10^{-5}$ m and the sizes of the non-dimensional numbers involved are $10^6$ for the Froude number and $10^4$ to $10^6$ for the global Reynolds number, with the representative particle Reynolds number being $10^2$ to $10^3$. Thus, concerning the Froude number, the typical body weight is very small compared with the pressure forces. The weight would be expected to have only a longer-term influence on the present solution behaviours. Further details on icing conditions and the ranges of parameters are given in Norde (2017), Palmer & Smith (2019), while, concerning the Reynolds number, viscous effects are considered in Appendix A.

There are no existing experiments, as far as we know, on bouncing or fly away phenomena, for example, or on crash responses in detail. Both categories of solution have important industrial applications in the field of aircraft icing. Understanding what kind of bodies under what circumstances are likely to collide with an airplane wing is crucial for the development of icing protection systems. Understanding the position and behaviour of those that do not crash is important as their presence may influence the behaviour of other bodies in multi-particle simulations, as well as the use of knowing which bodies do not crash for ice protection. In a similar vein, the behaviour of the body following a crash would be an interesting question for further work, as well as conducting multi-particle studies as in Smith & Ellis (2010) and Smith et al. (2019). Other extensions to the analysis not yet mentioned include accounting for three-dimensional bodies and flow and the inclusion of thermodynamic properties (melting/freezing) of ice.