3.1.1. Modelling assumptions and simulation set-up

The gas surrounding the cylinder is modelled as a perfect gas with heat capacity ratio $\gamma = 1.4$, so that the compressible, unsteady Navier–Stokes equations govern the behaviour of the flow. These were solved in their Reynolds-averaged form (URANS) using the solver rhoCentralFoam (Greenshields et al. Reference Greenshields, Weller, Gasparini and Reese2010), which is part of the Open Source CFD software package OpenFoam (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998), version 5.x. The choice of this numerical approach is driven by its simplicity and the ready availability of open-source code. The viscosity of the gas was determined using Sutherland's law, with coefficients changed from case to case to achieve flow situations of different Reynolds numbers, while keeping the cylinder diameter D constant for all simulations to facilitate mesh generation. The $k - \omega - \textrm{SST}$ turbulence model (Menter Reference Menter1994) was used as closure for the URANS equations in all conducted simulations in this study. This model has been used by other authors to investigate similar flow scenarios (Catalano & Amato Reference Catalano and Amato2003; Benim, Pasqualotto & Suh Reference Benim, Pasqualotto and Suh2008; Rosetti, Vaz & Fujarra Reference Rosetti, Vaz and Fujarra2012; Stringer, Zang & Hillis Reference Stringer, Zang and Hillis2014). While the inherent drawbacks of URANS modelling become apparent mostly in the critical flow regime (Stringer et al. Reference Stringer, Zang and Hillis2014), the $k - \omega - \textrm{SST}$ model was found to outperform other two-equation models for this type of flow situation (Catalano & Amato Reference Catalano and Amato2003; Benim et al. Reference Benim, Pasqualotto and Suh2008).

The boundary conditions were assigned to be of the symmetry type for the top and bottom boundaries; zero-gradient boundaries were assigned at the left and right end of the domain and a no-slip condition was enforced on the velocity field on the cylinder surface. Non-physical wave reflections from the boundaries were precluded by choosing a sufficiently large domain size. In order to decrease the computational effort, the assumption of two-dimensional flow was employed. Whereas the assumption of two-dimensionality is accurate for the initial wave reflection and diffraction (Sun et al. Reference Sun, Saito, Takayama and Tanno2005), resolving the flow structures in the wake of an object would necessitate a three-dimensional approach. The computational effort to conduct a three-dimensional parametric study was, however, deemed prohibitively expensive, in consideration of the objectives of the study and of the fact that other researchers have found reasonable agreement, in terms of drag, comparing two-dimensional URANS simulations to experiments, three-dimensional URANS and large eddy simulations for high-Reynolds-number flows past bluff bodies (Rodi Reference Rodi1997; Iaccarino et al. Reference Iaccarino, Ooi, Durbin and Behnia2003; Meliga, Pujals & Serre Reference Meliga, Pujals and Serre2012; Stringer et al. Reference Stringer, Zang and Hillis2014). It is clear, however, that limitations arise from employing a two-dimensional URANS approach. The fidelity of predicting turbulent transition, flow detachment and vortices is not comparable to more detailed approaches, such as large eddy simulations or direct numerical simulations. The implications on the predictions will be discussed later.

The pressure wave or shock wave was modelled as an initial field of pressure, particle velocity and temperature immediately adjacent to the object. The equations defining the spatial distribution of these quantities as a function of the wave parameters ${p_i},\;{\alpha _r},\;{\tau _i}$ are given in appendix A. The rest of the fluid domain was assigned homogeneous initial conditions of $p = {p_0}$, $T = {T_0}$, $v = 0\,\textrm{m}\,{\textrm{s}^{ - 1}}$. Figure 3 shows the computational domain and the initial pressure contours for the case ${\alpha _r} = 0.5$, ${\tau _i} = 10$, ${p_i}/{p_0} = 1$. It can be seen that the front of the pressure wave is initially placed almost in contact with the cylinder, such that the wave distortion before arrival at the object is minimised.

Figure 3.

Initial field of pressure for the case ${\alpha _r} = 0.5,\;{\tau _i} = 10,\;{p_i}/{p_0} = 1$.

