1. Introduction
Understanding the nature and severity of transient loads exerted on objects by a surrounding fluid has been a concern for researchers over the past decades, due to the numerous safety-relevant applications of such knowledge in the defence, transport, energy and processing industries. Compared to loads resulting from steady-state flows of similar particle velocity, the loads exerted by a shock wave sweeping over a body can be up to one order of magnitude greater in amplitude (Tanno et al. Reference Tanno, Itoh, Saito, Abe and Takayama2003; Sun et al. Reference Sun, Saito, Takayama and Tanno2005). When a shock wave encounters a solid body, high pressure gradients are caused around the body due to the finite speed of wave propagation and due to the reflection of the wave from the forward-facing surfaces of the body. Transient loads due to shock waves passing over cylinders and spheres have been previously measured experimentally (Takayama & Itoh Reference Takayama and Itoh1985; Tanno et al. Reference Tanno, Itoh, Saito, Abe and Takayama2003; Sun et al. Reference Sun, Saito, Takayama and Tanno2005), computed numerically (Drikakis et al. Reference Drikakis, Ofengeim, Timofeev and Voionovich1997; Ofengeim & Drikakis Reference Ofengeim and Drikakis1997; Zółtak & Drikakis Reference Zółtak and Drikakis1998; Sun et al. Reference Sun, Saito, Takayama and Tanno2005; Luo et al. Reference Luo, Luo, Jin and Fan2017a,Reference Luo, Luo, Jin and Fanb) and modelled analytically (Friedman & Shaw Reference Friedman and Shaw1960; Shaw Reference Shaw1975; Parmar, Haselbacher & Balachandar Reference Parmar, Haselbacher and Balachandar2009).
Many existing studies were motivated by defence applications and focused on shock waves; at the other end of the spectrum, extensive literature exists on the propagation of sound waves and their interactions with solid bodies. The intermediate regime has received very little attention and will be the focus of this study, which is motivated by the growing need for predictive approaches to determine the forces exerted in deflagrations of mixtures of air and gas. In accidental deflagration events in congested environments, such as hydrocarbon processing plants, pressure waves of considerable duration and rise time (of the order of 0.1 s) can be emanated (American Petroleum Institute 2006; Det Norske Veritas 2010) and impinge on surrounding structures. Some of these structures can be assimilated to cylinders, for instance piping racks, pressure vessels or centrifugal compressor casings. At the moment, no accurate analytical methods exist to predict loads on structures by deflagration events (The Steel Construction Institute 2018). Here we aim at filling this gap.
The physical process of shock wave interaction with circular cylinders has been thoroughly investigated by other researchers. Due to the change in incidence angle, the initially regular reflection transforms into a Mach reflection during the progression of the shock front over the cylinder surface. The transition angle between regular and irregular reflection was shown to depend on the incident Mach number (Ben-Dor, Takayama & Kawauchi Reference Ben-Dor, Takayama and Kawauchi1980) as well as on the Reynolds number (Takayama & Sasaki Reference Takayama and Sasaki1983). In a combined experimental and numerical investigation Tanno et al. (Reference Tanno, Itoh, Saito, Abe and Takayama2003) and Sun et al. (Reference Sun, Saito, Takayama and Tanno2005) determined the transient forces on spheres of different diameters by sustained shock waves, i.e. waves of rectangular evolution in pressure. After an initial peak of high amplitude, which the authors attributed to the initial reflection of the shock wave, the drag was shown to rapidly decrease with time due to the equilibration of the pressure around the cylinder. After approximately 10–15 non-dimensional units of time (defined as sphere radius divided by the ambient speed of sound), the drag loads were found to agree well with reference values for the same particle velocity and Reynolds number in steady-state flow (Sun et al. Reference Sun, Saito, Takayama and Tanno2005).
For the case of arbitrary incident wave shapes and wavelengths, the evolution of overpressure and incident particle velocity can be of a highly unsteady nature. In addition to the reflection and diffraction of the incident wave and the establishment of an inertial flow field over time, the fluid experiences acceleration relative to the solid body. This acceleration is known to cause a significant contribution to the load in some cases, and is termed added-mass force in the case of incompressible flow. In this case the force on an object is linearly related to the relative acceleration between flow and solid, which was demonstrated to hold for a wide range of Reynolds numbers (Chang & Maxey Reference Chang and Maxey1995; Magnaudet & Eames Reference Magnaudet and Eames2000; Wakaba & Balachandar Reference Wakaba and Balachandar2007). In compressible flow this simple relation was shown to be inapplicable due to the finite speed of wave propagation (Miles Reference Miles1951; Longhorn Reference Longhorn1952), and the resulting force amplitudes were found to be significantly higher for finite Mach numbers (Parmar, Haselbacher & Balachandar Reference Parmar, Haselbacher and Balachandar2008).
It can thus be inferred from previous research that multiple physical phenomena cause force contributions to the loading of cylinders by arbitrary pressure waves. Magnaudet & Eames (Reference Magnaudet and Eames2000) suggested that the force on a particle immersed in unsteady flow can be categorised into five contributions, namely: quasi-steady, inviscid unsteady, viscous unsteady, lift and buoyancy-gravity, i.e.
\begin{equation}{\boldsymbol{F}}(t) = {\boldsymbol{F}}_{qs}(t) + {\boldsymbol{F}}_{iu}(t) + {\boldsymbol{F}}_{vu}(t) + {\boldsymbol{F}}_l(t) + {\boldsymbol{F}}_{bg}(t).\end{equation}It seems obvious to neglect the viscous unsteady and buoyancy-gravity driven forces in unsteady, high Reynolds number, high speed flow. However, simplified modelling techniques will need to consider at least quasi-steady and inviscid unsteady contributions. Parmar et al. (Reference Parmar, Haselbacher and Balachandar2009) presented a simple model for spheres, including pressure gradient, acceleration reaction and quasi-steady contributions, but neglected the reflection and diffraction of the wave and the effects of changing Mach numbers. Other simple methods, for example those widely used in industrial design guidelines (The Steel Construction Institute 2018), make first-order estimates of the individual force contributions and identify the dominant one as a function of the object's size. These methods are, however, extrapolated from simplified methods used in defence applications and yield inaccurate results in many cases, as we have recently shown elsewhere (Gauch et al. Reference Gauch, Bisio, Rossin, Montomoli and Tagarielli2019b).
In a parallel paper (Gauch et al. Reference Gauch, Bisio, Rossin, Montomoli and Tagarielli2019a), we focused on the transient loading of two-dimensional box-like objects loaded by the passage of pressure waves of arbitrary shape, amplitude and time duration. Here we aim at extending this investigation to the case of two-dimensional circular cylinders. We will develop analytical predictions, validate them by numerical calculations and present the results in the form of non-dimensional design maps of immediate use to design engineers.
2. Problem definition
A planar pressure wave of length 
${\lambda _i}$ in space, rise coefficient 
${\alpha _r}$ and maximum overpressure 
${p_i}$ is incident upon a rigid, fixed circular cylinder of diameter D, as in figure 1. The initial overpressure distribution along the wave is assumed to be defined piecewise linear, to give a triangular wave profile. The triangular shape is chosen for its capability of approximating both shock waves and pressure waves originated by deflagration events, but the models developed in this study are applicable to pressure wave of arbitrary shape. The surrounding medium, air, is characterised by the heat capacity ratio 
$\gamma = 1.4$, the specific gas constant 
$R = 287\,\textrm{Jk}{\textrm{g}^{ - 1}}\,{\textrm{K}^{ - 1}}$ and the ambient pressure and temperature 
${p_0}$ and 
${T_0}$, respectively. The objective of this study is to determine the transient load on the cylinder. As the tail of the incoming wave travels at the ambient speed of sound, 
${c_0}$, a loading duration can be quantified as
\begin{equation}{t_i} = \frac{{{\lambda _i}}}{{{c_0}}},\quad {c_0} = \sqrt {\gamma R{T_0}} .\end{equation}
Figure 1. Definition of the pressure profile impinging on a circular cylinder.
Dimensional analysis dictates that the problem at hand depends on the following non-dimensional groups:
\begin{equation}\left. {\begin{array}{c@{}} {\dfrac{{{p_i}}}{{{p_0}}},\quad {\alpha_r},\quad
{\tau_i} = \dfrac{{{t_i}}}{{D/{c_0}}} =
\dfrac{{{\lambda_i}}}{D},}\\ {R{e_i} =
\dfrac{{{\rho_i}{v_i}D}}{{{\mu_i}}},\quad \gamma ,\quad \tau
= \dfrac{t}{{D/{c_0}}},} \end{array}} \right\}\end{equation}
where 
${\rho _i}$, 
${v_i}$ and 
${\mu _i}$ denote the maximum density, particle velocity and dynamic viscosity of the fluid within the incident wave, respectively, and t denotes time. Functional relationships between 
${\rho _i}$, 
${v_i}$ and the pressure wave coefficients 
${p_i},\;{\alpha _r}$ are given in appendix A.
With reference to the non-dimensional groups defined in (2.2), we can visualise the wide range of problem parameters we aim to analyse in this study and mark domains covered by existing studies. Figure 2(a) depicts the range of possible shock wave cases and the domains covered by existing studies as well as by the models presented in this paper. It can be seen that most studies on shock waves correspond to very large wavelengths 
$({\tau _i} \to \infty )$, as the waves in those studies were defined as rectangular. Analytical models for the interaction of a pressure wave with a cylinder were proposed by Friedman & Shaw (Reference Friedman and Shaw1960) and Shaw (Reference Shaw1975). These models are limited to the acoustic region 
$({p_i}/{p_0} \to 0)$. In comparison, the models proposed here cover the whole range of wavelengths and we will prove their validity from the acoustic region to significant overpressures 
$(0 \lt {p_i}/{p_0} \lt 3)$.
Figure 2. Domain of validity of proposed model in comparison to existing studies for (a) shock waves 
$({\alpha _r} = 0)$ and (b) pressure waves of finite rise time 
$({\alpha _r} \gt 0)$.
Similarly, in figure 2(b) we compare ranges of applicability of results available in literature to the domain of validity of models developed in this study, for pressure waves of non-zero risetime. For very large non-dimensional rise times 
$({\alpha _r}{\tau _i} \to \infty )$ abundant literature is available, as this corresponds to the case of steady-state flow around a cylinder. The previously introduced models by Friedman & Shaw (Reference Friedman and Shaw1960) and Shaw (Reference Shaw1975) are, again, only applicable for very low overpressure ratios. Additionally, the range of validity of these models narrows for longer non-dimensional rise times, as the contribution of the wave diffraction to the maximum drag on the cylinder decreases, as will be discussed in detail below.
In figures 2(a) and 2(b) ranges of Reynolds number have not been included for simplicity. As for overpressure ratio, non-dimensional wavelength and risetime, to our knowledge there is no model available in literature that covers wide ranges of Reynolds numbers. Developing such models is, therefore, the main objective of this study.
3. Numerical and analytical models
In this section we present both the numerical and the analytical modelling approaches developed in this study. Firstly, the numerical methodology and a mesh convergence study are presented, then the new analytical model is described in detail and validated against the numerical predictions.
3.1. Numerical model
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
\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:
\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)
\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
\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 3. 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:
\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. Analytical model
As a complement to the detailed numerical model, we now present a semianalytical model which is able to capture the most important physics without the significant computational effort of the CFD model and fosters understanding of the nature of the loading in a wide range of scenarios.
The outset of the new model is the force superposition after Magnaudet & Eames (Reference Magnaudet and Eames2000), given in (1.1). Along the lines of a drag model for spheres, proposed by Parmar et al. (Reference Parmar, Haselbacher and Balachandar2009), the viscous unsteady, lift and buoyancy-gravity force contributions are neglected. The inviscid unsteady force 
${\boldsymbol{F}}_{iu}$ consists of two parts, namely a pressure gradient and a history term. Parmar et al. (Reference Parmar, Haselbacher and Balachandar2009) define the former to be due to pressure gradients in the flow which exist neglecting the presence of the object, and the latter to be due to the acceleration of the ambient fluid. We argue, however, that, in the present case, the pressure gradient force is more accurately described as force resulting from the reflection and diffraction of the incoming pressure wave. The force parametrisation therefore reads
\begin{align} \begin{array}{ll} \boldsymbol{F}(t) & = {\boldsymbol{F}}_{qs}(t) + {\boldsymbol{F}}_{iu}(t)\\ & = {\boldsymbol{F}}_{qs}(t) + {\boldsymbol{F}}_{diff}(t) + {\boldsymbol{F}}_{hist}(t), \end{array}\end{align} with 
$\boldsymbol{F}_{diff}$ denoting the force resulting from reflection and diffraction and 
$\boldsymbol{F}_{hist}$ denoting the history force term. In figure 6 we provide an overview of the different parts of the analytical model and the flow of the overall calculation. From a set of chosen input parameters, the propagation of the wave under the influence of compressibility is predicted, yielding temporal and spatial evolution of the relevant flow quantities. These are used as input to the three separate models for the three force contributions defined in (3.6). Finally, the overall force on the cylinder is obtained by superposing these three force contributions.
Figure 6. Flowchart of calculations for the proposed analytical model.
We will proceed by introducing individual modelling approaches for each of the parts depicted in figure 6.
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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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
\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.
4. Results and discussion
In the following, we construct non-dimensional maps from the results of numerical simulations. Subsequently, the analytical predictions are used to analyse the influence of the three force contributions (3.6) on the overall load, to shed light on the nature of the loading in different regimes of response (compare figure 2).
4.1. Numerical predictions
We start by analysing shock wave cases at various wavelengths, 
${\tau _i}$, pressure ratios, 
${p_i}/{p_0}$, and maximum Reynolds numbers, 
$R{e_i}$. In figure 11(a) we present the maximum non-dimensional force on a cylinder at 
$R{e_i} = {10^6}$ for a range of incident wavelengths. It can be seen that after a rapid increase in the range 
${\tau _i} \lt 10$, the maximum force values approach a pressure ratio dependent asymptote at approximately 
${\tau _i} = 40$. These values can, therefore, be seen as representative for triangular shock waves of very long wavelength, corresponding to rectangular waves. It is evident that the average pressure on the cylinder front 
${F_{max}}/({p_i}D)$ increases with increasing overpressure ratio due to the effects of compressibility (compare with (3.7)).
Figure 11. Maximum load on the cylinder for the shock wave cases 
$({\alpha _r} = 0)$. (a) Dependency of maximum load on the load duration 
${\tau _i}$ for 
$R{e_i} = {10^6}$. (b) Dependency of maximum load on Reynolds number for 
${\tau _i} = 60$.
Figure 11(b) illustrates the dependency of the maximum drag force on the maximum incident Reynolds number. It can be seen that for small Reynolds numbers the maximum force is greater than for higher Reynolds numbers, which has previously been shown experimentally (Takayama & Sasaki Reference Takayama and Sasaki1983) and can be attributed to an influence of viscous effects on the shock reflection pattern (Kleine et al. Reference Kleine, Timofeev, Hakkaki-Fard and Skews2014). On the other hand, in the regime 
$R{e_i} \gt {10^4}$, the Reynolds number can be seen to have only a small influence on the load, in accordance with findings by Kleine et al. (Reference Kleine, Timofeev, Hakkaki-Fard and Skews2014).
Next, in figure 12 we analyse the influence of the non-dimensional rise time 
${\alpha _r}{\tau _i}$ on the maximum load for various overpressure ratios and Reynolds numbers. In the range 
${\alpha _r}{\tau _i} \lt 1$, the contours can be seen to approach horizontal asymptotes corresponding to the maximum loads for a shock wave (dependent on pressure ratio and Reynolds number), which were given in figure 11. For intermediate rise times, 
$1 \lt {\alpha _r}{\tau _i} \lt 10$, the maximum load is found to be strongly dependent on the rise time. As the pressure gradient in the incoming wave becomes less severe, the maximum drag load diminishes. It can be seen that for most of the values of overpressure ratio and non-dimensional rise time shown, the lowest Reynolds number yields the highest drag loads on the cylinder, whereas the contours for 
$R{e_i} = {10^4}$ and 
$R{e_i} = {10^6}$ correspond closely. The contours can be seen to flatten out again for 
${\alpha _r}{\tau _i} \gt 10$. However, for the two higher Reynolds numbers, the general trend of decreasing maximum load with increasing non-dimensional rise time can be observed to reverse. This can be explained by the development of a quasi-steady flow field around the cylinder for long non-dimensional wavelengths. Periodic vortex shedding is triggered after a while, which leads to increased maximum drag loads, especially for higher Mach numbers (Rodriguez Reference Rodriguez1984; Xu, Chen & Lu Reference Xu, Chen and Lu2009; Xia et al. Reference Xia, Zuoli, Yipeng and Shiyi2016).
Figure 12. Contours of the maximum force on the cylinder as a function of the rise time, pressure ratio and Reynolds number.
In figures 13–15 the same data as in figure 12 is presented in alternative form for closer examination. Figure 13 depicts the normalised maximum force on a cylinder for a maximum Reynolds number of 
${10^2}$. It can be seen that up until 
${\alpha _r}{\tau _i} = 100$ the force amplitudes are decreasing monotonically for all investigated pressure ratios. Interestingly, the maximum force changes by a factor of more than 10 in the investigated range of rise times for the lower pressure ratios, whereas for the highest pressure ratio this factor is up to two. As a reference, the mean steady-state drag value is given for every overpressure ratio as 
$0.5{c_D}{\rho _i}v_i^2D$ with drag coefficients 
${c_D}$ as defined in figure 8. The maximum forces seem to approach these values asymptotically. It is to be expected that for higher wavelengths, the above described vortex shedding will develop and that higher maximum force values would, therefore, be recorded.
Figure 13. Variation of the maximum force on the cylinder with rise time, for 
$R{e_{i}} = {10^2}$. Values at infinity are steady-state drag values.
Figure 14. Variation of the maximum force on the cylinder with rise time for 
$R{e_i} = {10^4}$. Values at infinity are steady-state drag values.
Figure 15. Variation of the maximum force on the cylinder with rise time for 
$R{e_i} = {10^6}$. Values at infinity are steady-state drag values.
Figure 14 shows similar information but for a maximum Reynolds number of 
$R{e_i} = {10^4}$. The higher Reynolds number gives rise to the earlier mentioned reversal of the downward trend of the maximum forces. This effect can only be observed here for 
${p_i}/{p_0} \ge 1$, which can be explained by the higher shedding frequencies due to higher incident particle velocities and the shock waves triggered by an oscillating wake for higher Mach numbers (Xu et al. Reference Xu, Chen and Lu2009).
Finally, in figure 15 we present the results for the highest investigated Reynolds number, 
$R{e_i} = {10^6}$. Comparing figures 14 and 15, it can be seen that the increase of the maximum force is triggered even earlier than for the case of 
$R{e_i} = {10^4}$, whereas similar values for the maximum force are obtained for low non-dimensional rise times.
We proceed by comparing the maximum impulse imparted on a cylinder to the incoming impulse for a wide range of wavelengths. Figure 16(a) illustrates the peak imparted impulse on a cylinder by shock wave loading. It can be seen that relative to the incident impulse the highest impulses are recorded for the shortest wavelengths. This can be explained by higher gradients in the incoming wave and a more prominent diffraction of waves of short wavelengths. For increasing wavelengths, the impulse can be seen to decrease in amplitude and only small changes are observed in the range 
${\tau _i} \gt 20$.
Figure 16. Variation of the maximum impulse exerted on the cylinder with wavelength.
In figure 16(b) the same information is shown for a rise coefficient of 
${\alpha _r} = 0.5$ and a wider range of wavelengths. It can be seen that for 
${\tau _i} \to 0$ the normalised imparted impulse attains values slightly above unity, whereas for large wavelengths, values one order of magnitude smaller are predicted. In contrast to the shock wave case, given in figure 16(a), the normalised impulses decrease until 
${\tau _i} \approx 50$ or longer, before approaching a steady value. This steady value can be approximated by calculating the imparted impulse on a cylinder due to purely quasi-steady drag. In first approximation, neglecting wave distortion, the impulse on a cylinder due to the quasi-steady flow caused by a pressure wave is
\begin{align}{I_{qs}} & = \int_0^{{t_i}} {\dfrac{1}{2}{c_D}(t)\rho (t){v^2}(t)D\,\textrm{d}t} \nonumber\\ & = \dfrac{D}{{2{c_0}}}\int_0^{{\lambda _i}} {{c_D}(\zeta )\rho (\zeta ){v^2}(\zeta )\,\textrm{d}\zeta } . \end{align}
In (4.1), 
${c_D}(\zeta )$ denotes the drag coefficient as in figure 8, with 
$M(\zeta )$ and 
$Re(\zeta )$, which can be evaluated using the equations given in appendix A. In contrast to the maximum forces, given in figures 13–15, the quasi-steady impulse can be seen to be a mostly conservative estimate of the imparted impulse predicted numerically for long wavelengths. This can be explained by the periodicity of the wake oscillations described earlier, which tend to balance out when integrated over time.
4.2. Analytical predictions
We now use the results obtained with the analytical model to draw conclusions on the importance of the individual force contributions (3.6) in different regimes of input parameters.
In figure 17 we present a contour map of the relative importance of the three contributions for a Reynolds number of 
${10^6}$ and the whole examined ranges of non-dimensional rise time and pressure ratio. It is evident that for all pressure ratios the contribution from diffraction and reflection is dominant for 
${\alpha _r}{\tau _i} \to 0$ and accounts for more than 60 % of the peak load, up until 
${\alpha _r}{\tau _i} = 10$. In this regime the contribution of the history force can be seen to reach a peak of over 25 % for the higher pressure ratios. For increasing wavelengths the quasi-steady contribution is found to gain more and more importance and accounts for over 75 % of the peak load at 
${\alpha _r}{\tau _i} \gt 100$, 
${p_i}/{p_0} \gt 1$. For lower pressure ratios the quasi-steady force is found to be less dominant.
Figure 17. Contours of the force contributions 
$F_{<>}$ to the maximum force according to the analytical model, for 
$R{e_i} = {10^6}$.
In order to examine the influence of the Reynolds number, figure 18 depicts a slice through figure 14 with an added line for 
$R{e_i} = {10^2},\;{10^4}$. The same trends as in figure 17 can be observed, with the diffraction contribution dominating for 
${\alpha _r}{\tau _i} \to 0$ and the quasi-steady contribution gaining more and more significance for rising non-dimensional rise times. Due to the higher drag coefficients for lower Reynolds numbers (compare to figure 8), the quasi-steady contribution is stronger for the lower Reynolds numbers.
Figure 18. Variation of the force contributions with rise time for 
${p_i}/{p_0} = 1.0$.
4.3. Discussion
The results presented in this section can directly be used in industrial design for arbitrary pressure wave loading, for example for the design of cylindrical structures common in hydrocarbon processing plants, opening new opportunities for structural optimisation. Figures 11–16 give access to peak loads for a wide range of input parameters without any further calculation. This constitutes a substantial improvement over the current industrial design practice, where methods and charts developed for shock waves are used (The Steel Construction Institute 2018), giving high levels of uncertainty and, in general, unnecessarily conservative predictions.
An increase in maximum drag load has been found for the higher examined Reynolds numbers at large non-dimensional rise times, which can be explained by the unsteady vortex shedding from the cylinder. It depends on the structural design case under assessment if these higher load frequencies need to be considered. Furthermore, it can be expected that the assumption of two-dimensionality of the flow and the use of the URANS equations in the numerical solver have a detrimental effect on the solution accuracy in these regimes. At high Reynolds numbers the cylinder wake becomes turbulent and transitions from laminar to turbulent flow appear around the cylinder. These can only be captured coarsely by the approach employed here.
A further investigation of this effect lies yet outside the scope of this work and studies have been published by other authors (e.g. Rodriguez Reference Rodriguez1984; Xu et al. Reference Xu, Chen and Lu2009; Xia et al. Reference Xia, Zuoli, Yipeng and Shiyi2016).
The newly proposed analytical model provides transient load histories with little computational effort and has been shown to agree well with the more detailed numerical model. As the load histories caused by a triangular wave can take complex evolutions in time, as can be seen in figure 9, this model can contribute valuable additional information in a design process, and costly numerical computations can be avoided.
It has to be borne in mind that the assumption of two-dimensionality was made throughout this study. Cylinders of finite width experience yet, in general, lower drag loads than cylinders of infinite width (assumption of two-dimensionality). The results presented here can, therefore, be seen as upper bounds to the loads experienced by structures of finite width.
Finally, we note that models similar to the analytical model proposed here can be developed for other simple geometrical shapes, as we have already shown for a two-dimensional box-like structure in a companion study (Gauch et al. Reference Gauch, Bisio, Rossin, Montomoli and Tagarielli2019a).
5. Conclusions
We have developed a numerical and an analytical modelling approach to predict the transient loading on circular cylinders by pressure waves of arbitrary shape, amplitude and time duration. The analytical modelling approach was validated using the more detailed numerical approach, and its asymptotic behaviour was explored. The main conclusions of this study are as follows:
(i) A predictive tool was developed to allow estimates of the transient loading histories induced on circular cylinders by arbitrary pressure waves.
(ii) Using the numerical model, a large parametric study was conducted and the results were synthesised in the form of non-dimensional design charts of immediate application.
(iii) The maximum drag loads exerted on cylinders by pressure waves of finite rise time were shown to depend primarily on the non-dimensional rise time
${\alpha _r}{\tau _i}$ and the overpressure ratio ${p_i}/{p_0}$.(iv) It was shown that for short non-dimensional rise times, the loads are most severe and mainly governed by reflection and diffraction, whereas for large non-dimensional rise times the loads approach quasi-steady nature. In between these two extremes, all three force contributions (3.6) were shown to be significant.
Acknowledgements
We are grateful to the helpful suggestions provided by Dr G. Rigas and to M. Ruggiero of B. Hughes for managing the financial side of the project.
Declaration of interests
The research was funded by the Turbomachinery Engineering team of Baker Hughes. https://www.bakerhughes.com/.
Appendix A
At 
$t = 0$, the given triangular pressure wave implies a distribution in terms of overpressure 
$p = {p_{abs}} - {p_0}$ is
\begin{equation}p(\zeta ) = {\begin{cases} {\dfrac{{{p_i}}}{{{\alpha_r}{t_i}{c_0}}}\zeta ,}&{0 \le \zeta \le {\alpha_r}{t_i}{c_0},}\\ {\dfrac{{{p_i}}}{{{t_i}(1 - {\alpha_r}){c_0}}}({t_i}{c_0} - \zeta ),}&{{\alpha_r}{t_i}{c_0} \lt \zeta \le {t_i}{c_0},}\\ {0,}&{\zeta \gt {t_i}{c_0},} \end{cases}}\end{equation}
where 
$\zeta $ is a spatial coordinate pointing from the wave front to the wave tail, with 
$\zeta = 0$ at the wave front. For the case of a negligible rise time coefficient 
${\alpha _r} = 0$, the pressure wave is a shock wave and the properties behind the shock front are determined by the Rankine–Huguinot equations for a perfect gas (e.g. Liepmann Reference Liepmann1957)
\begin{equation}\left. {\begin{array}{c@{}} {{\rho_{shock}} = {\rho_0}\dfrac{{2\gamma \textrm{ } + \textrm{ }(\gamma + 1)\dfrac{{{p_i}}}{{{p_0}}}}}{{2\gamma \textrm{ } + \textrm{ }(\gamma - 1)\dfrac{{{p_i}}}{{{p_0}}}}},\quad {T_{shock}} = \dfrac{{{p_i} + {p_0}}}{{R{\rho_{shock}}}},}\\ {{v_{shock}} = \dfrac{{\dfrac{{{p_i}}}{{{p_0}}}\sqrt {\dfrac{{{p_0}}}{{{\rho_0}}}} }}{{\sqrt {\dfrac{{\gamma + 1}}{2}\dfrac{{{p_i}}}{{{p_0}}} + \gamma } }},\quad {M_{shock}} = \dfrac{{{v_{shock}}}}{{\sqrt {\gamma R{T_{shock}}} }}.} \end{array}} \right\}\end{equation}Further behind the shock front it can be assumed that the gas undergoes isentropic expansion
\begin{equation}\frac{p}{{{\rho ^\gamma }}} = \textrm{const}\textrm{.}\end{equation}Therefore, the density, temperature and velocity fields can be calculated as
\begin{equation}\left. {\begin{array}{c@{}} {\rho (\zeta ) = {\rho_0}{{\left( {\dfrac{{p(\zeta )}}{{{p_i} + {p_0}}}} \right)}^{1/\gamma }},\quad T(\zeta ) = \dfrac{{p(\zeta )}}{{\rho (\zeta )R}},}\\ {v(\zeta ) = {v_{shock}} - \dfrac{2}{{\gamma - 1}}\left( {\sqrt {\gamma R{T_{shock}}} - \sqrt {\gamma R{T_{shock}}(\zeta )} } \right),\quad M(\zeta ) = \dfrac{{v(\zeta )}}{{\sqrt {\gamma RT(\zeta )} }}.} \end{array}} \right\}\end{equation}
Similarly, for the case of a non-negligible rise time coefficient 
${\alpha _r}$, the density, temperature and velocity fields can be calculated as
\begin{equation}\left. {\begin{array}{c@{}} {\rho (\zeta ) = {\rho_0}{{\left( {\dfrac{{p(\zeta )}}{{{p_0}}}} \right)}^{1/\gamma }},\quad T(\zeta ) = \dfrac{{p(\zeta )}}{{\rho (\zeta )R}},}\\ {v(\zeta ) = \dfrac{2}{{\gamma - 1}}\left( {\sqrt {\gamma RT(\zeta )} - {c_\textrm{0}}} \right),\quad M(\zeta ) = \dfrac{{v(\zeta )}}{{\sqrt {\gamma RT(\zeta )} }}.} \end{array}} \right\}\end{equation}
