%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
The problem of flexural-gravity wave propagation in ice covers has been actively investigated for the last decades in the context of the development of Polar territories and the possibility of long-term use of the Northern Sea Route. One of the first considerations of flexural-gravity waves was done by Greenhill (Reference Greenhill1886). Models of waves in ice cover were further developed by Kheysin (Reference Kheysin1963) where ice was modelled by a linear elastic plate. Reviews of works on hydroelastic waves in unbounded regions can be found in Kheysin (Reference Kheysin1963), Squire et al. (Reference Squire, Hosking, Kerr and Langhorne1996) and Squire (Reference Squire2020). The model of a linear elastic plate is still the main one, although viscoelastic (Hosking et al., Reference Hosking, Sneyd and Waugh1988; Wang & Shen Reference Wang1988) and poroelastic models (Meylan et al., Reference Meylan, Bennetts and Peter2017; Chanda et al., Reference Chanda, Barman, Sahoo and Meylan2024) plates are also considered. Nonlinear models (Plotnikov and Toland, Reference Plotnikov and Toland2011; Dinvay et al., Reference Dinvay, Kalisch and Parau2019), models that take into account ice compression (Stepanyants and Sturova, Reference Stepanyants and Sturova2020; Sturova, Reference Sturova2024; Stepanyants and Sturova, Reference Stepanyants and Sturova2021a ), the presence of shear flows (Tkacheva, Reference Tkacheva2024) or the presence of cracks in the plates (Tkacheva, Reference Tkacheva2019) also have been actively developed. The problems of propagation of hydroelastic waves in a channel were studied by the authors of the article in different formulations and for different ice conditions in Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014), Khabakhpasheva et al. (Reference Khabakhpasheva, Shishmarev and Korobkin2019) and Shishmarev et al., (Reference Shishmarev, Khabakhpasheva and Korobkin2016, Reference Shishmarev, Khabakhpasheva and Oglezneva2023, Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2024, Reference Shishmarev, Sibiryakova, Naydenova and Khabakhpasheva2022). The problems on methods of ice destruction, including methods involving hydroelastic waves generated by moving loads, were studied by Ni et al. (Reference Ni2024), Yuan et al., (Reference Yuan, Ni, Wu, Xue and Han2022, Reference Yuan, Ni, Wu, Lu and Xue2024) and Zemlyak et al., (Reference Zemlyak, Kozin and Vasilyev2022, Reference Zemlyak, Pogorelova and Kozin2023).
The main part of the flexural-gravity wave research was carried out for ice sheets of infinite extent. Studies on ice cover in channels belong to the broader category of ice–fluid–coastal structure interaction problems. The majority of the ice tanks, where scientific and technical experiments with ice cover are conducted, have finite dimensions and rectangular cross-sections. For this reason, studying the features of waves in channels and how they differ from waves in unbounded sheets is highly important. The motivation for investigating such problems also arises from the analysis of these experiments and from studying safe movement across frozen water reservoirs or intentional ice breaking to prevent flooding; see, for example, Kozin et al. (Reference Kozin, Chizhumov and Zemlyak2010, Reference Kozin, Rogozhnikova and Zemlyak2017), Pogorelova et. al. (Reference Pogorelova, Zemlyak and Kozin2019a , Reference Pogorelova, Zemlyak and Kozinb ) and Zemlyak et al. (Reference Zemlyak, Baurin, Kozin, Chizhiumov and Matiushina2019). The boundary conditions and parameters of heterogeneous ice significantly affect the results of these studies. Hydroelastic waves in a channel are a combination of waves propagating along the channel and sloshing waves across the channel due to multiple reflections from the channel walls. The distribution of ice deflections and deformations in a channel can differ significantly from those in unbounded ice sheets, see Shishmarev et al. (Reference Shishmarev, Khabakhpasheva and Korobkin2016). Channel problems can be divided into two categories: studies on periodic free waves and those on forced waves induced by external loads. The first category includes papers by Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014), Ren et al. (Reference Ren, Wu and Li2020), Shishmarev et al. (Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2022) and Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2015, Reference Batyaev and Khabakhpasheva2022, Reference Batyaev and Khabakhpasheva2024). These papers study wave characteristics (dispersion relations, phase and group velocities), as well as wave profiles and strain distributions in the ice cover. Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014) and Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2015) investigated ice covers of constant thickness, assuming clamped or free-edge boundary conditions at the channel walls, respectively. Shishmarev et al. (Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2022) and Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2022) studied the effect of the ice cover with a linearly varying thickness, with cases of symmetric linear variation relative to the channel or continuous linear change from one side to the other. Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2024) studied a similar problem for an ice cover with non-uniform compression. The second category involves hydroelastic waves in a channel induced by moving external loads, either underwater or along the upper surface of ice, see Shishmarev et al., (Reference Shishmarev, Khabakhpasheva and Korobkin2016, Reference Shishmarev, Khabakhpasheva and Oglezneva2023, Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2024), Khabakhpasheva et al. (Reference Khabakhpasheva, Shishmarev and Korobkin2019) and Zeng et al., (Reference Zeng, Korobkin, Ni and Xue2021, Reference Zeng, Korobkin, Ni and Xue2022). In Shishmarev et al. (Reference Shishmarev, Khabakhpasheva and Korobkin2016), a general solution was obtained for a load moving at constant speed along the central line of the channel. Shishmarev et al., (Reference Shishmarev, Khabakhpasheva and Oglezneva2023, Reference Shishmarev, Sibiryakova, Naydenova and Khabakhpasheva2024) analysed the effects of linear ice thickness variation and porosity during the movement of external loads or the oscillations of a submerged sphere with constant frequency. Zeng et al., (Reference Zeng, Korobkin, Ni and Xue2021, Reference Zeng, Korobkin, Ni and Xue2022) investigated the effects of leads and cracks on ice cover oscillations in a channel. The impact of cracks of arbitrary shapes and locations in an infinitely extended ice sheet was studied by Li et al. (Reference Li, Wu and Ren2020). All these findings demonstrated the significant influence of ice inhomogeneities.
Ice inhomogeneities, in general, lead to complex change in the ice behaviour. A significant part of the literature is devoted to studying non-uniform compression in different directions within an ice sheet, see Stepanyants & Sturova (Reference Stepanyants and Sturova2021a ,Reference Stepanyants and Sturova b ), Eyov et al. (Reference Eyov, Klar, Kadri and Stiassnie2013), Ilichev (Reference Ilichev2016) and Sturova (Reference Sturova2021a ). Tensile and compressive forces can develop within an ice cover due to influences such as wind, water currents and temperature fluctuations (Evans and Peter, Reference Evans and Peter2011). In three-dimensional scenarios, these forces are as follows: longitudinal; transverse; combined (Timoshenko and Woinowsky-Krieger, Reference Timoshenko and Woinowsky-Krieger1959). The behaviour of free hydroelastic waves under non-uniform compression conditions was analysed by Bukatov and Goncharov (Reference Bukatov and Goncharov1989), the steady wave motion under external loading was studied by Bukatov et al. (Reference Bukatov, Zharkov and Zav’yalov1991). Sturova (Reference Sturova2021b ) calculated the limiting values of compression parameters required for a correct mathematical model. For the frozen channel(see Batyaev and Khabakhpasheva Reference Batyaev and Khabakhpasheva2024), it was shown that the force of compression and compression direction both have a major effect on the propagation of the first two waves in terms of the lowest frequency. They aslo have a strong effect on both the shape and speed of propagation and on the strain distribution within the ice plate.
The importance of varying ice thickness has been demonstrated in many works, but most of these studies were focused on geophysical investigations (observations and descriptions) in the shelf zone, see a good review in Liang et al. (Reference Liang, Pitt and Bennetts2024) for example. Many studies are devoted to the propagation of waves along the direction of thickness variation in a one-dimensional plate (beam) in a two-dimensional problem. In our study, the thickness varies perpendicular to the wave direction, and we consider waves in a two-dimensional plate (three-dimensional formulation of the problem). The hydroelastic waves, propagating along the ice-covered channel with linearly varying thickness in three-dimensional formulation were studied by Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2022) and Shishmarev et al. (Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2022). It was found that, as in the case of a constant-thickness ice cover in a channel (see Korobkin et al. Reference Korobkin, Khabakhpasheva and Papin2014; Batyaev and Khabakhpasheva Reference Batyaev and Khabakhpasheva2015), an infinite number of hydroelastic waves can propagate along the channel with identical wavelengths but differing transverse profiles. These waveforms were obtained, and it was shown that the points of maximum ice deflection shift towards the thinner regions of the plate. For long waves the dispersion curves and, therefore, the phase and group velocities, remain almost indistinguishable from those of a plate with constant thickness. As the wavenumber increases, the corresponding frequencies decrease relative to the constant-thickness scenario. The results for considered channel and sets of ice properties showed that for long waves the maximum strain occurs near the channel walls, while for short waves it shifts towards the interior of the channel.
The problem of flexural-gravity wave scattering by an ice sheet of variable geometry in the presence of uniform currents was studied by Aggarwal et al. (Reference Aggarwal, Barman, Martha and Tsai2024). The numerical results described a resonance phenomenon, the effect of plate shape and the intricate relationship between wave characteristics and varying plate properties. The hydrodynamic behaviour of a floating flexible plate of variable flexural rigidity connected to the seabed by a spring–damper system was studied by Michele et al., (Reference Michele, Zheng, Renzi, Borthwick and Greaves2024a ,Reference Michele, Zheng, Renzi, Borthwick and Greaves b ). Influence of ice-sheet geometry and supraglacial lakes on seasonal ice-flow was studied on large scales by Joughin et al. (Reference Joughin, Das, Flowers, Behn, Alley, King, Smith, Bamber, van den Broeke and Van Angelen2013). The propagation of long waves through water of variable shallow depth with a surface covered by fragmented broken ice was studied by Dafydd and Porter (Reference Dafydd and Porter2024). It was proposed that random variations in ice thickness could be a mechanism for the attenuation of waves propagating through broken sea ice. The wave attenuation rate of ocean surface waves in the marginal ice zone was studied by Meylan et al. (Reference Meylan, Bennetts, Mosig, Rogers, Doble and Peter2018). It was also discussed that the difference in the magnitude of the attenuation rates could be caused by differences in ice properties, especially ice thickness. Wave scattering by an ice floe of variable thickness, modelled as an elastic plate, was studied in Smith and Meylan (Reference Smith and Meylan2011). The model of ice shelf flexural vibrations in response to ocean waves, in which both the ice shelf thickness and the seabed beneath the shelf vary over distance, was studied in Meylan et al. (Reference Meylan, Ilyas, Lamichhane and Bennetts2021).
The main focus of this paper is on the effect of variable thickness across the channel on the characteristics of periodic waves. Firstly, the ice thickness shape is approximated by linear functions. Then, a method of decomposing the vibrations of a beam with an approximated thickness into normal modes (eigenfunctions) is employed. The thickness of the ice plate is approximated ensuring that the thickness remains a continuous function. Normal modes of ice vibrations are derived using trigonometric, hyperbolic and Bessel functions for the approximated thickness and provide exact solution.
The method presented is semianalytical. For a given ice thickness across the channel, a piecewise linear approximation keeping its continuity is used. Normal modes of an elastic beam with the approximated thickness are introduced. The accuracy of this method for the piecewise linear geometry depends solely on the precision of the eigenvalue calculations. After these eigenvalues were calculated the remaining part of the problem was solved analytically. This article demonstrates that, due to the effect of variable ice thickness, the transverse wave profiles do not always coincide with the shape of the corresponding mode (as was the case for ice of constant thickness, see Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014)), but rather represent a complex combination of them, depending on the wavenumber.
Several approaches exist for describing the vibrations of a thin finite plate with variable parameters. In Michele et al., (Reference Michele, Zheng, Renzi, Borthwick and Greaves2024a
,Reference Michele, Zheng, Renzi, Borthwick and Greaves
b
), a plate with piecewise constant rigidity
$D_m$
was studied. In that method, continuity of deflections, slopes, bending moments and shear forces was maintained at the boundaries of the intervals with constant rigidity. However, this approach does not ensure the continuity of stresses or strains in a plate for any finite set of intervals with piecewise constant rigidity. The stress field in the ice plate generated by hydroelastic waves is a key parameter investigated in a substantial number of studies. Piecewise constant rigidity results in stress jumps, as demonstrated in Korobkin et al. (Reference Korobkin, Khabakhpasheva and Malenica2017), and continuity of both stresses and bending moments cannot be satisfied together. The stress should become continuous in the limit that the length of the intervals used to approximate the variable parameters tends to zero. However, this significantly increases the number of intervals and one needs to be careful with possible issues of calculations. In Smith and Meylan (Reference Smith and Meylan2011) and Meylan et al. (Reference Meylan, Ilyas, Lamichhane and Bennetts2021), the eigenfunctions for the vibrations of ice with variable thickness were derived using expansion of them on ‘classic’ eigenfunctions for constant thickness beams and also using variational methods for determination of expansion coefficients. Both deflections and strains were investigated.
This article presents a method for piecewise linear thickness that ensures not only the continuity of deflections, slopes, bending moments and shear forces, but also the continuity of stresses and strains. In these terms, the method is closer to that used by Smith and Meylan (Reference Smith and Meylan2011) and Meylan et al. (Reference Meylan, Ilyas, Lamichhane and Bennetts2021). In Meylan’s work, in contrast, the thickness can be arbitrary, but the modes themselves are approximated. In § 4.3, we compare our results with the results obtained by the direct use of normal modes of a plate of constant thickness and demonstrate that such approach cannot always provide accurate results of the vibrations of an ice cover with piecewise linear thickness in the considered problem.
Two different types of thickness variation are considered: thickening and thinning near the channel walls. A study has been carried out about how these thickness variations affect the waveforms, their frequency characteristics and the resulting stresses in the plate. Knowledge of the influence of the thickness non-uniformity on the wave characteristics makes it possible to solve the reverse problem as a practical application. It will be possible to predict the thickness of the ice cover from the non-uniformity of hydroelastic waves travelling through the ice field .
2. Formulation of the problem
Periodic hydroelastic waves propagating along a frozen channel with constant speed and constant amplitude are considered. The amplitude of the considered waves is much smaller than their wavelength and the characteristic size of the channel. The scheme of the problem is shown in figure 1. The channel is of rectangular section with a finite depth
$H$
(
$-H \lt \hat {z} \lt 0$
) and a finite width
$2b$
(
$-b \lt \hat {y} \lt b$
), the channel is of infinite extent in the
$\hat {x}$
direction. Here,
$O\hat {x}\hat {y}\hat {z}$
is a Cartesian coordinate system. The channel is filled with liquid of density
$\rho _\ell$
. The liquid is inviscid, incompressible and covered with a thin ice cover. The ice is modelled by a thin elastic plate with given constant density
$\rho _i$
and variable across the channel rigidity
$\hat {D}(\hat {y})$
, where
$\hat {D}(\hat {y}) = E \hat {h}_i^3(\hat {y}) / [12 (1- \mu ^2)]$
,
$E$
is the Young’s modulus for ice,
$\mu$
is the Poisson ratio for ice and
$ \hat {h}_i (\hat {y})$
is a thickness of the ice cover. In general, the thickness of the ice cover can vary arbitrarily, but one condition must hold: the thickness changes small relatively to the width of the channel. For any given thickness of the ice under this condition we shall propose a solution method. All calculations and the solution method will be presented for two specific cases of the ice thickness distribution, including six subcases, described further. The ice thickness is constant along the channel in all cases.
Scheme of the problem of hydroelastic wave propagation along an ice channel. Variable thickness of ice across the channel is considered.
The problem is studied within the linear theory of hydroelasticity. The ice cover is modelled by a thin elastic plate. Deflections of the ice plate are governed by the Kirchhoff–Love theory and are described by the vertical displacement of the midsurface of the plate,
$\hat {z} = \hat {w}(\hat {x}, \hat {y}, \hat {t})$
. The mathematical formulation of the problem is coupled and consists of elastic and hydrodynamic parts. The elastic part is related to the determination of the ice deflections. The second part of the problem is related to the description of the hydrodynamics in the channel and determination of the flow potential,
$\hat {\varphi }(\hat {x}, \hat {y}, \hat {z}, \hat {t})$
. The ice deflections
$\hat {w}(\hat {x}, \hat {y}, \hat {t})$
satisfy the equation of a thin elastic plate with variable thickness across the channel. We consider the ice frozen to the channel walls, which is modelled by clamped conditions at
$\hat {y} = \pm b$
. In this paper, we restrict ourselves to solving the problem with the clamped conditions at the channel walls. The potential
$\hat {\varphi }(\hat {x}, \hat {y}, \hat {z}, \hat {t})$
satisfies the Laplace equation in the channel domain and boundary conditions of impermeability at rigid boundaries of the channel (walls and bottom) and linearized kinematic condition at the ice–liquid interface. Governing equations of the elastic part of the problem read
\begin{equation} \hat {h}_i (\hat {y}) \rho _i \frac {\partial ^2 \hat {w}}{\partial \hat {t}^2} - \frac {\partial \hat {Q}_x}{\partial \hat {x}} - \frac {\partial \hat {Q}_y}{\partial \hat {y}} = \hat {p}(\hat {x},\hat {y},0,\hat {t}) \quad (-\infty \lt \hat {x}\lt \infty, -b\lt \hat {y}\lt b, z = 0), \end{equation}
where
$\hat {Q}_x(\hat {x},\hat {y})$
and
$\hat {Q}_y(\hat {x},\hat {y})$
are transverse shears
\begin{equation} \hat {Q}_x = \frac {\partial \hat {M}_x}{\partial \hat {x}} + \frac {\partial \hat {M}_{xy}}{\partial \hat {y}}, \qquad \hat {Q}_y = \frac {\partial \hat {M}_{yx}}{\partial \hat {x}} + \frac {\partial \hat {M}_y}{\partial \hat {y}}. \end{equation}
Corresponding bending
$\hat {M}_x(\hat {x},\hat {y})$
,
$\hat {M}_y(\hat {x},\hat {y})$
and twisting moments
$\hat {M}_{xy}(\hat {x},\hat {y})$
caused by elastic forces are
\begin{equation} M_x = -\hat {D} \left ( \frac {\partial ^2 \hat {w}}{\partial \hat {x}^2} + \mu \frac {\partial ^2 \hat {w}}{\partial \hat {y}^2} \right), \quad M_y = -\hat {D} \left ( \frac {\partial ^2 \hat {w}}{\partial \hat {y}^2} + \mu \frac {\partial ^2 \hat {w}}{\partial \hat {x}^2} \right), M_{xy} = -(1-\mu) \hat {D} \frac {\partial ^2 \hat {w}}{\partial \hat {x} \partial \hat {y_{_{_{_{\!}}}}}}. \end{equation}
Here,
$\hat {p}(\hat {x},\hat {y},0,\hat {t})$
is the liquid pressure at the ice–liquid interface, determined by the linearized Bernoulli integral
\begin{equation} \hat {p}(\hat {x},\hat {y},0,\hat {t}) = - \rho _\ell \frac {\partial \hat {\varphi }}{\partial \hat {t}}(\hat {x},\hat {y},0,\hat {t}) - \rho _\ell g \hat {w}, \end{equation}
where
$g$
is the acceleration of gravity. These equations are closed by the clamped conditions of the ice deflections at the walls
\begin{equation} \hat {w} = 0, \quad \hat {w}_{\hat {y}} = 0, \qquad (\hat {y} = \pm b). \end{equation}
There are no initial conditions for the ice deflections and no conditions at
$\hat {x} \to \infty$
as they are not required for the periodic waves.
Governing equations for the hydrodynamic part of the problem read
\begin{equation} \frac {\partial ^2 \hat {\varphi }}{\partial \hat {x}^2} + \frac {\partial ^2 \hat {\varphi }}{\partial \hat {y}^2} + \frac {\partial ^2 \hat {\varphi }}{\partial \hat {z}^2} = 0 \quad (-\infty \lt \hat {x}\lt \infty, -b\lt \hat {y}\lt b, -H\lt \hat {z}\lt 0), \end{equation}
including the boundary conditions of impermeability and the kinematic condition
\begin{equation} \frac {\partial \hat {\varphi }}{\partial \hat {y}} = 0 \;\; (\hat {y} = \pm b), \quad \frac {\partial \hat {\varphi }}{\partial \hat {z}} = 0 \;\; (\hat {z} = -H), \quad \frac {\partial \hat {\varphi }}{\partial \hat {z}} = \frac {\partial \hat {w}}{\partial \hat {t}} \;\; (\hat {z} = 0). \end{equation}
The kinematic boundary condition at the ice–liquid interface is given in the linearized form. The original condition states that the normal velocities of ice and liquid at the contact boundary must be equal
\begin{equation} \nabla \hat {\varphi } \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{v}_{plate} \boldsymbol{\cdot} \boldsymbol{n} \qquad (\hat {z} = \hat {w} - \hat {h}_i(\hat {y}) /2). \end{equation}
Note that we consider the thickness to be symmetrically distributed relative to the surface
$\hat {z} = 0$
, so the kinematic condition is written at
$\hat {z} = \hat {w} - \hat {h}_i(\hat {y}) /2$
. Small ice thickness allows us to assume that the neutral axis of the plate remains horizontal, which leads to the assumption of symmetry in the ice cross-section. If the cross-section is not symmetric, we can model it as symmetric while maintaining the thickness distribution, as we do not account for changes in the fluid thickness. Additionally, in the linear theory of elasticity, when deriving the vibration equations of the plate, equations of elasticity are integrated across plate’s thickness under the assumption that the midsurface is free from stresses, as described in the Kirchhoff–Love theory.
The equation of the plate describes its vertical displacement,
$\hat {z} = \hat {w}$
, so corresponding velocity vector is given by
$\boldsymbol{v}_{plate} = (0,0,\hat {w}_{\hat {t}})$
. The normal to the surface
$\hat {z} = \hat {w} - \hat {h}_i(\hat {y}) /2$
is
\begin{align} ( - \hat {w}_{\hat {x}}, -\hat {w}_{\hat {y}} + \hat {h}_{i,\hat {y}} /2, 1 ). \end{align}
Thus, in the nonlinear formulation, the kinematic condition (2.8) reads
\begin{align} - \frac {\partial \hat {\varphi }}{\partial \hat {x}} \frac {\partial \hat {w}}{\partial \hat {x}} - \frac {\partial \hat {\varphi }}{\partial \hat {y}} \frac {\partial \hat {w}}{\partial \hat {y}} + \frac {1}{2} \frac {\partial \hat {\varphi }}{\partial \hat {y}} \frac {\partial \hat {h}_i}{\partial \hat {y}} + \frac {\partial \hat {\varphi }}{\partial \hat {z}} = \frac {\partial \hat {w}}{\partial \hat {t}}. \end{align}
In the classic linear theory of hydroelasticity, it is assumed that all terms in this equation are of the same order of smallness, and in the leading-order approximation, the first two terms on the left-hand side of the equation are omitted. We restrict ourselves to the case where the third term is also neglected. This term can be omitted if, for example,
$\hat {h}_{i,\hat {y}}$
is of the same order with
$\hat {w}_{\hat {y}}$
or any other derivative appearing in the equation. In all subsequent results, we have maintained the condition that
$\hat {h}_{i,\hat {y}}$
is of the same order or smaller than
$\max \limits _{\hat {y}} (\hat {w}_{\hat {y}})$
in both dimensional and non-dimensional variables.
The dynamic condition (2.1) and the kinematic condition (2.7) are both given at
$\hat {z} = 0$
. In the original nonlinear formulation, these conditions are given at
$\hat {z} = \hat {w} - \hat {h}_i(\hat {y}) / 2$
. Within the classic linear theory, in the leading order and in the absence of the variable thickness, this condition is rewritten at
$\hat {z} = 0$
. In our case, we solve the problem also in the leading order and under the assumption that
$(\hat {h}_i(\hat {y})) / 2$
is of the same order of
$\hat {w}$
or smaller.
Unknown functions
$\hat {w}$
and
$\hat {\varphi }$
are sought in the form of periodic hydroelastic waves propagating along the channel in the positive
$\hat {x}$
-direction with constant amplitude
$A$
, wavenumber
$k$
and frequency
$\omega$
,
\begin{equation} \hat {w}(\hat {x},\hat {y},\hat {t})= \mathrm{Re}\big( A F(\hat {y}) e^{i(k \hat {x} - \omega \hat {t})} \big), \end{equation}
\begin{equation} \hat {\varphi }(\hat {x},\hat {y},\hat {z},\hat {t})= \mathrm{Re}\big( \!- i \omega A \Phi (\hat {y},\hat {z}) e^{i(k \hat {x} - \omega \hat {t})}\big). \end{equation}
Using this form, we shall determine the dispersion relations
$\omega (k)$
and other characteristics of the periodic waves. The main focus of the study is on the effect of variable ice thickness
$\hat {h}_i(\hat {y})$
on the characteristics of the waves.
The problem is solved in dimensionless variables. Further, all variables and functions are written in dimensionless form. Dimensional variables are marked with
$\hat{}$
– dimensional, dimensionless variables have conventional notations. Dimensional wavenumber is denoted by
$k$
, dimensionless by
$\kappa$
,
$\kappa = k b$
. Length scale is taken to be half the channel width
$b$
, the time scale is
$1 / \omega$
, the amplitude
$ A$
is the scale of the ice deflections,
$\rho _\ell g A$
is the scale of the liquid pressure,
$Ab \omega$
is the scale of the flow velocity potential. We introduce an average value of the ice thickness across the channel,
$h_{avg}$
, and use it as a scale of the function
$\hat {h}_i (\hat {y})$
. Note that the scales of the ice deflections and the flow velocity potential are proportional to the amplitude
$A$
due to the linearity of the problem, and the parameter
$A$
can be arbitrary within conditions of the linear theory. In dimensionless variables, the channel cross-sectional boundaries are (
$ - 1 \lt y \lt 1$
), (
$ - h \lt z \lt 0$
), where
$ h = H / b$
is the dimensionless channel depth. We will consider only solutions with positive and real
$k$
and real
$\hat {F}(\hat {y})$
.
In the dimensionless variables equation of the plate reads
\begin{align} \beta \bigg[h^3_i(y)\left (\frac {\partial ^4{w}}{\partial {x}^4}+ 2\frac {\partial ^4{w}}{\partial {x}^2\partial {y}^2}+ \frac {\partial ^4{w}}{\partial {y}^4}\right) &+ 6 h^2_i(y) \frac {\textrm{d} h_i}{\textrm{d}y} \left (\frac {\partial ^3{w}}{\partial {x}^2\partial {y}} + \frac {\partial ^3{w}}{\partial {y}^3}\right) \nonumber \\ &+ 6 h_i(y) \left (\frac {\textrm{d} h_i}{\textrm{d}y}\right)^2 \left (\frac {\partial ^2{w}}{\partial {y}^2}+ \mu \frac {\partial ^2{w}}{\partial {x}^2}\right) \bigg] + \delta \gamma \frac {\partial ^2{w}}{\partial {t}^2}\nonumber \\ &\quad = {p(x,y,0,t)}, \end{align}
where
$\beta ={D_{avg}}/[{\rho _{\ell }gb^4}]$
,
$D_{avg}= E h^3_{avg} /[12(1-\mu ^2)]$
,
$\gamma ={b\omega ^2}/{g}$
and
$\delta ={h_{sc} \rho _i}/[{b\rho _{\ell }}]$
. The problem is rewritten in terms of
$F(y)$
and
$\Phi (y,z)$
. Substituting (2.11) and (2.12) into (2.13) gives
\begin{align} \beta \bigg[ h^3_i(y) ( \kappa ^4 F - 2\kappa ^2 F^{\prime \prime } + F^{\prime \prime \prime \prime } ) &+ 6h^2_i(y)h_i^{\prime } ( -\kappa ^2 F^{\prime } + F^{\prime \prime \prime } ) \nonumber \\ &+ 6h_i(y)(h_i^{\prime })^2 ( F^{\prime \prime } - \mu \kappa ^2F) - \frac {\delta \gamma }{\beta } h_i(y) F \bigg]\nonumber \\ &\quad = \gamma \Phi (y,0) - F, \end{align}
where prime stands for the derivative with respect to
$y$
. The rest of the equations of the problem in the dimensionless variables are
\begin{equation} F = 0, \quad F_y = 0 \qquad (y = \pm 1), \end{equation}
\begin{equation} \Phi _{yy} + \Phi _{zz} = \kappa ^2 \Phi \qquad (-1 \lt y\lt 1, -h\lt z\lt 0), \end{equation}
\begin{equation} \Phi _y = 0 \quad (y= \pm 1), \qquad \Phi _z = 0 \quad (z = -h), \qquad \Phi _z = F \quad (z = 0). \end{equation}
The solution of the problem (2.14)–(2.17) depends on the three dimensionless parameters
$\beta$
,
$\gamma$
and
$\delta$
. We shall find dispersion relations
$\omega (k)$
for some typical values of parameters of the problem. Then we can calculate phase and group velocities, as well as wave profiles and strain distributions in the ice cover.
In this paper, we use clamped boundary conditions to describe the case where the ice sheet is frozen to the channel walls. These conditions were observed in experiments in an ice tank, the parameters of which we use, see Pogorelova et. al. (Reference Pogorelova, Zemlyak and Kozin2019a ). However, the method can be adapted to other boundary conditions that may be required for specific practical problems. The problem with constant ice thickness and clamped conditions was studied by Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014). Same problem within another boundary conditions (free–free conditions) was investigated in Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2015). The latter study provided a comparison between both cases, showing that, for ice with free–free conditions, periodic waves generally propagate along the channel at lower speeds and with lower frequencies.
In our paper we further consider two different types of ice thickness variation across the channel. In the first case, the ice thickness is maximal at the walls and has its minimum value near the centre of the channel. The second case describes to the opposite distribution. As will be shown further in § 4, for the first type of ice thickness shape, the change in wave characteristics, such as dispersion relations, phase and group velocity, is monotonic and predictable as the rate of thickness variation increases. In general, all characteristics decrease. No abnormal behaviour was observed. Based on our experience, we may expect a similar reduction in frequencies and wave speeds for free–free boundary conditions within the first case. However, this still requires verification. In the second case, the change in wave characteristics is non-monotonic, and the conclusions of Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2015) cannot be used for a preliminary analysis. Other boundary conditions, such as simply supported conditions, can also be considered. These conditions may also have specific physical applicability. From a mathematical point of view, other boundary conditions do not introduce any new significant complexities into the solution of the problem, and the solution method can be repeated with minor adjustments.
3. Method of solution
The solution of (2.14) is sought in the form of an infinite series of eigenmodes of oscillations of an elastic beam with variable thickness
\begin{equation} F(y)=\sum _{n=1}^{\infty }a_n\psi _n(y). \end{equation}
These spectral functions are introduced for each case of variable thickness; however, in general, the algorithm for determining these functions is the same and consists of repeated steps. To determine the functions
$\psi _n(y)$
, the channel width,
$y \in [-1, 1]$
, is divided into small segments
$[b_k, b_{k+1}]$
, where
$1 = b_0 \lt b_1 \lt \ldots \lt b_{N_k-1} \lt b_{N_k} = 1$
. On each of these segments, the ice thickness
$h_i(y)$
is approximated by a linear function
$h^{(k)}_i (y)$
, see figure 2. In total there can be three different cases of thickness variation shown in figure 2(a–c).
Example of thickness shape approximations. (a)–(c) Possible ice thickness shapes for its piecewise approximation by linear functions.
On the segments with constant thickness (see figure 2a
) the function
$h^{k}_i(y)$
is equal to
$h^{(k)}_{avg}$
. Note that any of these values are not equal to
$h_{avg}$
. The equality is possible only in the case of a plate with constant thickness over a channel. The functions
$\psi _n$
are solutions of the spectral problem describing the oscillations of an elastic beam with a constant thickness
\begin{equation} \left (h^{(k)}_{avg}\right)^2 \psi ^{(k)}_{n,yyyy}=\theta _n^4 \psi ^{(k)}_n(y) \quad (b_k \leqslant y \leqslant b_{k+1}). \end{equation}
The solution of this problem is
\begin{align} \psi ^{(k)}_n(y) &= A^{(k)}_n \cos (\theta ^*_n y) + B^{(k)}_n \sin (\theta ^*_n y) + C^{(k)}_n \cosh (\theta ^*_n y) + D^{(k)}_n \sinh (\theta ^*_n y),\nonumber \\ \theta ^*_n & = \frac {\theta _n}{\sqrt {h^{(k)}_{avg}}}. \end{align}
The expansion coefficients
$A^{(k)}_n, B^{(k)}_n, C^{(k)}_n$
and
$D^{(k)}_n$
, as well as the spectral parameter
$\theta _n$
, which is the same for all solutions on each segment, are determined from the orthogonality condition for these modes and the boundary conditions at the boundary of the segments. These conditions are specified further. For the case with linearly varying ice thickness (see figure 2b
,c), the values of the thickness on the left-hand,
$h^{(k)}_{left}$
, on the right-hand,
$h^{(k)}_{right}$
, and the average value of the thickness on this segment,
$h^{(k)}_{avg}$
, are introduced. Then the approximated ice thickness
$\tilde {h}_i(y)$
is described by the function
\begin{equation} \tilde {h}_i (y) =\alpha ^{(k)}_0 + \alpha ^{(k)} y, \quad \tilde {h}_i(b_k) = h^{(k)}_{left}, \quad \tilde {h}_i(b_{k+1}) = h^{(k)}_{right}, \end{equation}
where
$\alpha ^{(k)}_0 = h^{(k)}_{right} - \alpha ^{(k)} b_{k+1}$
,
$\alpha ^{(k)} = (h^{(k)}_{right} - h^{(k)}_{left})/(b_{k+1} - b_k)$
. The most important parameter is
$\alpha ^{(k)}$
which describes the tangent of the angle’s inclination for the linear variation in ice thickness at given segment. Note that these angles are for the shape of the dimensionless ice thickness. Dimensional angles
${\zeta }^{(k)}$
can be calculated by formula
$\arctan (\alpha ^{(k)} h_{avg} /(2 b))$
. These angles are small to keep correctness of the linear model and usability of the boundary kinematic and dynamic conditions, see discussion after (2.7). Functions
$\psi _n(y)$
are solutions of the spectral problem describing the oscillations of an elastic beam with a linearly varying thickness
\begin{equation} \!\!\tilde {h}^2_i(y)\psi ^{(k)}_{n,yyyy} + 6\tilde {h}_i(y)\tilde {h}_{i,y}\psi ^{(k)}_{n,yyy} + 6(\tilde {h}_{i,y})^2 \psi ^{(k)}_{n,yy} =\theta _n^4\psi {(k)}_n \quad (b_k \leqslant y \leqslant b_{k+1}).\,\, \end{equation}
The solution of this problem is
\begin{align} \psi ^{(k)}_n(y) & =\frac {1}{\xi }\big[A^{(k)}_n\,J_1(\eta ^{(k)}_n\xi)+B^{(k)}_n\,Y_1(\eta ^{(k)}_n\xi)+C^{(k)}_n\,I_1(\eta ^{(k)}_n\xi)+D^{(k)}_n\,K_1(\eta ^{(k)}_n\xi)\big] \nonumber \\ & \quad n=1,2,3\ldots, \end{align}
where
$\eta ^{(k)}_n={2\theta _n}/{\alpha ^{(k)}}$
,
$\xi =\sqrt {\alpha _0^{(k)}+\alpha ^{(k)} y}$
and
$J$
,
$Y$
,
$I$
,
$K$
are Bessel functions. The same spectral problem but for linear asymmetric and symmetric cases of thickness distributions was considered in Batyaev and Khabakhpasheva (Reference Batyaev and Khabakhpasheva2022)and Shishmarev et al. (Reference Shishmarev, Zavyalova, Batyaev and Khabakhpasheva2022). A similar spectral problem for a free-free beam was studied in Korobkin et al. (Reference Korobkin, Khabakhpasheva and Shishmarev2023).
When the solutions on each segment are determined, the resulting function
$\psi _n(y)$
is introduced as
\begin{align} \psi _n(y)=\left \{ \begin{array}{ll} \psi ^{(1)}_n(y), \quad (-1 \lt y \leqslant b_1),\\ \ldots \\ \psi ^{(k)}_n(y), \quad (-b_k \leqslant y \leqslant b_{k+1}),\\ \ldots \\ \psi ^{(N_k)}_n(y), \quad (b_{N_k - 1} \leqslant y \lt 1).\\ \end{array} \right. \end{align}
We shall determine the boundary conditions at the channel walls and between each segment
$[b_k, b_{k+1}]$
. Since we are considering ice frozen to the channel walls, it gives
\begin{equation} \psi ^{(1)}_n(-1) = 0, \quad \psi ^{(1)}_{n,y}(-1) = 0, \quad \psi ^{(N_k)}_n(1) = 0, \quad \psi ^{(N_k)}_{n,y}(1) = 0. \end{equation}
At the boundaries of the segments,
$y = b_{k}$
,
$k = 1, 2, \ldots, N_k - 1$
, it is necessary that the continuity conditions of ice deflections, slope, bending moments and shear forces are satisfied. Taking into account the same average ice rigidity
$D_{avg}$
on every segment, these conditions read
\begin{equation} \psi ^{(k)}_n(b_k) = \psi ^{(k+1)}_n(b_k), \end{equation}
\begin{equation} \psi ^{(k)}_{n,y}(b_k) = \psi ^{(k+1)}_{n,y}(b_k), \end{equation}
\begin{equation} \psi ^{(k)}_{n,yy}(b_k) = \psi ^{(k+1)}_{n,yy}(b_k), \end{equation}
\begin{equation} \mbox{shear force}^{(k)}(b_k) = \mbox{shear force}^{(k+1)}(b_k), \end{equation}
where expression
$\mbox{shear force}^{(k)}(y)$
in the last condition equals to
$\psi ^{(k)}_{n,yyy} + 3 \tilde {h}_{i,y} \psi ^{(k)}_{ n,yy}(y)$
or to
$\psi ^{(k)}_{n,yyy}(y)$
depending on the case, linear or constant, for the corresponding segment. In total, these conditions add up to
$4 N_k$
algebraic equations. The last needed equation is the orthogonality condition of the modes
$\psi _n(y)$
with a weight on each segment
\begin{align} \int \limits _{-1}^{1}\tilde {h}_i(y)\psi _n(y)\psi _m(y)\textrm{d}y=\delta _{nm} \qquad (n\neq m). \end{align}
The kinematic condition (2.17) shows that the profile
$\Phi (y, z)$
of the flow potential in the cross-section of the channel can be sought in the form
\begin{equation} \Phi (y,z)=\sum _{n=1}^{\infty }a_n\Phi _n(y,z) \end{equation}
and
$\Phi _n$
are the solutions of the boundary value problems
\begin{equation} \frac {\partial ^2\Phi _n}{\partial y^2}+ \frac {\partial ^2\Phi _n}{\partial z^2}=\kappa ^2\Phi _n, \end{equation}
\begin{align} \frac {\partial \Phi _n}{\partial y}=0\quad (y=\pm 1),\qquad \frac {\partial \Phi _n}{\partial z}=0\quad (z=-h),\qquad \frac {\partial \Phi _n}{\partial z}=\psi _n\quad (z=0). \end{align}
The solution of the last problem is
\begin{align} \Phi _n(y,z)=\frac {1}{2}\Phi _{n0}^c(z)+\sum _{k=1}^{\infty }\Phi _{nk}^c(z)\cos (\pi k y)+ \sum _{l=1}^{\infty }\Phi _{nl}^s(z)\sin [(\pi l-\pi /2) y], \end{align}
where
\begin{equation} \Phi _{nk}^c=\frac {\psi _{nk}^c}{\mu _k^c}\frac {\cosh \mu _k^c(z+h)}{\sinh \mu _k^c h},\qquad \qquad \Phi _{nl}^s=\frac {\psi _{nl}^s}{\mu _l^s}\frac {\cosh \mu _l^s(z+h)}{\sinh \mu _l^s h}, \end{equation}
\begin{equation} \psi _{nk}^c=\int \limits _{-1}^1\psi _n(y)\cos (\pi k y)\textrm{d}y,\qquad \psi _{nl}^s=\int \limits _{-1}^1\psi _n(y)\sin [(\pi l-\pi /2) y]\textrm{d}y, \end{equation}
\begin{equation} \mu _k^c=\sqrt {\kappa ^2+(\pi k)^2}\qquad \qquad \mu _l^s=\sqrt {\kappa ^2+(\pi l-\pi /2)^2}. \end{equation}
Integrals
$\Phi _{nk}^c$
and
$\Phi _{nl}^s$
are calculated numerically.
Substituting (3.1) and (3.14) into (2.14), we arrive at the infinite system of equations to find the principal coordinates
$a_n$
,
\begin{align} \sum _{n=1}^{\infty }a_n\left \{\beta \tilde {h}_i\theta _n^4\psi _n + \kappa ^4\beta \tilde {h}^3_i\psi _n - 2\kappa ^2\beta \left ( \tilde {h}^3_i \psi _n^{\prime \prime } + 3\tilde {h}_{i,y} \tilde {h}^2_i\psi _n^{\prime } \right) \right. \nonumber \\ \left. - 6 \kappa ^2 \mu \beta (\tilde {h}_{i,y})^2 \tilde {h}_i\psi _n - \delta \gamma \tilde {h}_i \psi _n - \gamma \Phi _n(y,0) + \psi _n \right \} = 0. \end{align}
Note that
$\tilde {h}_{i,y}(y)$
is equal to
$\alpha ^{(k)}$
depending on
$y$
and, in general, is piecewise given step function with discontinuities of the first kind at the boundaries between segments
$[b_k, b_{k+1}]$
. Multiplying both sides of (3.21) by
$\psi _m(y)$
, integrating the result from
$-1$
to
$1$
in
$y$
and reducing the number of the equations to a finite value
$N$
, we arrive at the following algebraic problem:
\begin{equation} \sum \limits _{n=1}^N a_n \left \{\beta D_n \delta _{nm} - 6 \beta \mu \kappa ^2 T_{nm} + \kappa ^4 \beta K_{nm} + 2 \kappa ^2 \beta S_{nm} - \delta \gamma \delta _{nm} + M^{(1)}_{nm} - \gamma M^{(2)}_{nm}\right \} = 0, \end{equation}
where
\begin{align} D_n = \theta _n^4, \quad K_{nm} = \int \limits _{-1}^1 \tilde {h}_i^3 \psi _n \psi _m \textrm{d}y, \quad S_{nm} = \int \limits ^1_{-1} \tilde {h}_i^3 \psi _n^{\prime }\psi _m^{\prime }\textrm{d}y, \end{align}
\begin{align} T_{nm} = \int \limits ^1_{-1} (\tilde {h}_i^\prime)^2 \tilde {h}_i \psi _n\psi _m\textrm{d}y, \quad M^{(1)}_{nm} =\int \limits ^1_{-1} \psi _n\psi _m\textrm{d}y, \quad M^{(2)}_{nm}= \int \limits ^1_{-1} \Phi _n(y,0)\psi _m \textrm{d}y. \end{align}
The resulting system of linear equations for the principal coordinates
$a_n$
in the series (3.1) can be written in matrix form,
\begin{equation} \{\beta [\boldsymbol {D} - 6 \mu \kappa ^2 \boldsymbol {T}+ \kappa ^4\boldsymbol {K} +2\kappa ^2\boldsymbol {S}] +\boldsymbol{M}_{\textbf{1}} - \gamma [\delta \boldsymbol {I}+\boldsymbol{M}_{\textbf{2}}]\}\boldsymbol {a} = \boldsymbol {R}\boldsymbol {a} = 0, \end{equation}
where all matrices
$\boldsymbol {D} = diag\{D_{n}\}$
,
$\boldsymbol {T} = \{T_{nm}\}$
,
$\boldsymbol {K} = \{K_{nm}\}$
,
$\boldsymbol {S} = \{S_{nm}\}$
,
$\boldsymbol{M}_{\textbf{1}}=\{M_{nm}^{(1)}\}$
,
$\boldsymbol{M}_{\textbf{2}}=\{M_{nm}^{(2)}\}$
are symmetrical,
$\boldsymbol {I}$
is an identity matrix. Non-trivial solutions of the resulting system of equations exist if the determinant of the matrix in curly brackets in (3.22) is equal to 0. Frequencies
$\omega _n$
are calculated for a given interval of wavenumbers
$\kappa$
. In general, for unreduced system of equations (3.22) there will be infinitely many such relations,
$\omega _n(\kappa)$
. In the reduced case, we can calculate only first
$N$
frequencies. Dispersion relations are numbered being the smallest for
$n=1$
and increasing with the number
$n$
for the same wavenumber
$\kappa$
. Each calculated dispersion relation
$\omega _n(k)$
gives unique solution of the problem in the form (2.11)–(2.12) with unique wave profile
$F_n(y)$
across the channel.
In addition to the wave characteristics of the studied waves, we will also investigate the distribution of stresses and strains in the ice cover. In the linear theory, the stress is proportional to the strain and varies linearly through the thickness of the ice, being zero at the middle surface (Timoshenko Reference Timoshenko1951; Timoshenko and Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1959). For this reason we restrict ourselves to the study of the strains only.
The solution to the problem is sought in the form of waves propagating along the channel. We consider periodic waves, so the wave is expressed through its phase
$\theta (x,t) = \kappa x - t$
,
$0 \leqslant \theta (x,t) \lt 2 \pi$
,
\begin{align} w(x,y,t) = F(y) \cos (\theta). \end{align}
The absolute maximum value of the dimensionless strain along the phase,
$\varepsilon _{{ABS}}$
, is calculated as
\begin{align} \varepsilon _{{ABS}} = \max _{-1 \leqslant y \leqslant 1} \left [\max _{0 \leqslant \theta \lt 2\pi } \varepsilon (y, \theta)\right], \end{align}
where
$\varepsilon (y, \theta)$
is the strain at point
$y$
and phase
$\theta$
. Note that the distribution of these strains is repeated periodically along the channel. The strain field in the plate is described by the strain tensor
\begin{align} T_{\varepsilon } = - \zeta \begin{bmatrix} w_{xx} & w_{xy} \\ w_{xy} & w_{yy} \end{bmatrix}, \end{align}
where
$\zeta$
is the coordinate through the ice plate thickness,
$\zeta \in [(-{\tilde {h}_i(y)})/{2}, ({\tilde {h}_i(y)})/{2}]$
.
To compute the maximum strain at the point
$(\theta, y)$
, we need to find the principal strains, which are determined from the eigenvalues of the strain tensor
$T_{\varepsilon }$
. The principal strains are given by
\begin{align} \varepsilon ^{\pm }(y,\theta) = \frac {\zeta }{2} \left [-a \cos {\theta } \pm \sqrt {(b^2 - c^2)\cos ^2{\theta } + c^2}\right], \end{align}
where
\begin{align} a(y) = F^{\prime\prime}(y) - \kappa F(y), \quad b(y) = F^{\prime\prime}(y) + \kappa ^2 F(y), \quad c(y) = 2\kappa F^{\prime}(y). \end{align}
Here,
$a$
,
$b$
and
$c$
are functions depending on the wave profile
$F(y)$
and the wavenumber
$\kappa$
.
The maximum strains occur at the ice plate’s surface, where
$\zeta = \pm ({\tilde {h}_i(y)})/{2}$
. Since the functions
$\varepsilon ^{\pm }(y, \theta)$
are smooth in
$\theta$
, their maximum is reached at the extrema, which are found from the condition
\begin{align} \frac {\textrm{d} \varepsilon ^{\pm }}{\textrm{d} \theta } = 0. \end{align}
The maximum principal strain values along the wavelength are determined as follows:
\begin{align} \varepsilon _{{max}} = \begin{cases} \varepsilon ^{(1)}_{{max}}, & \text{if } c^2 - b^2 \leqslant |ab|, \\[5pt] \varepsilon ^{(2)}_{{max}}, & \text{if } c^2 - b^2 \geqslant |ab|, \end{cases} \end{align}
where
\begin{align} \varepsilon ^{(1)}_{{max}} = \frac {\tilde {h}_i(y)}{2} (|a| + |b|), \quad \varepsilon ^{(2)}_{{max}} = \frac {\tilde {h}_i(y) |c|}{2} \sqrt {\frac {b^2 - c^2 - a^2}{b^2 - c^2}}. \end{align}
In the case of
$\varepsilon ^{(1)}_{{max}}$
, the maximum strain is either transverse (across a channel),
$w_{yy} = F^{\prime\prime}$
, or longitudinal (along a channel),
$w_{xx} = \kappa ^2 F$
. In the case of
$\varepsilon ^{(2)}_{{max}}$
, it corresponds to combined maximum strains that are neither strictly transverse nor longitudinal, and they form an angle with the channel axis.
4. Numerical results and discussion
Calculations of periodic hydroelastic waves and their characteristics were carried out for parameters of the channel corresponding to the experimental ice tank at the Sholem Aleichem Amur State University in Birobidzhan:
$ H = 1$
m,
$ 2b = 3$
m. Average ice thickness is assumed to be 0.0035 m, which was used in most of the published experiments in the same pool (see Zemlyak et al., (Reference Zemlyak, Kozin, Baurin and Petrosyan2014, Reference Zemlyak, Baurin, Kozin, Chizhiumov and Matiushina2019, Reference Zeng, Korobkin, Ni and Xue2022), among others). The density and parameters of rigidity of ice are
$\rho _i = 917$
kg m
$^{-3}$
,
$\mu = 0.3$
, and
$E = 4.2 \times 10^9$
Pa. These values were used in our previous researches, see, e.g. Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014) and Shishmarev et al., (Reference Shishmarev, Khabakhpasheva and Oglezneva2023, Reference Shishmarev, Sibiryakova, Naydenova and Khabakhpasheva2024). The Young’s modulus used here is within the range of values reported by various researchers, see table 3.1 in Liang et al. (Reference Liang, Pitt and Bennetts2024).
The algorithm of calculations has a few steps: first, the values of
$\theta _n$
are calculated from the boundary conditions (3.8)–(3.12); then, the functions
$\psi _n$
are determined and the dispersion relations
$\omega _n(k)$
are calculated as values for which the determinant of the matrix in (3.22) is equal to 0. Hydroelastic waves with wavelengths
$\lambda \gt 0.42$
m,
$0 \lt k \lt 15$
m
$^{- 1}$
, are considered. Frequencies
$\omega _n(k)$
are determined for finite number of wavenumbers
$k$
from this interval with a constant step. The profiles of the hydroelastic modes
$ F_n (y)$
are calculated using (3.1) with the coefficients
$ a_n$
calculated as eigenvectors of the matrix equation (3.22), and normalized so that
$\max {|F_n(y)|} = 1$
.
To compute
$\theta _n$
, it is necessary to calculate the determinant of the diagonally structured matrix given by (3.8)–(3.12). This parameter describes piecewise normal vibrations of a beam for a given piecewise ice thickness in the channel. Each
$\theta _n$
corresponds to a unique vibration profile and depends on the plate parameters, including the variable ice thickness. In the case of a piecewise-linear approximation with three segments, the system (3.8)–(3.12) consists of 12 algebraic equations. The piecewise-defined modes
$\psi _n$
are constructed to have continuous derivatives up to the third order (the third derivative is ensured to be continuous by the continuity condition of shear forces). These conditions prevent any singularities at the segment boundaries for the considered derivative orders.
The parameters
$\theta _n$
increase with
$n$
, which limits the number of computable modes
$\psi _n(y)$
for small variations in plate thickness. For an extremely small angle
$\alpha ^{(k)}$
in a segment with linearly varying piecewise approximated thickness
$\tilde {h}_i$
, the computation of modes becomes impossible because the modified Bessel function
$I_0$
depends on
$\eta ^{(k)}_n = {2\theta _n}/{\alpha ^{(k)}}$
, which diverges as
$\alpha ^{(k)} \to 0$
. The main difficulty in computing the eigenvalues of the matrix (3.8)–(3.12) arises from the third derivatives of the Bessel functions
$I$
and
$K$
at the segment boundaries
$b_k$
. Each derivative is evaluated at the left- and right-hand boundaries of a segment, which can lead to excessively large or small values. To mitigate these issues, the functions
$I$
and
$K$
are scaled as
$I \sim e^{{arg}} / \sqrt {2\pi \; \text{arg}}$
and
$K \sim e^{-{arg}} \sqrt {2 / (\pi \; \text{arg})}$
, where
$\text{arg}$
is the maximum argument of the Bessel functions at the segment boundaries. This method enables an accurate computation of up to 40–50 modes, depending on the piecewise approximated ice thickness profile. Note that each new segment gives four new conditions and correspondingly four new equations and new set of Bessel functions
$I$
and
$K$
, which values can be large. One needs to be careful with a large number of segments in order to ensure the accuracy of the calculations.
The dispersion relations
$\omega _n(k)$
are computed as the values at which the determinant of the matrix
$\boldsymbol{R}$
in (3.25) changes sign. These values are determined with an accuracy of
$10^{-12}$
. The advantage of the constructed piecewise-defined normal modes
$\psi _n(y)$
lies in their orthogonality. All submatrices in
$\boldsymbol{R}$
are symmetric, which simplifies the computation of eigenvalues and eigenvectors. All coefficients and matrices in (3.25) are real. For the undamped ice plate model (where damping due to viscosity or ice porosity is neglected), we did not find any complex frequencies
$\omega _n$
for the corresponding real-valued wavenumbers
$\kappa$
. However, when damping is taken into account, complex frequencies may exist even for real wavenumbers. It should be noted that classic periodic waves propagating along the channel with a constant amplitude cannot exist in an ice cover with damping. The accuracy of computations
$\omega _n(k)$
also increases with the number of modes considered and is determined by the absolute difference in frequencies calculated for different numbers of modes. Typically, the first four smallest frequencies could be computed with an accuracy of 0.1 % considering only 10 normal modes
$\psi _n(y)$
for the taken range of wavenumbers. For frequencies with larger indices, a greater number of considered modes
$\psi _n(y)$
is required. Generally, we used 30–40 modes in presented calculations. Note that the elements
$R_{nn}$
of the final matrix in (3.25) behave like
$\theta _n^4$
, where
$\theta _n \gt n$
. To avoid computational problems with large values, a new vector
$a_n^* = \theta _n^2 a_n$
can be introduced, and the
$m$
th equation in the system (3.22) can be divided by
$\theta _m^2$
. The resulting matrix
$\boldsymbol {R*}$
with elements
$R_{nm}/(\theta _n^2 \theta _m^2)$
is symmetric, and all its elements are bounded as
$n,m \to +\infty$
. For the plate with constant thickness, analysis of the values
$a_n$
showed that the main contribution to each wave with frequency
$\omega _n$
was made by a classic normal mode with the same index
$n$
. The contribution of other modes quickly diminishes as one moves away from this index, see Korobkin et al. (Reference Korobkin, Khabakhpasheva and Papin2014). These contributions are significantly different for plates of variable thickness.
The results are presented for two different types of ice plate’s thickness shapes. These shapes are shown in figure 3. For these shapes simple piecewise-approximated forms will be provided. For each piecewise-approximated form we consider two additional piecewise forms with large variation in thickness. In the first case, the ice thickness has large values near the walls of the channel and a small value at the middle segment of the channel. In the second case, the opposite distribution of thickness is considered, with small values near the walls and large value closer to the centre of the channel. The asymmetric variation in ice thickness is studied. The main goal of the study is to analyse the influence of approximated ice plate’s thickness shapes on the behaviour and characteristics of periodic hydroelastic waves propagating along the channel. The thickness of naturally frozen ice varies due to multiple factors, primarily temperature. Lower temperatures result in more uniform ice formation, while higher temperatures lead to greater thickness variation. Near walls regions tend to be thicker than the centre, especially at moderate freezing temperatures. The figure 3(a) shape was chosen to model that scenario. Note that obtaining experimental data about the thickeness of natural ice and its correlations remains challenging.
As mentioned previously, there will be a countable number of such waves, each with a unique frequency
$\omega _n(k)$
and a unique wave profile across the channel
$F_n(y)$
. It is important to note that the wave profile
$F_n(y)$
also depends on the wavenumber
$k$
. For convenience, the dispersion relations are numbered in ascending order, i.e. the wave with frequency
$\omega _1(k)$
is the wave with the lowest frequency in the considered channel. We will refer to each such set of waves with a single index for different wavenumbers as a ‘waveform’. Accordingly, the first waveform corresponds to the waves with frequency
$\omega _1(k)$
and wave profile
$F_1(y)$
. In further sections, we consider piecewise given ice thickness only, so word ‘thickness’ will usually refer to its piecewise approximation
$\tilde {h}_i(y)$
.
Types of ice thickness shapes which piecewise-linear approximations used in the calculations.
4.1. Ice thickness being large near the channel walls and small near the centre of the channel
In this section we considered the shape approximating the original shape shown in figure 3(a). The shapes of the dimensionless piecewise-approximated ice thickness across the channel used in the calculations are shown in figure 4. For all shapes, the channel width is divided into three segments. The left- and right-hand segments are segments with linear thickness, while the middle segment has a constant ice thickness. All segments have its average thickness,
$h_{avg}^{(i)}$
,
$i = 1,2,3$
, and linear segments have its angles of linear inclination,
$\alpha _{1,3}$
. The shape given in figure 3(a) corresponds to the piecewise approximation of the original thickness. Shapes given by figure 3(a–c) are similar to each other, differing in the ratio of the average ice thickness on the left- and right-hand linear segments: for case figure 3(a),
$h_{avg}^{(1)} /h_{avg}^{(3)} = 1$
, for figure 3(b),
$h_{avg}^{(1)} /h_{avg}^{(3)} = 1.5$
and for figure 3(c),
$h_{avg}^{(1)} /h_{avg}^{(3)} = 2$
. By changing the ratio of the average ice thickness on the middle segment to the average ice thickness on the right-hand (or left-hand) linear segment, we control the rate of variation in ice thickness across the channel and could set it to be nearly constant or highly variable. Note that nearly constant ice thickness can only be achieved in the figure 3(a) case due to the relation between the average thicknesses on the left- and right-hand linear segments. In all calculations, the average ice thickness in the ice cover remains unchanged and independent of the shape, and in dimensionless variables, it is equal to unity. The dimensionless ice thickness is scaled by the value
$h_{avg}=0.0035$
m. Therefore, the calculation results will be compared with results for an ice plate of constant thickness
$0.0035$
m, this shape is illustrated by figure 4(d). The segment boundaries are
$-0.7$
and
$0.3$
in dimensionless variables and remain unchanged in the calculations. The dimensionless functions
$\tilde {h}_i(y)$
for each shape are shown in figure 5. The values of average thicknesses on each segments are: for the blue shape
$h^{(1)}_{avg} = h^{(3)}_{avg} \approx 1.42$
,
$h^{(2)}_{avg} = 0.57$
; for the red shape
$h^{(1)}_{avg} \approx 2.14$
,
$h^{(2)}_{avg} = 0.36$
,
$h^{(3)}_{avg} \approx 1.42$
; for the green shape
$h^{(1)}_{avg} \approx 2.8$
,
$h^{(2)}_{avg} = 0.14$
,
$h^{(3)}_{avg} \approx 1.42$
.
The piecewise ice thickness shapes across the channel used in the calculations.
Dimensionless piecewise ice thickness distribution
$\tilde {h}_i(y)$
for each case from figure 4 for maximum considered variation in the thickness.
4.1.1. Convergence while decreasing the rate of thickness variation
In this section, we assess the method’s accuracy by comparing the wave characteristics with known values for constant ice thickness. We gradually reduce the variation in ice thickness, starting from an already small variation. In all cases, the average ice thickness remains constant. The convergence of results (difference of values for constant and variable thickness ice) as the thickness variation decreases is demonstrated for the ice thickness profile from figure 3(a), shown in figure 4(a). Figure 6(a) displays the ratios
$ [|\overline {\omega }_n(k) - \omega _n(k)|/\overline {\omega }_n(k) \times 100]$
for the first four waveforms with frequencies from
$\omega _1(k)$
to
$\omega _4(k)$
, as functions of the wavenumber
$k$
. Here,
$\overline {\omega }_n(k)$
represents the frequencies for an ice cover of constant thickness. Figure 6(b) shows the corresponding wave profiles
$F_1(y)$
to
$F_4(y)$
across the channel for
$k = 10.45\, \textrm{m}^{-1}$
. The minimum considered average thickness on the right-hand linear segment,
$0.003541$
m, allows computing at least 10 modes
$\psi _n(y)$
. In the blue case (see legend in figure 6),
$h_{avg}^{(3)} / h_{avg}^{(2)} \approx 1.024$
, in the red case,
$h_{avg}^{(3)} / h_{avg}^{(2)} \approx 1.035$
, and in the yellow case,
$h_{avg}^{(3)} / h_{avg}^{(2)} \approx 1.047$
. The case
$h_{avg}^{(3)} / h_{avg}^{(2)} = 1$
corresponds to constant ice thickness (black lines in figure 6
b).
The differences in dispersion relations depend on both the index and the wavenumber. The largest difference is observed for the wave with the first number and lowest frequency,
$\omega _1(k)$
. For the considered range of wavenumbers, the largest difference occurs for short waves in the yellow case, but it does not exceed 3% and decreases as the thickness variation reduces. For long waves, the ratios
$|\overline {\omega }_n(k) - \omega _n(k)|/\overline {\omega }_n(k)$
are nearly independent of the wavenumber
$k$
. These ratios decrease from
$k = 0$
to a point
$k_*$
where the frequencies match, then increase with
$k$
. The value of
$k_*$
rises with the waveform number
$n$
. Near
$k = 0$
, these ratios are not exceeding 1 %, and decrease as thickness variation decreases. For
$k \gt k_*$
, the ratios
$|\overline {\omega }_n(k) - \omega _n(k)|/\overline {\omega }_n(k)$
decrease with increasing wavenumber
$n$
and can become smaller than for long waves,
$k \lt k_*$
.
Wave profiles in figure 6(b) confirm that waves with lower frequencies are more influenced by variable ice thickness. The largest visual difference is seen for the wave with
$\omega _1(k)$
. These visual differences in wave profiles decrease as thickness variation decreases and with increasing waveform number.
(a) The ratios
$(|\overline {\omega }_n(k) - \omega _n(k)|/\overline {\omega }_n(k) \times 100)$
,
$n = 1,2,3,4$
. (b) The corresponding wave profiles across the channel for
$k = 10.45\,\textrm{m}^{-1}$
. Here
$\overline {\omega }_n(k)$
are for an ice cover with constant thickness,
$\omega _n(k)$
– for the shapein panel (a), see figure 4.
4.1.2. Dispersion relations and corresponding characteristics
The effect of variable thickness on wave characteristics, i.e. dispersion relations, phase velocity, and group velocity, for the first four waveforms is shown in figures 7 and 8. Each colour in these figures corresponds to a specific ice thickness shape shown in figure 4. The dispersion relations,
$\omega _n(k)$
, are presented in figure 7.
Dispersion relations
$\omega _n(k)$
of the first four waveforms. Frequencies of the even waveforms are shown by dot–dashed lines, of the odd waveforms – by solid lines. See figure 4 for the colour–shape relations.
Phase
$C_n(k)$
and group
$C_n^g(k)$
velocities of the first four waveforms. Phase velocities are shown by solid lines, group velocities – by dotted lines. Panels (a), (b), (c) and (d) correspond to the first, second, third and fourth waveforms, respectively.
For long waves, the influence of variable thickness on frequencies is nearly imperceptible. However, for short waves, the effect becomes more noticeable as the wavenumber increases.
Figure 6(a) shows the calculated point
$k_*$
in
$k$
axis, where the frequency difference between constant and variable thickness plates equals zero. Such point also exists for current results. For
$k \lt k_*$
, the frequencies for variable thickness are slightly higher than for constant thickness plates. For
$k \gt k_*$
, the frequencies are lower and the differences grow as the wavenumber increases.
As the variation in thickness increases, the frequencies for the variable thickness plate decrease, as seen in the comparison of the blue, red and green lines. The effect of variable thickness is most noticeable for the green shape, where frequencies for the first and second waveforms, and the third and fourth waveforms, become nearly identical for short waves.
Phase velocities,
$C_n(k) = \omega _n/k$
, and group velocities,
$C_n^{(g)} = \textrm{d}\omega _n /\textrm{d}k$
, for the first four waveforms are shown in figure 8. These velocities follow the same trend as the frequencies. For
$k \lt k_*$
, phase velocities are slightly higher for the variable thickness plate. For
$k \gt k_*$
, these velocities decrease rapidly compared with the constant thickness plate.
This pattern holds for all waveforms. The phase and group velocities intersect at the point of minimal phase velocity. For variable thickness, this occurs at shorter waves and lower velocities. As with dispersion relations, the effect becomes more pronounced with greater variation in ice thickness. For very short waves, the phase and group velocities for the variable thickness plate are much lower than for the constant thickness plate. This effect strengthens as the waveform number increases.
4.1.3. Wave profiles and strains for long waves
The wave profile and strain distribution in the ice cover depend on the wavenumber
$k$
. As shown in figures 7 and 8, wave characteristics of long waves are almost unaffected by the ice thickness shape.
Figure 9 shows the profiles
$F_n(y)$
of long waves for the first four waveforms. Considered ‘long’ waves have
$k = 0.4\,\textrm{m}^{-1}$
, corresponding to a wavelength
$\lambda = 15.7$
m, approximately five times the channel width
$2b = 3$
m. The wave profiles in the variable thickness channel generally follow those of constant thickness ice, with a slight shift of oscillations into the constant thickness area. This shift is most noticeable in the left-hand segment with linear thickness,
$ y \in [-1, -0.7]$
, where the largest slope occurs. At the boundary between the left-hand and middle segments (
$y = -0.7$
), the wave profile curvature increases, especially for the green shape.
Wave profiles
$F_n(y)$
across the channel of the long waves,
$k = 0.4\,\textrm{m}^{-1}$
, for the first four waveforms. See figure 5 for the colour–shape relations. Grey area corresponds to the middle segment with thinnest and constant part of the ice cover.
Dimensionless maximum strains along the phase,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, for the long waves,
$k = 0.4\,\textrm{m}^{-1}$
, for the first four waveforms. See legends for strain direction.
Figure 10 shows the dimensionless maximum strains along the phase,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, corresponding to the wave profiles in figure 9. For constant thickness, the maximum strains for long waves occur at the channel walls and are transverse. These strains increase with the waveform number due to more oscillation peaks and greater profile curvature. For variable thickness shapes, strains at the channel walls behave naturally, decreasing as thickness variation increases and makes ice thicker at the walls. In the middle segment (the grey area), where the ice thins, strains do not increase. For moderate thickness variation (blue and red shapes), strains decrease across the channel, with the most significant reduction at the walls. Local peaks of maximum strains shift towards the centre of the channel as the waveform number increases.
At the left-hand boundary of the segments (
$y = -0.7$
), the maximum strain increases significantly for the red and green shapes. For the green shape, these strains are largest across the entire phase of a wave. These strains are nearly twice those at the walls for the constant thickness ice in the first waveform, 1.5 times larger for the second and third waveforms, and slightly more for the fourth waveform. The maximum strain direction remains unchanged.
Summarizing, while the variable thickness shape has little effect on wave characteristics and weakly influences wave profiles across the channel, it significantly affects the distribution of maximum strains in the ice cover.
4.1.4. Wave profiles and strains for short waves
Figure 11 shows the profiles
$F_n(y)$
of short waves across the channel for the first four waveforms. Considered ‘short’ waves have
$k = 10.45\,\textrm{m}^{-1}$
, corresponding to a wavelength
$\lambda = 0.6$
m, which is five times smaller than the channel width
$2b = 3$
m. The short-wave profiles follow the shape of the ice thickness. The amplitudes shift into the grey area, where the ice is thinnest. In the segments with linear thickness (
$ y \in [-1, -0.7]$
and
$ y \in [0.3, 1]$
), the amplitudes are nearly zero. Similar to long waves, an increase in curvature is observed at
$ y = -0.7$
for the green shape. In the grey area, at a slight distance from the boundaries, the wave profiles are visually nearly identical for all studied variable thickness shapes.
Figure 12 shows the dimensionless maximum strains,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, corresponding to the short-wave profiles in figure 11. For all waveforms in the channel with constant thickness, maximum strains for short waves occur between the channel walls, and these are longitudinal strains. These strains are much larger than for long waves in both constant and variable thickness cases (blue and red shapes).
The maximum strain distribution follows the ice thickness shape. Maximum strains shift into the grey segment, where the ice is thinnest, and decrease to nearly zero in the thicker segments. The strain direction remains the same. For variable thickness, strains at the left-hand boundary (
$y = -0.7$
) increase for the red and green shapes, but not as much as for long waves. For the first waveform, strains at this point are nearly unchanged. Only for the green shape the maximum strain occurs for long waves, not short waves.
Wave profiles
$F_n(y)$
across the channel of the short waves,
$k = 10.45\,\textrm{m}^{-1}$
, for the first four waveforms. See figure 5 for the colour–shape relations. Grey area corresponds to the middle segment with thinnest and constant part of the ice cover.
The dimensionless maximum strains along the phase,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, for the short waves,
$k = 10.45\,\textrm{m}^{-1}$
, for the first four waveforms. See legends for strain direction.
4.1.5. Absolute maximum strain analysis
Figure 13(a) shows the absolute maximum dimensionless strain,
$\varepsilon _{ABS} = \max {(\varepsilon (x,y))}$
, as a function of wavenumber
$k$
for the first four waveforms. The locations of these absolute maximums across the channel are shown in figure 13(b). Note that these locations for constant thickness (black line) are symmetric about the channel’s centre and can be reflected symmetrically.
The absolute maximum of dimensionless strain,
$\varepsilon _{ABS} = \max {(\varepsilon (x,y))}$
, for the first four waveforms as functions of wavenumber
$k$
(a). See figure 5 for the colour–shape relations. Locations of these absolute maximums across the channel as functions of wavenumber
$k$
(b). Grey area corresponds to the middle segment with thinnest and constant part of the ice cover.
For longitudinal strains (
$w_{xx}$
), the maximum strains are nearly independent of the waveform number. For transverse strains, the maximum strains strongly depend on the waveform number. The variable ice thickness affects the wavenumber at which the strain transitions from transverse to longitudinal. This transition can occur sooner (blue shape) or later (red and green shapes) as the wavenumber increases. For example, for the first waveform, the strain becomes longitudinal at
$k \approx 3\,\textrm{m}^{-1}$
for constant thickness, at
$k \approx 2.6\,\textrm{m}^{-1}$
for the blue shape, at
$k \approx 4.7\,\textrm{m}^{-1}$
for the red shape, and at
$k \approx 8\,\textrm{m}^{-1}$
for the green shape.
As shown in figure 13(a), increasing the rate of thickness variation causes a consistent decrease in longitudinal strains for each waveform. Transverse strains for the blue and red shapes are usually lower than those for constant thickness, but for the green shape, the strains are larger for long waves. This effect is most pronounced for the first waveform.
The locations of the maximum strains, shown in figure 13(b), are linked to the type of strain (see figure 13
a). For constant thickness ice, the maximum strains are transverse at the channel walls and longitudinal between the walls. For the red and green shapes, the maximum strain locations behave similarly. These strains are located at the boundary of segments (
$y = -0.7$
) for long waves, as shown in figure 10. As the wavenumber increases, the locations shift inward to the local maximum strain peaks, which are compressed in the middle segment (showed by grey area). The wavenumber at which this change occurs increases with both the waveform number and the rate of variation in ice thickness for the red and green shapes.
The behaviour for the blue shape is more complex. Initially, the maximum strains are located at the right-hand wall for waveforms 1–3. As
$k$
increases, the location shifts sharply to the left-hand boundary at
$k \approx 2.3\,\textrm{m}^{-1}$
for the first waveform, and
$k \approx 2.6\,\textrm{m}^{-1}$
for the second waveform, while remaining transverse. As
$k$
continues to increase, the locations change direction to longitudinal, moving to the local peak in the middle segment for the first waveform, or to one of the local peaks for subsequent waveforms. For the third waveform, the location shifts sharply to a peak in the middle segment at
$k \approx 1.6\,\textrm{m}^{-1}$
, then moves to the right-hand boundary at
$y = 0.3$
at
$k \approx 2.05\,\textrm{m}^{-1}$
, before changing direction to longitudinal at a local peak in the middle segment. For the fourth waveform, the absolute maximum strains are never at the channel walls but are always at one of the local peaks inside the channel or at the segment boundaries. Initially transverse, these strains change to longitudinal at
$k \approx 6.2\,\textrm{m}^{-1}$
.
4.1.6. Critical speeds
The critical speed, defined as the minimum phase velocity, is crucial for the movement of loads on the ice cover. When a load moves, it generates waves propagating at the same phase velocity as the load’s speed. If the load moves at the minimum phase speed, the energy of the wave, propagating at group speed, accumulates beneath the load, causing large ice deflections and strains. The motion of a load at critical or near-critical speeds on an unbounded ice cover has been studied (e.g. Hosking et al. Reference Hosking, Sneyd and Waugh1988; Kozin et al., 2017). As shown in figure 8, phase velocities are highly sensitive to the ice thickness shape. We shall study the dynamics of ice thickness variation, ranging from nearly constant to highly variable shapes. The rate of variation is described by parameter
$r$
, where
$r = 1$
represents minimal thickness variation and
$r = 151$
corresponds to maximum variation, the last one is shown in figure 4.
Figure 14(a) shows the relationship between
$r$
and the ratio of average thicknesses between the left-hand and middle segments,
$ h_{avg}^{(1)} /h_{avg}^{(2)}$
. For the blue shape, this ratio varies from 1.03 to 2.5; for the red shape, from 1.8 to 6; for the green shape, from 2.9 to 20.
The ratios
$ h_{avg}^{(1)} / h_{avg}^{(2)}$
as functions of the parameter
$r$
describing the variation of the ice plate’s thickness shape (a). See figure 4 for the colour–shape relations. The values of critical speeds for the first three waveforms (b). Black dots indicate the values of the corresponding critical speeds for the ice plate with constant thickness.
Figure 14(b) shows the change in critical speeds as a function of
$r$
. All critical speeds decrease monotonically as
$r$
increases. A higher rate of variation in thickness results in a lower critical speed. For example, for
$r = 151$
, the critical speed of the first waveform decreases by a factor of 1.22 for the blue shape, 1.48 for the red shape and 1.86 for the green shape compared with
$r = 1$
. As
$r$
increases, the difference between critical speeds for different thickness shapes grows, meaning that the reduction in critical speed becomes more pronounced. For the green shape, the critical speeds are nearly the same at large
$r$
. In this case, the ice’s response to a moving load will be complex and harder to analyse. For variable ice thickness, the critical speeds of some waveforms can be lower than the first critical speed for constant thickness. This implies that a load moving at subcritical speed on a constant-thickness ice plate may generate multiple hydroelastic waves in areas with variable thickness.
4.2. Ice thickness being large near the centre of the channel and small near the channel walls
The second series of results shows the influence of an ‘opposite’ ice plate’s thickness shape on the wave characteristics, wave profilesand the distribution of maximum strains. Similarly to the previous section, three similar to each other ice thickness shapes are considered, each with different ratios of average thicknesses on the left- and right-hand edge segments with linearly decreasing thickness (figure 15). The ratio of the average ice thickness on the left- and right-hand segments is as follows: for the figure 3(a) shape,
$ h_{avg}^{(1)} /h_{avg}^{(3)} = 1$
; for the figure 3(b) shape,
$ h_{avg}^{(1)} /h_{avg}^{(3)} = 1.15$
; and for the figure 3(c) shape,
$ h_{avg}^{(1)} /h_{avg}^{(3)} = 1.3$
. For these shapes, the ratio between the average thicknesses is not very large due to the rapid decrease in ice thickness at the channel’s walls. The corresponding considered dimensionless ice thicknesses
$ \tilde {h}_i(y)$
are shown in figure 16. In all cases, the minimum thickness is reached at the right-hand wall:
$\tilde {h}_i(1) = 0.47$
for the blue shape;
$ \tilde {h}_i(1) = 0.3$
for the red shape;
$ \tilde {h}_i(1) = 0.14$
for the green shape. Note that studying shapes with further decrease of the ice thickness at the walls,
$ \tilde {h}_i(\pm 1) \to 0$
, is meaningless under clamped conditions. As will be shown further, as the thickness at the wall decreases, the maximum strains will be located there and will increase rapidly with decreasing thickness. When the plate breaks away from the wall at this point, the edges become free of stresses, meaning that the transition from a plate with clamped conditions to a plate with a free edge occurs through a jump discontinuity in strains. The positions of the boundaries of segments are the same to those considered previously,
$ b_1 = -0.7$
,
$ b_2 = 0.3$
, and remain unchanged in the calculations.
The piecewise ice thickness shapes across the channel used in the second series of calculations.
Non-dimensional piecewise ice thickness distribution
$\tilde {h}_i(y)$
for each case from figure 15 for maximum considered variation in the thickness.
4.2.1. Dispersion relations, phase and group velocities
Figure 17 shows the effect of variable ice thickness on the dispersion relations of the first four waveforms. For the first waveform, the behaviour mirrors that from § 4.1, where the dispersion relations consistently decrease with increasing ice thickness variation, with a pronounced effect for shorter waves. For waveforms starting from the second, the behaviour becomes more complex. The dispersion relations for variable ice thickness can be greater than those for constant thickness, as seen for the blue shape (waveforms 2–4) and the red shape (third waveform). In some cases, the dispersion relations stay close to each other over a wide range of wavenumbers.
Dispersion relations
$\omega _n(k)$
of the first four waveforms. See figure 15 for the colour–shape relations.
In § 4.1, the wave characteristics for the studied ice shapes followed a uniform, monotonic trend. However, for the ‘opposite’ thickness shape (figure 15), the wave characteristics become less uniform. Figure 17 shows dispersion relations for the first six waveforms (figure 17
a,c) and corresponding group velocities (figure 17
b,d) for each ice thickness shape. A key difference here is that the dispersion relations for waveforms 2–5 do not behave uniformly as the wavenumber increases. They can come close together, then separate. For example, for the blue shape,
$\omega _2(k)$
approaches
$\omega _3(k)$
, while
$\omega _4(k)$
approaches
$\omega _5(k)$
, before moving apart.
The group velocities, the derivatives of the dispersion relations, clearly show how the rate of change varies. The lines of groups speeds are not monotonic and can ‘swap’ between waveforms. For instance, for the blue shape, the group speeds for waveforms 2 and 3 swap around
$k \approx 8.5\,\textrm{m}^{-1}$
. The same behaviour is seen for the red and green shapes. The number of waveforms affected and the corresponding wavenumbers also depend on the ice thickness shape and can differ.
For stronger variations in ice thickness (red and green shapes), these transitions are sharper, suggesting a non-smooth behaviour of the dispersion relations at these points. It may appear that the frequencies of lower-numbered waveforms surpass those of subsequent waveforms, as seen for waveforms 4 and 5 of the green shape (figure 18(a
iii) at
$k \approx 12.5\,\textrm{m}^{-1}$
). However, no intersections of the dispersion relations were found; they remain very close in these regions. The group velocity change is so sharp that it requires a finer calculation step in
$k$
, as shown in figure 19. This indicates a strong influence of the ice thickness shape combined with the wavenumber. Hydroelastic waves in the ice cover may behave very differently for wavenumbers before and after such points.
Dispersion relations
$\omega _n(k)$
(a) and corresponding group speeds (b) of the first six waveforms. Panels (a i) and (b i) represent results for the blue shape, panels (a ii) and (b ii) for the red shape and panels (a iii) and (b iii) for the green shape, see figure 15.
Zoomed-in figure 18(a
i
ii
,b
iii) for the fourth and fifth waveforms with a small step in
$k$
.
4.2.2. Absolute and maximum strains analysis
Figure 20(a) shows the absolute maximum dimensionless strain,
$\varepsilon _{ABS} = \max {(\varepsilon (x,y))}$
, for the first four waveforms. For the ice plate with constant thickness, the locations of the absolute maximums are symmetrical about the channel’s central line and can be reflected symmetrically.
For variable thickness, all strains are transverse and occur at the channel walls. For the green and red shapes, these maximums may not be exactly at the right-hand wall, where the minimum thickness across the shape is located. Generally, strain values increase with higher wavenumber, higher waveform number and greater variation in thickness. However, for waveforms starting from the second, strains may first decrease and then increase within certain wavenumber ranges. For the green shape, this happens smoothly for the third waveform and sharply for the fourth waveform. These fluctuations correspond to sharp changes in the locations of the absolute maximums from one wall to the other and align with wavenumbers where the dispersion relations are very close, as discussed for figure 18.
Overall, for all studied ice thickness shapes, maximum strains occur at the walls and are greater than those for the constant thickness ice plate.
The absolute maximum of dimensionless strain,
$\max {(\varepsilon (x,y))}$
, for the first four waveforms as functions of wavenumber
$k$
(a). See figure 15 for the colour–shape relations. Locations of these absolute maximums across the channel as functions of wavenumber
$k$
(b). Grey area corresponds to the middle segment with thickest part of the ice cover with constant thickness.
4.2.3. First mode wave profile and strains analysis
The behaviour of the first waveform is the most uniform and monotonic. The only exception is for the blue shape, where the absolute maximum of long-wave strains initially lies at the left-hand wall but sharply shifts to the right-hand wall as the wavenumber increases.
The wave profiles and maximum strain distribution across the channel for the first waveform follow the expected pattern, consistent with the ice thickness shape. Figure 21 shows the wave profiles
$F_1(y)$
and corresponding maximum strains,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, for long waves (
$k = 0.4\,\textrm{m}^{-1}$
) and short waves (
$k = 10.45\,\textrm{m}^{-1}$
). For long waves, the profiles are nearly unaffected by variable thickness, with maximum strains increasing at the walls, especially at the right-hand wall.
For short waves, the profiles shift towards the region of smallest thickness, which is the right-hand segment with linear thickness. Oscillation amplitudes vanish in the middle and left-hand segments. As a result, strains at the left-hand wall approach zero, while those at the right-hand wall increase rapidly. These strains are transverse. The maximum longitudinal strains, initially at the channel centre for the constant-thickness case, shift towards the right-hand wall and decrease in amplitude.
Wave profiles
$F_1(y)$
across the channel for the long waves (a
i)
$k = 0.4\,\mathrm{m^{-1}}$
, and for the short waves (a
ii)
$k = 10.45\,\mathrm{m^{-1}}$
, of the first waveform. Corresponding maximum strains along the phase,
$\max _x(\varepsilon (x,y))$
(b). See figure 15 for the colour–shape relation. Grey area corresponds to the middle segment with thickest part of the ice cover with constant thickness.
4.2.4. Higher-order wavemodes profiles and strains analysis
The behaviour of wave profiles and the distribution of maximum strains becomes more complex for waveforms starting from the second one as the wavenumber
$k$
increases. Figures 22 and 23 show the wave profiles
$F_n(y)$
and the corresponding maximum strains,
$\varepsilon _{max}(y) = \max _x(\varepsilon (x,y))$
, for waveforms 2–4. The strains are normalized with respect to the maximum value on each plot,
$\varepsilon _m$
(for non-normalized strains, see figure 20).
Wave profiles
$F_n(y)$
across the channel for wavenumbers shown in the legends. See figure 15 for the colour–shape relations. Panel (a) is the second waveform, (b) is the third and (c) is the fourth. Grey area corresponds to the middle segment with thickest part of the ice cover with constant thickness.
Normalized maximum strains along the phase,
$\max _x(\varepsilon (x,y)) /{\varepsilon _m}$
, for wavenumbers shown in the legends. See figure 15 for the colour–shape relations. Panel (a) is the second waveform, (b) is the third and (c) is the fourth. For the non-normalized values of maximum strains, see figure 20. Grey area corresponds to the middle segment with thickest part of the ice cover with constant thickness.
For the second waveform, the behaviour is unique only for the blue shape. For the red and green shapes, the wave profiles and strains monotonically shift towards the right-hand segment, concentrating near the wall with the smallest thickness. For the blue shape, wave profiles and strains initially ‘occupies’ the left-hand segment with linear thickness (see figure 22
a for
$k = 6\,\textrm{m}^{-1}$
). At
$k \approx 8.5\,\textrm{m}^{-1}$
, when the dispersion relations
$\omega _2(k)$
and
$\omega _3(k)$
nearly coincide (figure 18
a
i), the wave profile shifts to the right-hand segment, and the maximum strains move from the left-hand wall to the right-hand wall. This transition occurs rapidly over a small range of
$k$
.
A similar behaviour is observed for the third waveform across all shapes. Figures 22(b)and 23(b) show the corresponding results for wavenumbers that cover all positions of the maximum strains in figure 20(b
iii). Generally, at a certain wavenumber, the wave profile shifts towards the left-hand segment. This segment doesnot allow for multiple peaks (e.g. for the third waveform, three peaks are expected), so only one peak is observed, as seen in the wave profiles for the red shape in figure 23(b) for
$k = 10\,\textrm{m}^{-1}$
and the blue shape for
$k = 14\,\textrm{m}^{-1}$
. As
$k$
increases, the right-hand segment with linear thickness becomes more dominant, and wave profiles shift to it. These transitions align with the anomalous behaviour of dispersion relations described earlier.
For the fourth waveform, a similar pattern is seen, with one exception for the blue shape: one of the four oscillation peaks in the wave profile ‘dissolves’ in the left-hand segment, leaving three peaks on the right. The wavenumbers at which these transitions occur depend on the shape and the rate of thickness variation. In some cases, the change in behaviour may be symmetric. For instance, for the red shape, the third waveform’s wave profiles move to the right-hand segment, while the fourth waveform’s wave profiles move to the left.
4.2.5. Critical speeds
We studied the dynamics of the rate of variation in the ice plate’s thickness shape, starting from nearly constant and ending with the extremely variable shape. The rate of variation is described by a natural parameter
$r$
, where
$r = 1$
represents an ice shape with minimal variation and
$r = 151$
with maximum considered variation (see figure 15). Figure 24(a) shows the relationship between the parameter
$r$
and the ratio of average thicknesses between the middle and right-hand segments,
$h_{avg}^{(2)} / h_{avg}^{(3)}$
. The ice thickness for the blue shape varies from nearly constant,
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.01$
, to
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.42$
. For the red shape, it changes from
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.13$
to
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.6$
, and for the green shape, from
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.26$
to
$h_{avg}^{(2)} / h_{avg}^{(3)} = 1.78$
. The change of critical speeds as a function of the parameter
$r$
is shown in figure 24(b). All considered speeds decrease monotonically, except for the first one. The first critical speed initially increases, then decreases. For the blue shape it has a practically unchanged final value for
$r = 151$
. The largest increase in the first critical speed is observed for the blue shape, where the ratio of maximum critical speeds is
$U^{blue}_{1,crit}/U^{constant}_{1,crit} \approx 1.017$
, while for the other two shapes, it is
$U^{red}_{1,crit}/U^{constant}_{1,crit} \approx 1.012$
and
$U^{green}_{1,crit}/U^{constant}_{1,crit} \approx 1.006$
(here, superscripts correspond to the type of shape). Critical speeds for the second waveform decrease very slowly, regardless of the considered ice thickness shape. The behaviour of critical speeds for the third waveform is similar to the first waveform, but only with a monotonically decreasing trend. Note that compared with the critical speeds for the ice thickness from § 4.1, the current critical speeds do not intersect and, overall, change much more slowly and to a lesser extent.
The ratios
$h_{avg}^{(2)} / h_{avg}^{(3)}$
as functions of the parameter
$r$
, which describes the rate of variation in the ice plate’s thickness shape (a). The values of critical speeds for the first three waveforms (b). See figure 15 for the colour–shape relations. Black markers indicate the corresponding critical speeds for a plate with constant thickness.
4.3. Usability of classic normal vibration modes for constant thickness in the considered problem with variable thickness
One might consider that the equation for the wave profile across the channel
\begin{align} \beta \bigg[ \tilde {h}^3_i(y) ( \kappa ^4 F - 2\kappa ^2 F^{\prime \prime } + F^{\prime \prime \prime \prime } ) &+ 6\tilde {h}^2_i(y)\tilde {h}_i^{\prime } ( -\kappa ^2 F^{\prime } + F^{\prime \prime \prime } ) \nonumber \\ &+ 6\tilde {h}_i(y)(\tilde {h}_i^{\prime })^2( F^{\prime \prime } - \mu \kappa ^2F) - \frac {\delta \gamma }{\beta } \tilde {h}_i(y) F \bigg]\nonumber \\ &\quad= \gamma \Phi (y,0) - F, \end{align}
could be solved by decomposing the wave profile
$F$
into classic normal vibration modes of an elastic beam with constant thickness
\begin{equation} F(y)=\sum _{n=1}^{\infty }\widehat {a}_n {\widehat {\psi }}_n(y), \end{equation}
where the normal modes
${\widehat {\psi }}_n(y)$
are solutions of the spectral problem
\begin{equation} {\widehat {\psi }}_{n, yyyy} = \lambda ^4_n {\widehat {\psi }}, \qquad {\widehat {\psi }}_n = {\widehat {\psi }}_{n, y} = 0 \quad (y = \pm 1). \end{equation}
These modes are orthogonal
\begin{equation} \int \limits _{-1}^1 {\widehat {\psi }}_n {\widehat {\psi }}_m \textrm{d}y = \delta _{nm}. \end{equation}
The solutions of (4.2)–(4.4) are even and odd modes not interacting with each other. Even modes have the form
\begin{align} {\widehat {\psi }}^c_n (y) = A_n^c \left (\cos ( \lambda ^c_n y) -B_n^c \cos h (\lambda ^c_n y) \right), \quad B^c_n = \frac {\cos \lambda ^c_n}{\cos h \lambda ^c_n}, \quad A^c_n = \frac 1{\sqrt {1+(B^c_n)^2}}, \end{align}
where
$\lambda _j^c$
are the solutions of the equation
\begin{align} \tan \lambda _n^c= - \tanh \lambda _n^c, \qquad \lambda _n^c=\pi n-\pi /4+\varDelta ^c_n, \end{align}
where
$\varDelta ^c_n\to 0$
as
$n\to \infty$
. Odd modes have the form
\begin{align} {\widehat {\psi }}_n^s (y) = A_n^s \left ( \sin (\lambda _n^s y) - B_n^s \sinh (\lambda _n^s y) \right), \qquad B_n^s = \frac {\sin \lambda ^s_n}{\sinh \lambda ^s_n}, \quad A^s_n = \frac 1{\sqrt {1 - (B^s_n)^2}}, \end{align}
where
$\lambda _n^s$
satisfy the equation
\begin{align} \tan \lambda _n^s=\tanh \lambda _n^s \end{align}
and
$\lambda _n^s = \pi n + \pi /4+\varDelta ^s_n$
, where
$\varDelta ^s_n\to 0$
with
$n\to \infty$
.
Taking into account the separation into even and odd modes, the wave profile
$F$
is sought in the form
\begin{equation} F(y) = \sum _{n=1}^{\infty }\widehat {a}_n {\widehat {\psi }}_n(y), \qquad {\widehat {\psi }}_{2 n - 1}(y) = {\widehat {\psi }}^c_n(y), \quad {\widehat {\psi }}_{2 n}(y) = {\widehat {\psi }}^s_n(y). \end{equation}
The profile of the potential across the channel
$\widehat {\Phi }(y, z)$
is sought using the same method as used in the main problem
\begin{equation} \widehat {\Phi }(y,z) = \sum _{n = 1}^{\infty }\widehat {a}_n\widehat {\Phi }_n(y,z) \end{equation}
and
$\widehat {\Phi }_n$
are the solutions of a boundary value problems similar to (3.15).
Substituting (4.9) and (4.10) into the plate equation, multiplying by
${\widehat {\psi }}_m$
, and integrating the result over
$y$
from
$-1$
to 1, we obtain an infinite algebraic system of equations
\begin{align} &\sum \limits _{n = 1}^{\infty }\widehat {a}_n \left (\beta \left [ \kappa ^4 K^{(3)}_{nm} + K^{(34)}_{nm} + 6 \kappa ^2 S^{(3)}_{nm} + 6 K^{(23)}_{nm} + 6 K^{(12)}_{nm} - 6 \mu \kappa ^2 K^{(1)}_{nm}\right]\right. \nonumber \\ &\qquad \left. - \delta \gamma K^{(0)}_{nm} + \delta _{nm} - M^{P}_{nm}\right) = 0, \end{align}
where
\begin{align} & K^{(3)}_{nm} = \int \limits ^1_{-1} \tilde {h}_i^3 \; {\widehat {\psi }}_n \; {\widehat {\psi }}_m \textrm{d}y, \quad K^{(34)}_{nm} = \int \limits ^1_{-1} \tilde {h}_i^3 \; {\widehat {\psi }}^{\prime \prime \prime \prime }_n \; {\widehat {\psi }}_m \textrm{d}y, \quad S^{(3)}_{nm} = \int \limits ^1_{-1} \tilde {h}_i^3 \; {\widehat {\psi }}^{\prime }_n \;{\widehat {\psi }}^{\prime }_m \textrm{d}y, \end{align}
\begin{align} & K^{(23)}_{nm} = \!\int \limits ^1_{-1}\! \tilde {h}_i^2 \; \tilde {h}^{\prime }_i \; {\widehat {\psi }}^{\prime \prime \prime }_n \; {\widehat {\psi }}_m \textrm{d}y, \!\quad K^{(12)}_{nm} = \!\int \limits ^1_{-1}\! \tilde {h}_i \; (\tilde {h}^{\prime }_i)^2 \; {\widehat {\psi }}^{\prime \prime }_n \; {\widehat {\psi }}_m \textrm{d}y, \!\quad K^{(1)}_{nm} = \!\int \limits ^1_{-1}\! \tilde {h}_i \; (\tilde {h}^{\prime }_i)^2 \; {\widehat {\psi }}_n \; {\widehat {\psi }}_m \textrm{d}y, \end{align}
\begin{align} K^{(0)}_{nm} = \int \limits ^1_{-1} \tilde {h}_i \; {\widehat {\psi }}_n \; {\widehat {\psi }}_m \textrm{d}y, \quad M^{P}_{nm} = \int \limits ^1_{-1} \widehat {\Phi }_n(y, 0) \; {\widehat {\psi }}_m \textrm{d}y. \end{align}
Note that in the general case, when
$\tilde {h}_i(y)$
is an arbitrary piecewise linear function, the integrals in (4.11) will not possess orthogonality and symmetry with respect to the indices
$n$
and
$m$
. One needs to be careful when computing such systems and thorough convergence checks are necessary.
The comparison of dispersion relations
$\omega _n(k)$
calculated using (3.25) and relations
$\widehat {\omega }_n(k)$
calculated as eigenvalues of the system (4.11) with finite values of number of modes is conducted for the previously considered shapes of ice thickness (see figures 4 and 15). The calculations of
${\omega }_n(k)$
were performed for a number of modes from
$N = 10$
to
$N = 30$
depending on the rate of thickness variation, and the calculations of
$\widehat {\omega }_n(k)$
were performed for
$\widehat {N} = 20$
and
$\widehat {N} = 60$
.
First, we studied a small rate of ice thickness variation: the blue shape of the ice plate’s thickness with
$h^{(1)}_{avg} = h^{(3)}_{avg} \approx 0.0036$
m,
$h^{(2)}_{avg} \approx 0.0034$
m, see figure 4(a). The ratios
$ [\max (|\widehat {\omega }_n(k) - \omega _n(k)|/\widehat {\omega }_n(k)) \times 100]$
are shown in figure 25 as functions of the waveform number
$n$
. The calculations of
$\widehat {\omega }_n(k)$
for this figure were performed for
$\widehat {N} = 20$
. The difference in the dispersion relations calculated by both methods is not uniform; furthermore, the largest difference occurs for the fourth waveform, but it does not exceed 0.6 %. These results confirm that the solution for a plate with an extremely small rate of thickness variation can be approximated by classic normal vibration modes with a not large quantity of their number.
The ratios
$(\max (|\widehat {\omega }_n(k) - \omega _n(k)|/\widehat {\omega }_n(k)) \times 100)$
as functions of number of the waveform, n. The blue shape of ice plate’s thickness (see figure 4) with the low variation,
$h^{(1)}_{avg}=h^{(3)}_{avg}\approx 0.0036$
m, is considered.
We now compare frequencies calculated using different methods for large variations in ice thickness. The ratios
$ [\max (|\widehat {\omega }_n(k) - \omega _n(k)|/\widehat {\omega }_n(k)) \times 100]$
for the shapes in figure 5 are shown in figure 26 as functions of the waveform number,
$n$
. The legends provide values of
$\widehat {N}$
. Average thickness values on each segment can be found at the end of § 4.1 and in figure 5. As in the previous figure, the difference in dispersion relations calculated by both methods varies non-uniformly and strongly depends on
$\widehat {N}$
. For
$\widehat {N} = 20$
, the difference reaches 40 % for the 10th waveform of the green shape, while for the blue and red shapes, the differences are around 6%–7 %. Increasing the number of classic vibration modes to 60 showed the following. (i) A larger variation in ice thickness does not always lead to the smallest difference. In our case, the smallest difference was found for the sixth waveform of the blue shape and the tenth waveform of the red shape. (ii) The dispersion relations of the first few waveforms can be calculated with an accuracy of up to 1 %, but the frequencies of higher waveforms are less accurate. We did not find a reasonable number of classic vibration modes (between 60 and 150) that would reduce the difference between the methods below 4 % for all waveforms considered. Similar ratios for the ice thickness shapes in figure 15 are shown in figure 27. Values of average thicknesses can be found in figure 16. In the same way,
$\widehat {N} = 20$
is insufficient for accurate calculations using classic vibration modes. However, for these shapes,
$\widehat {N} = 60$
is sufficient to calculate the dispersion relations of the first 10 waveforms with a difference of up to 1 %. Further increases in the number of vibration modes reduced the calculations difference.
The ratios
$(\max (|\widehat {\omega }_n(k) - \omega _n(k)|/\widehat {\omega }_n(k)) \times 100)$
as functions of number of the waveform, n. See figure 5 for the colour–shape relations.
The ratios
$( \max (|\widehat {\omega }_n(k) - \omega _n(k)|/\widehat {\omega }_n(k)) \times 100)$
as functions of number of the waveform, n. See figure 16 for the colour–shape relations.
In general, it can be concluded that solving the considered problem using classic normal vibration modes of an elastic beam with constant thickness leads to accurate results for ice covers with thickness close to constant or for certain specific shapes of ice cover. However, it is necessary to compute a large number of these modes to achieve required accuracy. The example was provided where a very large number of modes is required, but even then, exact matching of dispersion relations is not guaranteed. Note that while the accuracy of the calculations of dispersion relations may be good, the wave profiles can be more different. We recommend using the constructed vibration modes
$\psi _n(y)$
of a beam with variable thickness for modelling ice covers of variable thickness with complex shapes.
5. Conclusions
The problem of periodic hydroelastic waves in a frozen channel with variable ice thickness was studied within the linear theory of hydroelasticity. The normal modes of the ice plate were computed using a piecewise-linear approximation, ensuring continuity of deflections, slopes, stresses/strains, bending moments and shear forces. In regions with linear thickness variation, modes were expressed via Bessel functions, while in regions with constant thickness, trigonometric and hyperbolic functions were used. Used linearized boundary conditions at the ice–liquid interface in the considered method are correct when, in addition to the classical assumptions of linear theory, the derivative of the ice thickness profile
$h_{i,y}(y)$
is of the same order of magnitude as the first spatial derivatives of the ice deflections,
$w_x$
or
$w_y$
, or smaller.
Two cases of ice thickness distribution across the channel were analysed. In Case I, the thickness was constant in the middle and increased linearly at the edges. In Case II, the middle thickness was larger and constant, while the edges had a decreasing linear profile. Each case included three subcases with different thickness variation rates. The average thickness remained constant.
A method of computation of the eigenvalues of the normal modes was constructed taking into account large values of Bessel functions and their derivatives. Using scaled modes, up to 40–50 modes (depending on the case) were accurately calculated, ensuring convergence to the constant-thickness solution as variation decreased.
For both cases, periodic hydroelastic waves propagating with constant speed and amplitude were found, ordered by increasing frequency. The first wave exhibited predictable behaviour, while higher-order waves showed complexity depending on the thickness profile. Contributions from different normal modes to the wave profiles depend significantly on the shape of the ice thickness and vary from the case with ice of constant thickness.
Hydroelastic waves with wavelengths
$\lambda \gt 0.42$
m and wavenumbers
$0 \lt k \lt 15$
m
$^{-1}$
were analysed. Long-wave characteristics were nearly independent of thickness variation, whereas short-wave properties were strongly affected. In Case I, increasing the thickness variation rate reduced all wave characteristics. In Case II, a similar reduction was observed, but certain wavenumber intervals showed an increase. Dispersion relations behaved monotonically in Case I but exhibited complex non-monotonic behaviour in Case II, with frequencies of adjacent waves converging at specific wavenumbers. However, no wave overlap was observed.
The interplay between ice thickness and wavenumber was a key factor. For long waves, thickness variation had little effect on wave profiles but significantly influenced maximum strain locations and values. Short waves were highly sensitive to thickness shape. In Case I, oscillations and maximum strains shifted to thinner regions, with strain values decreasing as variation increased. For high variation rates, strain peaks appeared at segment boundaries, with transverse strains exceeding longitudinal ones. Long-wave strains increased at segment boundaries where the thickness changed from linear to constant.
In Case II, long-wave profiles remained unchanged, but strain increased at channel walls near thin ice regions. Short-wave profiles varied with wavenumber, concentrating in thin regions. Oscillations grouped on either the left- or right-hand segment, sometimes switching sharply from one to another at certain wavenumbers. This transition can occur simultaneously in opposite directions for adjacent frequency waves.
Critical velocities decreased in all cases except for the first wave in Case II, where a slight increase was observed.
Using classic normal modes for a beam with constant thickness showed comparable to our method results for nearly constant thickness or specific profiles but required large number of modes. The normal modes constructed in this paper are recommended for modelling ice covers with complex thickness distributions. The results demonstrated that periodic hydroelastic wave characteristics strongly depend on the ice thickness distribution across the channel.
Understanding how thickness variation affects all wave characteristics enables a practical application in solving the inverse problem. Specifically, by analysing the non-uniformity of hydroelastic waves propagating through the ice field, it becomes possible to determine variations in the ice cover thickness.
Results will be useful for laboratory experiments with ice frozen in natural environments without additional installations, conducted in a similar way as it was done by Pogorelova et. al. (Reference Pogorelova, Zemlyak and Kozin2019a ) and with same-order parameters of the problem.
Funding
This work is supported by the National Natural Science Foundation of China (nos. 52350410468, W2433140, 52192693, 52192690, 52371270, U20A20327), to which the authors are most grateful.
Declaration of interests
The authors report no conflict of interest.