The chosen integration schemes were of first order in time and second order in space. Due to the explicit prediction of the fluxes in rhoCentralFoam (Greenshields et al. Reference Greenshields, Weller, Gasparini and Reese2010), the maximum time step was determined by enforcing a Courant–Friedrichs–Lewy number less than 0.2. Interpolation of the convective terms was accomplished with the scheme by Kurganov & Tadmor (Reference Kurganov and Tadmor2000), employing flux limiters after van Leer (Reference van Leer1974), as recommended by Greenshields et al. (Reference Greenshields, Weller, Gasparini and Reese2010). Zółtak & Drikakis (Reference Zółtak and Drikakis1998) have compared various computational schemes as well as static and adaptive meshing techniques for the simulation of the interaction of a shock wave with a cylinder. It was concluded that, whereas small differences between the results obtained with different computational schemes exist, very good agreement between static and adaptive mesh techniques was observed.

Due to the wide range of cases to be examined and the analysis of non-shock waves with the same computational scheme, a static meshing approach is employed here. In order to efficiently simulate cases for a wide range of wavelengths, i.e. ${\tau _i} = [1,200]$, the domain size and cell distribution need to adapt, however, to the individual cases. Due to the high local gradients of the flow, the mesh was successively refined towards the region closest to the cylinder to a side length ${\rm \Delta} x_{{cyl}}, \ {\rm with} \ D/\mathrm{\Delta }{x_{cyl}} = 200$. To fully resolve the boundary layer, the mesh was further refined in the direction normal to the cylinder surface to guarantee that for the non-dimensional wall distance y ⁺ holds $ {y^ + } \lt 1$ in the first cell, in line with recommendations by Menter (Reference Menter1994) and findings by Benim et al. (Reference Benim, Pasqualotto and Suh2008) for turbulent flow past circular cylinders. In order to ensure sufficient resolution of the wave, the maximum cell size in the whole domain was limited to

(3.1)\begin{equation}\mathrm{\Delta }{x_{max}} = \frac{{{\lambda _i}}}{{200}}.\end{equation}

In figure 4, an example of the mesh for the case ${\alpha _r} = 0.5$, ${\tau _i} = 10$, ${p_i}/{p_0} = 1$ is depicted. The successive refinement of the cells towards the object surface leads to a sufficient resolution of the zones with the highest gradients due to wave diffraction, the influence of the boundary layer, flow separation and vortex shedding. We note that a coarsened mesh is depicted in figure 4 for visualisation purposes.

Figure 4.

Detail of the mesh for cases with ${\tau _i} = 10$, displayed with fourfold coarsened grid.

3.1.2. Mesh convergence study

An extensive mesh convergence study was conducted to estimate the spatial and temporal discretisation errors of the CFD simulations. As wide ranges of the non-dimensional parameters ${p_i}/{p_0},\;{\alpha _r},\;{\tau _i},\;R{e_i}$ are of interest in this study, numerous cases with different parameter combinations needed to be examined. The following parameter combinations were considered:

(3.2)\begin{equation}\left. {\begin{array}{c@{}} {{\alpha_r} = \{ 0,0.5\} ,\quad {p_i}/{p_0} = \{ 0.1,1,3\} ,}\\ {{\tau_i} = \{ 1,10,50\} ,\quad Re_{i} = \{ {{10}^2},{{10}^4},{{10}^6}\} ,} \end{array}} \right\}\end{equation}

where the highest pressure ratio was omitted for the finite rise time case, ${\alpha _r} = 0.5$, as pressure waves of this amplitude turn into shock waves almost immediately. Therefore, spatial and temporal convergence were investigated for a total of 45 cases following the widely used methodology by Roache (Reference Roache1994). We quantify convergence using the grid convergence index (GCI) of the form (Roache Reference Roache1994)

(3.3)\begin{equation}\textrm{GC}{\textrm{I}_{12}} = \frac{{{F_s}}}{{{r^{\hat{p}}} - 1}}\left|{\frac{{{f_2} - {f_1}}}{{{f_1}}}} \right|,\end{equation}

where ${F_s}$ denotes a factor of safety, $\hat{p}$ the observed order of convergence, ${f_1},\;{f_2}$ denote scalar solution values obtained on the finest and second finest grid and r is the refinement factor. The observed order of convergence can be computed as

(3.4)\begin{equation}\hat{p} = \frac{{\ln \left( {\dfrac{{{f_3} - {f_2}}}{{{\,f_2} - {f_1}}}} \right)}}{{\ln (r)}},\end{equation}

making use of a solution value obtained on a third grid f ₃. Typically, values for $\hat{p}$ between one and two were obtained, which is to be expected as the spatial discretisation scheme reduces to first order in the vicinity of shocks (Banks, Aslam & Rider Reference Banks, Aslam and Rider2008). Three meshes were used for each case, with an isotropic refinement factor of $r = 1.5$. Applying the recommendations proposed by Roy (Reference Roy2010) for the factor of safety and the limits of $\hat{p}$, we obtained the maximum and mean GCI values across all investigated cases listed in table 1. The chosen solution variables were the maximum drag force on the cylinder, ${F_{max}}$, the maximum imparted impulse on the cylinder, ${I_{max}}$, and the time at which the maximum drag value is reached, ${t_{max}}$. It can be seen that the maximum GCI values in both force and time are found to be around 10 %; these were found for the cases employing the lowest pressure ratio of ${p_i}/{p_0} = 0.1$, whereas the maximum GCI value in terms of maximum impulse was slightly lower. The mean values across all 45 cases were significantly smaller and were deemed satisfactory.

Table 1.

Results of the mesh convergence study in terms of the GCI.

3.1.3. Validation of numerical model with results from literature

We now compare results obtained with the proposed numerical model to results available in literature (see figure 2). As can be seen in figure 5, the results of the present numerical model are in very good agreement with those obtained by Drikakis et al. (Reference Drikakis, Ofengeim, Timofeev and Voionovich1997). Experimentally obtained results published by Takayama & Itoh (Reference Takayama and Itoh1985) are also included, which are in broad agreement with the two sets of numerical results.

Figure 5.

Comparison of simulation results with numerical model to results from literature for an example case with ${p_i}/{p_0} = 0.805$, ${\alpha _r} = 0$, $R{e_i} = 7 \times {10^5}$, ${\tau _i} \to \infty$.

By means of figure 5 we introduce a scaling of the force on the cylinder different to most studies in literature. Whereas most existing studies use the classic drag coefficient of the form $F/(0.5{\rho _i}v_i^2D)$ (left ordinate), we choose to employ a scaling of the form $F/({p_i}D)$ (right ordinate). This is due to the fact that drag coefficients for pressure waves with small overpressures $({p_i} \to 0)$ tend to infinity, whereas $F/({p_i}D)$ allows compact representation of results for wide ranges of overpressure.

3.1.4. Parametric study

An extensive parametric study was conducted exploring the following parameter combinations:

(3.5)\begin{equation}\left. {\begin{array}{c@{}} {R{e_i} = \{ {{10}^2},{{10}^4},{{10}^6}\} ,\quad {p_i}/{p_0} = \{ 0.1,0.5,1,1.5,2,3\} ,}\\ {{\tau_i} = \{ 1,5,10,20,30,40,50,60,100,200\} ,\quad {\alpha_r} = \{ 0,0.25,0.5\} .} \end{array}} \right\}\end{equation}

It was deemed unrealistic to encounter finite rise time pressure waves of amplitude ${p_i}/{p_0} \gt 2$ in practice, as these evolve rapidly into shock waves. The highest pressure ratio was, therefore, only used for ${\alpha _r} = 0$, whereas the two longest wavelengths were only combined with ${\alpha _r} = 0.5$, leading to a total of 414 cases. The simulations were run on a high-performance computing cluster, using 32 processors. A significant amount of explicit time steps was necessary due to the locally refined mesh and long simulation times, leading to run times for the individual cases of up to six days.

3.2.1. Propagation of a finite amplitude wave

The analytically formulated models described in the following make use of time-dependent values of the flow variables pressure, particle velocity, density, Reynolds number and Mach number, which is denoted as M. As defined in § 2, we assume an incoming wave of triangular overpressure evolution. At $t = 0$, the front of the wave has reached the front surface of the structure. Subsequently, the wave propagates along the object, inducing diffraction and a transient flow field. As the pressure waves of interest in this study are of significant amplitude, the wave shape distorts during propagation due to the local differences in speed of sound and particle velocity (see e.g. Liepmann Reference Liepmann1957). This effect is further intensified by the wave diffraction and reflection, which increase the differences in the local properties across the wave.

A semianalytical approach, based on the method of characteristics, is used to compute the time-dependent flow variables around the structure: the given wave is split into 100 individual wavelets, each possessing an individual particle velocity and local speed of sound (see appendix A). After a time step the jth wavelet has advanced by the distance $\Delta {x_j} = ({c_j} + {v_j})\Delta t$, leading to a distortion of the initial wave shape. The cases of a shock wave and of a pressure wave that develops into a shock wave need special treatment, as the velocity of a shock front is not equal to that of a simple (non-shock) wave of the same amplitude. In first-order approximation, the shock front propagates at the mean value of the velocities of the simple waves in front of and behind the shock front (Courant Reference Courant1948). The wavelets behind the shock front, therefore, propagate faster and, by catching up with the wave front, continuously change the pressure and velocity of the shock front.

The procedure is illustrated in figure 7, which shows the distortion of a finite amplitude wave, the transition to a shock wave once the wave front is overtaken by the wavelets behind it, and the decay of the maximum pressure after the peak of the wave has overtaken the front. At every point in space or time, the arrival of the individual wavelets yields a discrete distribution of pressure, which was interpolated linearly to approximate the distorted wave shape. This simple procedure yields time-dependent flow properties of good accuracy at different positions (e.g. front $(x = 0)$ and back edge $(x = D)$) around the cylinder. Average or ‘effective’ quantities are then calculated by averaging the flow variables over a number of points on the cylinder surface with uniform angular spacing of 5°. These are denoted in the following with an overbar and the index ‘cyl’.

Figure 7.

Prediction of the distortion of a finite amplitude wave via the method of characteristics.

3.2.2. Diffraction model

When a pressure wave encounters an impermeable solid object, the wave is subject to reflection and diffraction which cause a highly transient pressure distribution on the object surface. With the assumption of small disturbances, the well-known linearised equations of motion of acoustics can be deduced. The procedure outlined in Friedman & Shaw (Reference Friedman and Shaw1960) and Shaw (Reference Shaw1975) is followed, but amended by the introduction of the reflection coefficient

(3.7)\begin{equation}{C_R}(t) = \frac{{({3\gamma - 1} )\left( {\dfrac{{{p_{in}}(\theta ,t)}}{{{p_0}}}} \right) + 4\gamma }}{{({\gamma - 1} )\left( {\dfrac{{{p_{in}}(\theta ,t)}}{{{p_0}}}} \right) + 2\gamma }},\end{equation}

which permits an approximation of the effect of wave reflection from the cylinder surface under the influence of compressibility. The definition of the reflection coefficient in (3.7) corresponds to the normal reflection of a shock wave from a rigid surface (see e.g. Courant Reference Courant1948). While the reflection coefficient for an isentropic wave $({\alpha _r} \gt 0)$ is higher for large pressure ratios (Gauch, Montomoli & Tagarielli Reference Gauch, Montomoli and Tagarielli2018), this difference is negligible for the isentropic waves treated in this study $({p_i}/{p_0} \le 2)$. In (3.7) ${p_{in}}(\theta ,t)$ denotes the incident pressure at time t and angular position $\theta $, with $\theta = 0$ at the point that the wave encounters first. The time retarded integral equation for the pressure on the cylinder surface

(3.8)\begin{align} p(r,t) & = {p_r} + \dfrac{1}{{2{\rm \pi} }}\int_S {{\left\{ {\left( {\dfrac{p}{{{R^2}}} + \dfrac{1}{{{c_0}R}}\dfrac{{\partial p}}{{\partial {t_0}}}} \right)\dfrac{{\partial R}}{{\partial {n_0}}}} \right\}}_{{t_0} = t - R/{c_0}}} \ \textrm{d}{S_0} \nonumber\\ & \quad + \dfrac{1}{{\rm \pi} }\int_0^{2{\rm \pi} } {\dfrac{{\partial {z_{ou}}}}{{\partial t}}{{\left[ {\dfrac{p}{{{c_0}R}}\dfrac{{\partial R}}{{\partial {n_0}}}} \right]}_{t= {t_0} - R/{c_0}\atop {z_0} = {z_{ou}}}}\textrm{d}\theta } \end{align}

is solved at points with equal angular spacing of 10° over the cylinder surface at every time step. The time step has to be chosen such that points influence neighbouring points only with their past values, i.e. $\Delta t \lt \Delta \theta D/2{c_0}$ (Shaw Reference Shaw1975). In (3.8) ${S_0}$ denotes the cylinder surface and R the distance between a point on the cylinder surface ${\textbf {r}} = {(r,\theta ,z)^\textrm{T}}$ and the source point (i.e. integration variable) ${{\textbf {r}}_0}$, whereas ${n_0}$ denotes the inward surface normal direction at ${{\textbf {r}}_0}$. The last term in (3.8) is due to the pressure directly behind a possible shock front, and ${z_{ou}}$ thus denotes the axial distance to a source point at which the shock front arrives at the delayed time ${t_0} = t - R/c_{0}$. Finally, ${p_r}$ models the effect of the reflected incoming pressure, which is approximated as

(3.9)\begin{equation}{p_r}(\theta ,t) = \begin{cases} {\left( {{C_R}(\theta ,t) - \dfrac{{{C_R}(\theta ,t) - 2}}{{{t_{fade}}}}t} \right){p_{in}}(\theta ,t),}&{\cos \,\theta \gt 0 \cup t \lt {t_{fade}},}\\ {2,}&{\textrm{elsewhere}\textrm{.}} \end{cases}\end{equation}

It can be seen in (3.9) that a linear ‘fade-out’ function, using the parameter t _fade, was applied to the reflection coefficient ${C_R}$ to account for the set-up of an inertial flow over time. A linear form was chosen for the sake of simplicity.

Equation (3.8) can be solved with little computational effort by direct numerical approximation (Shaw Reference Shaw1975), yielding a value for the pressure at discrete points on the cylinder surface at every time step and thus direct information about the force due to diffraction and reflection.

3.2.3. History force model

In this section we will adapt the results published by Parmar et al. (Reference Parmar, Haselbacher and Balachandar2008) to yield an estimate for the force contribution on a circular cylinder due to flow acceleration. Parmar et al. (Reference Parmar, Haselbacher and Balachandar2008) computed the time-dependent history forces on cylinders and spheres at finite Mach numbers numerically, and provided a functional relationship using a convolution integral

(3.10)\begin{equation}{\boldsymbol{F}}_{hist} = - \int_{ - \infty }^t {K\frac{{\textrm{d}({m_{df}}\boldsymbol{v})}}{{\textrm{d}t}}\,\textrm{d}\left( {\frac{{{c_0}\chi }}{{D/2}}} \right)} .\end{equation}

Here $K = K({c_0}(t - \chi )/(D/2);\;M)$ denotes Mach number dependent kernel functions which were published in graphic form by Parmar et al. (Reference Parmar, Haselbacher and Balachandar2008). The quantities ${m_{df}}$ and ${\boldsymbol{v}}$ denote the time-dependent mass of the displaced fluid and the particle velocity of the incident flow, respectively, and $\chi$ is an integration variable. As the forces in Parmar et al. (Reference Parmar, Haselbacher and Balachandar2008) were computed for a previously fully developed flow and constant Mach numbers, several changes are necessary in order to predict forces during the highly transient scenario of a passing pressure wave.

Firstly, the force kernel K changes with Mach number, and thus, to account for this fact, the convolution integral is evaluated multiple times at different Mach numbers (e.g. for 20 time intervals) to account for the change of the incident flow conditions. Secondly, as the flow does not start from a fully developed state, the forces are multiplied by a ‘fade-in’ function $\beta $, defined, in linear form for the sake of simplicity, as

(3.11)\begin{equation}\beta (t) = \begin{cases} {\dfrac{t}{{{t_{fade}}}},}&{t \lt {t_{fade}},}\\ {1,}&{t \ge {t_{fade}}.} \end{cases}\end{equation}

Finally, the resulting forces are averaged over a small time span, $\Delta {t_{avg,hist}}$, to account for the inertia of the flow field for changing incident conditions. The history force is thus defined as

(3.12)\begin{equation}{F_{hist}}(t) = \frac{\beta }{{\mathrm{\Delta }{t_{avg,hist}}}}\int_{t - \mathrm{\Delta }{t_{avg,hist}}}^t {\int_{ - \infty }^{\tilde{t}} {K({{{\bar{M}}_{cyl}}} )\frac{{\textrm{d}({m_{df,cyl}}{{\bar{v}}_{cyl}})}}{{\textrm{d}\tilde t}}\,\textrm{d}\left( {\frac{{{c_0}\chi }}{{D/2}}} \right)\textrm{d}\tilde{t}} } ,\end{equation}

with

(3.13)\begin{equation}{m_{df,cyl}} = \frac{{{\rm \pi} {{\bar{\rho }}_{cyl}}{D^2}}}{4}.\end{equation}

Equation (3.12) can readily be integrated numerically with time steps small enough to sample the force kernel with sufficient accuracy. This was assured by limiting the time step to a one-hundredth of the wave duration and by sampling the non-zero portion of the force kernel with at least 150 points. The flow quantities $\bar{v}_{\textit{cyl}}, \bar{\rho}_{\textit{cyl}}, \bar{M}_{\textit{cyl}}$ are ‘effective’ flow quantities averaged over the whole cylinder surface, as explained in § 3.2.1.

3.2.4. Quasi-steady force model

The drag force on a circular cylinder due to a steady-state flow depends on both the governing Reynolds number and Mach number, such that

(3.14)\begin{equation}{F_{qs}}(t) = \frac{1}{2}{c_D}({Re,M} )\rho {v^2}D.\end{equation}

Ample literature exists on Reynolds and Mach number dependent drag coefficients $c_{D}$ (summarised in e.g. Hoerner Reference Hoerner1965; Blevins Reference Blevins1984). In this context, a simplified representation of the drag coefficient is used, as given in Blevins (Reference Blevins1984). Figure 8 depicts the assumed dependency of the drag coefficient on the governing Reynolds and Mach number. In this context, the drag forces are averaged over a short period of time, $\Delta {t_{avg,qs}}$, to account for the non-instantaneous change of the flow field with changing incoming particle velocity. The forces due to quasi-steady flow thus read

(3.15)\begin{equation}{F_{qs}}(t) = \frac{1}{{\Delta {t_{avg,qs}}}}\int_{t - \Delta {t_{avg,qs}}}^t {\frac{1}{2}{c_D}({{{\overline {Re} }_{cyl}},{{\overline{M}}_{cyl}}} ){{\bar{\rho }}_{cyl}}\bar{v}_{cyl}^2D\,\textrm{d}\tilde{t}} .\end{equation}

Figure 8.

Definition of steady-state drag coefficient for a circular cylinder for varying Mach and Reynolds numbers (Blevins Reference Blevins1984).

The averaging employed in (3.15) is executed numerically and the used flow quantities ${\bar{v}_{cyl}}$, ${\bar{\rho }_{cyl}}$, ${\overline {Re} _{cyl}}$, ${\overline{M}_{cyl}}$ are once again to be interpreted as ‘effective’ averaged flow quantities, as described in § 3.2.1.

3.2.5. Choice of parameters

The previously introduced parameters $\Delta {t_{avg,hist}}$, $\Delta {t_{avg,qs}}$, and ${t_{fade}}$ are meant to account for the inertia of the transient flow field with respect to the incident flow. Both the quasi-steady and history force models were initially developed for fully developed flow, and are rendered more flexible using the time averaging and fade-in functions.

Similarly, the reflection coefficient ${C_R}$ only holds for the reflection of a shock wave before an inertial flow has developed. Once the particles navigate along the object instead of being brought to an abrupt halt, this coefficient quickly diminishes before attaining the acoustic value of 2. In a preliminary parametric study, it was found that

(3.16)\begin{equation}\Delta {t_{avg,hist}} = \Delta {t_{avg,qs}} = {t_{fade}} = \frac{{{\rm \pi} D}}{{{c_0}}}\end{equation}

gives good agreement between the analytically and numerically obtained force histories. We note that ${\rm \pi} D/{c_0}$ equals the time a sound wave takes to orbit the cylinder once.

3.2.6. Validation of the analytical model

Results obtained using the new analytical modelling technique are now compared to results obtained with the numerical model. In figure 9 we present non-dimensional force histories predicted by numerical simulations and analytical calculations, as well as pressure contours at the moment of maximum drag load (as predicted by CFD).

Figure 9.

Comparison of the numerical and analytical force histories (a,c,e,g) and numerically obtained pressure contours at maximum load (b,d,f,h). Input parameters: ${\tau _i} = 30$, $R{e_i} = {10^4}$; (a,b) ${\alpha _r} = 0.0$, ${p_i}/{p_0} = 0.1$; (c,d) ${\alpha _r} = 0.5$, ${p_i}/{p_0} = 0.1$; (e,f) ${\alpha _r} = 0.0$, ${p_i}/{p_0} = 1.0$; (g,h) ${\alpha _r} = 0.5$, ${p_i}/{p_0} = 1.0$.

Figures 9(a) and 9(b) illustrate the loading of a circular cylinder by a shock wave of small overpressure. It can be seen that the maximum load occurs before the wave front has reached the midplane of the cylinder. The load amplitude is, therefore, mostly determined by the reflection and diffraction of the wave. We observe very good agreement between the analytical and the numerical model in terms of total drag load on the cylinder. Figure 9(a) also shows the contributions of the three force terms defined in (3.6), which confirms the dominance of the reflection–diffraction term. Figure 9(a) further shows the evolution of the impulse imparted on the cylinder

(3.17)\begin{equation}I = \int_0^t {F\,\textrm{d}\tilde{t}}\end{equation}

normalised by the incident impulse on the cross-section area of the cylinder

(3.18)\begin{equation}{i_i}D = D\int_0^{{\lambda _i}} {\rho (\zeta )v(\zeta )\,\textrm{d}\zeta } .\end{equation}

The initial distributions of density and particle velocity over the wavelength are given in appendix A. It can be seen that the maximum impulse imparted on the cylinder is only ~1/10 of the incident impulse. Further, we note that the analytical model predicts a slightly higher impulse compared to the numerical model.

Figures 9(c) and 9(d) show a case that differs from the previous case only in the rise coefficient ${\alpha _r}$, which is now 0.5, representing a pressure profile in the form of an isosceles triangle; this causes a reduction in maximum drag by approximately one order of magnitude. As can be seen in figure 9(d), the pressure gradients are much milder for this case, corresponding to a long non-dimensional rise time ${\alpha _r}{\tau _i}$. We note again good agreement between analytical and numerical model, although the peak force is overestimated by the analytical model in this case.

In figure 9(e–h) we present results for two cases that differ from the previous two cases in the pressure amplitude of the incoming wave. First, we note very good agreement in both cases between analytical and numerical model. Comparing the individual force contributions in figures 9(e) and 9(g), we find that while the shock wave case is again dominated by wave diffraction and reflection, in the finite rise time case all three force contributions have significant influence on the load amplitude. This can also be seen comparing figures 9(f) and 9(h), with the latter showing pressure contours similar to those encountered in steady-state flow, with vortex structures in the cylinder wake.

We proceed by assessing the loading intensity, in terms of force and impulse, obtained with the analytical and numerical models for the large data set defined in (3.5). It can be seen in figure 10(a) that the analytical model is in good agreement with the numerical model in terms of maximum drag load for the whole range of parameters explored, with predicted normalised peak loads spanning two orders of magnitude. Similarly, in figure 10(b) we observe broad agreement between the two approaches in terms of peak imparted impulse, with a slightly larger scatter. As can be seen in both figures 10(a) and 10(b), the analytical model tends to overpredict the load amplitudes, such that our analytical estimate can be considered slightly conservative.

Figure 10.

Correlation between numerical (subscript CFD) and analytical predictions (subscript an) of (a) the maximum force on the cylinder F _max and (b) the maximum transmitted impulse I _max.

We note that the proposed analytical model also needs a discretisation in space and time, as equations are not obtained in closed form. A comparison of the computational times using the set of simulations presented in figure 10 showed that the analytical calculation is much faster than the CFD simulations, with savings in computation time of several orders of magnitude. Solution times on a commercial workstation were of a few seconds in the case of the analytical model, while they were of several hours, or of a few days in the most demanding cases, for the CFD simulations.