Wave Resistance Caused by a Point Load Steadily Moving on the Surface of a Floating Viscoelastic Plate

: The wave resistance caused by a point load steadily moving on an inﬁnitely extended viscoelastic plate ﬂoating on an inviscid ﬂuid is analytically studied, which can be used to describe the response due to the motion of amphibious air-cushion vehicles on the continuous ice sheet on the ocean. The action of concentrated and distributed point loads are both considered. Under the assumptions that the ﬂuid is incompressible and homogeneous and the motion of the ﬂuid is irrotational, the Laplace equation is taken as the governing equation. For the ﬂoating plate, the Kelvin–Voigt viscoelastic model is employed. At the plate–ﬂuid interface, linearized boundary conditions are used when the wave amplitude generated is less than its wavelength. The Fourier integral transform is performed to achieve the formal solution. The residue theorem is applied to derive the response of ﬂexural–gravity wave resistance. It is indicated that for a point load with a uniform rectilinear motion, the wave resistance shows a sharp decrease with the increase in the moving speed when the load velocity is greater than the minimum phase velocity. There is no steady wave resistance when the load velocity is smaller than the minimum phase velocity. The effects of different parameters are obtained. Wave resistance decreases with the increasing plate thickness, viscoelastic parameter, and Poisson’s ratio, especially for a small value of viscoelastic parameter.


Introduction
To make use of polar resources, it is necessary to strengthen the monitoring and research of polar hydrology, meteorology, sea ice, etc. Due to the complex and changeable environment in polar regions, ships passing through polar regions must be upgraded to be suitable for ice navigation. In nature, the interaction between waves and sea ice is one of the main dynamic reasons for the evolution of polar waves and for the generation or movement of sea ice (Ni et al., 2020) [1].
In studying the hydroelasticity of sea ice, the wave-plate and the wave-beam models are generally used (Amouzadrad et al., 2022) [2], combining with the potential flow theory to solve an initial-boundary value problem related to the Laplace equation for an incompressible fluid. Considering the hydroelastic response of ices, most researchers have assumed that the ice sheet is an infinitely long elastic plate (Kheysin 1963, Nevel 1970, Davys et al. 1985, Schulkes andSneyd 1988, Pogorelova 2011, Tkacheva 2019a, Tkacheva 2019b, Stepanyants and Sturova 2021) [3][4][5][6][7][8][9][10]. The takeoff and landing experiment of the LC-130 Hercules transport aircraft was carried out by Squire (1988) [11] on the ice surface of McMurdo, Antarctica, which verified the correctness of the uniform elastic plate theory. Calculations by Strathdee et al. (1991) [12] showed that for distances 20 times greater than the ice thickness, the accuracy of sheet theory is within 5%. Amouzadrad et al. 2022 [2] presented, on basis of the Timoshenko-Mindlin beam theory, a simplified two-dimensional mathematical model of wave interaction with a moored flexible floating rectangular beam in finite water depth.
Ice sheets can also be regarded as more complex viscoelastic models. There are three commonly used viscoelastic models: Maxwell's model (a spring and a damper in series), Kelvin's model (a spring and a damper in parallel), and Maxwell-Kelvin's model (the concatenation of the above two models). Some scholars have applied the viscoelastic model to floating ice. Hutter (1975) [13] introduced a general viscoelastic plate theory for floating ice. Squire and Allan (1980) [14] mentioned that the incident ocean wave action is one of the main reasons for ice breaking, and they used Maxwell-Voigt's model for their study. Hosking et al. (1988) [15] considered a two-parameter memory function which was employed to describe the viscoelastic properties of ice, which were assessed via Maxwell's model. Wang et al. (2004) [16] investigated the transient response caused by a sudden starting point load on a floating viscoelastic plate, for which the Boltzmann hereditary delay integral approach was used to predict the dynamic response. Kozin and Pogorelova (2009) [17] compared the results of the three kinds of viscoelastic models with the experimental results. Shishmarev et al. (2023) [18] considered the viscoelastic plate in terms of Kelvin-Voigt's model.
Beyond linear models, nonlinear ones have been considered. Agarwala and Nair (2004) [19] analyzed the floating airport with the aid of a one-dimensional Timoshenko-Mindlin plate equation, which was solved via Fourier transform in the wave number domain. Dinvay et al. (2019) [20] applied a fully dispersive weak nonlinear model to describe the motion of an ice sheet with moving loads, and the model is general enough to study the motion of arbitrary shape loads. Dinvay et al. (2022) [21] investigated the waves generated by the moving loads on ice plates floating on an incompressible fluid. Two different viscoelastic approximations were considered for the ice cover: a model depending on the strain relaxation time, and another including a hereditary delay integral. The problem was formulated in terms of the exact dispersion relation and the Dirichlet-Neumann operator connected to the fluid motion. Pogorelova et al. (2023) [22] touched upon an unsteady rectilinear motion of a submarine in a liquid under an ice cover, using the Kelvin-Voigt model for a viscoelastic plate.
Moving loads forced on an ice sheet induce steady-state wave modes and pressures. Earlier studies assumed that moving loads move in a straight line with a uniform velocity on a linear elastic sheet. Kheysin (1963) [3] studied the deflection of the plate forced by a line load and a point load moving uniformly on an elastic plate. Nevel (1970) [4] revised the conclusion about point loads proposed by Kheysin (1963) [3], considering the displacement at an infinite time and ignoring the transient movement, and the conclusion can be extended to the case of circular distributed load. Davys et al. (1985) [5] studied the steady-state wave mode of an elastic plate under a moving point load with a uniform speed, which does not exist at the critical velocity. Pogorelova (2011) [7] studied the threedimensional unsteady problem of a floating elastic plate under a uniformly moving point load in an infinitely deep fluid. There are also scholars who have investigated the response to a sudden start of a moving source. Schulkes and Sneyd (1988) [6] studied the elastic plate under suddenly starting line load, using the stationary phase method for a large time domain, and the results show that the deflection of the plate grows by t 1/3 when the moving velocity is equal to the maximum of group velocity, where t is the time. Nugroho et al. (1999) [23] considered the elastic plate under suddenly starting point loads. Shishmarev et al. (2023) [18] showed that the variation of the ice thickness significantly changes the characteristics of the hydroelastic waves in a channel.  [20], exploring the response of floating ice sheets to loads moving along curved paths.
Some scholars considered the shape or thickness of the plate. Pogorelova et al. (2016) [26] considered the hydroelastic response due to a moving load acting on a thin plate with variable thickness. Meylan (2021) [27] analyzed the time-dependent motion of a circular elastic plate floating on the surface.
The response of the floating ice sheet to the moving load will lead to the generation of structural waves with added drag and undulations, which may permit the safe operation of aircraft or amphibious air-cushion vehicles (AACVs). Yeung and Kim (2000) [28] described the wave resistance as a kind of force of drag. In order to ensure the safety of the transportation, whatever the aircraft or AACV, eliminating the wave resistance has attracted attention. The formulation of wave resistance can be obtained by integrating the deflection expression of the plate or directly using semi-empirical formulae after experiments. Doctors et al. (1972) [29] dealt with the theoretical wave resistance of an air-cushion vehicle traveling over water of uniform finite or infinite depth in a steady or unsteady motion. Liu and Papanikolaou (2016) [30] combined the speed of the ship and the wave flow with the main characteristic parameters to quickly estimate the wave resistance of the ship in the head wave. Lang and Mao (2020) [31] predicted the loss of ship speed in sea areas with the aid of a semi-empirical model. Yu et al. (2021) [32] proposed a semi-empirical formula that can directly calculate wave resistance.
Nowadays, researchers have made some progress in the field of wave resistance led by the load on an ice sheet in deep water. Kozin and Milovanova (1996) [33] considered the stationary problem of wave resistance to AACV in broken ice. Kozin and Pogprelova (2003) [34] discussed the wave resistance caused by the rectangular load moving at variable speeds on the ice sheet and obtained the numerical results. Lu and Zhang (2013) [35] obtained the analytical solution to hydroelastic wave resistance due to the steady movement of a linear source over an infinitely deep fluid surface within the framework of linear potential theory. Li and Lu (2017) [36] simplified the elastic sheet into a viscoelastic version. The asymptotic solution of the wave resistance generated by the line point source was given. Pogorelova and Kozin (2019) [37] developed the influence of the current on wave resistance generated by an air-cushion vehicle that starts suddenly or moves at a constant speed and proved that the transverse current basically has no effect on the wave resistance, whereas the longitudinal wave flow has a significant effect. Some scholars use the stationary phase method or the steepest descent method to expand the integral equation asymptotically. Nugroho et al. (1999) [23] used the asymptotic solution of Cooke (1982) [38]. Schulkes and Sneyd (1988) [6] gave the asymptotic expansion of integral in the large time domain with the aid of the steepest descent method, transforming the integral curve. The asymptotic expansion of the double integral was obtained. The form of the final results obtained in Nugroho et al. (1999) [23] is the same as that in Schulkes and Sneyd (1988) [6].
Yeung and Kim (2000) [28] considered the response which included the flection and wave resistance of an aircraft moving steadily on a floating mat-flat airport and investigated the analytical results of a point load on an elastic plate. Li and Lu (2017) [36] gave the analytical solution of wave resistance generated by a line load. In order to identify some differences between two-and three-dimensional models and the necessity of considering the parameter of viscoelasticity, the point load on a viscoelastic plate in terms of Kelvin-Voigt's model is employed to study the wave resistance. The mathematical model is established with linear potential theory under some basic assumptions. The formal solution is derived with the aid of the Fourier integral transform method. Then, the residue theorem and numerical calculation are employed to simplify the solutions as asymptotic results. A flowchart that briefly shows the process of the methodology in the present work is given in Figure 1.

Mathematical Formulation
The hydroelastic problem of AACV moving on an infinite floating ice sheet is considered. The mathematical model diagram is shown in Figure 2. AACV is modeled as an external point load. The ice sheet is deemed to be an infinitely long viscoelastic plate floating on water of finite depth. The isotropic homogeneous hypothesis is used for the plate. The initially stationary fluid has a uniform depth H. It is assumed that the fluid is inviscid, incompressible, and homogeneous, and the motion of the fluid is irrotational. The Cartesian coordinate system oxyz is taken. The oxy plane coincides with the stationary horizontal plane, and oz is the vertical axis. The plate has a thickness of d and a draft of zero. Let ρ, g, and t represent the density of the fluid, the gravitational acceleration, and the time, respectively.  Under the irrotational motion hypothesis on the fluid, the Laplace equation with the velocity potential φ is established to describe the flow, namely, where ∇ 2 is the Laplace operator and ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 + ∂ 2 /∂z 2 .
For a fluid of finite depth, the velocity potential satisfies the non-penetrating boundary condition at the bottom: Let ζ(x, y, t) indicate the vertical deflection of the plate under an external load P ext (x, y, t) moving on a plate. Assuming that the amplitude of the generated wave is smaller than the wavelength, we use linearized boundary conditions at the junction of the plate and the fluid. The kinematic boundary condition at the interface (z = 0) between the viscoelastic thin plate and the fluid is as follows: Equation (3) shows that in the process of fluctuation, the fluid always adheres to the viscoelastic thin plate without separation. By using the linear Kelvin-Voigt's viscoelastic model, which was mentioned in Kozin and Pogorelova (2009) [17], very large floating structures (VLFS) are simplified as a viscoelastic thin plate. Then the dynamic boundary condition at z = 0 is as follows: where with E, ν, ρ e , and µ representing Young's modulus, Poisson's ratio, the density, and the viscosity parameter the of viscoelastic plate, respectively. D is the bending stiffness coefficient of the viscoelastic thin plate, τ the viscoelastic relaxation time of Kelvin-Voigt's model, and M the mass of the plate. In Equation (4), ∇ 4 = ∇ 2 ∇ 2 = ∂ 4 /∂x 4 + 2∂ 4 /∂x 2 ∂y 2 + ∂ 4 /∂y 4 for the plane of the plate. From Equation (4), it can be seen that the balance among the viscoelastic force and the inertial force of the thin plate, the hydrodynamic force of the fluid, and the external load leads to the deflection equation of the plate. In order to obtain the solution in integral form, we introduce the Fourier integral transform as follows: Applying the transform (6) on Equations (1) and (2), we obtain the following: where k 2 = α 2 + β 2 . By the transform (6), Equations (3) and (4) After solving the above differential equations, the vibration equation of the plate is obtained; for z = 0, Equation (13) represents the dispersion relationship of the flexural-gravity wave caused by the moving source on the surface of viscoelastic plate. According to the wave resistance formula in Cohen et al. (2001) [39], the relationship between the wave resistance R, the wave amplitude ζ, and external loadP ext can be obtained: whereP * (k, t) is the conjugate function ofP ext .

Analytical Results for the Contributed Uniform Load
Consider a point load steadily moveing along the positive x-axis at a uniform speed V on z = 0; the loaded pressure function is expressed as where P 0 is a constant strength and δ(·) Dirac's function. The Fourier transform of P ext (x, y, t) reads as follows:P Substitution of Equation (17) into Equation (11) yields the deflection equation of the viscoelastic thin plate subjected to the steadily moving load, for which the steady-state solution, neglecting the transient effects, reads as follows: Since the wavelength of viscoelastic thin plate is greater than the thickness of the thin plate, the inertial force term of the thin plate in the dynamic boundary condition in the plate-fluid interface could be ignored. Thus σ → 0 is taken hereinafter.
The polar coordinate α = k cos θ, β = k sin θ is introduced. By substituting Equation (18) into Equation (15), the steady-state wave resistance R is represented as where Obviously, c(k) is the phase speed of the hydroelasitc waves, of which the minimum is denoted by c min .
With the change of ς = e iθ , Equation (20) is deduced as where b = A/B. Obviously, there are three zeros (denoted by ς 0 , ς 1 , and ς 2 ) of the denominator of the integrand in Equation (22): Appendix A demonstrates that |ς 1 | < 1 and |ς 2 | > 1 hold for V > c min , while |ς 1 | > 1 and |ς 2 | < 1 hold for V < c min . The value of the integral in Equation (22) depends on the residues contributed from ς 0 and ς 1 when V > c min and from ς 0 and ς 2 when V < c min .
With the aid of the Cauchy residue theorem, the integral in Equation (22) is performed. Upon some straightforward manipulation, Equation (19) becomes where with L 0 = V 2 /g, and sgn(·) denoting the sign function. The residue theorem is still used to derive the integration of k. Then, the final expression of the wave resistance can be given as follows: Res k=k j F(k), k j (28) where k j is the root of P 2 + Q 2 = 0 and n is the number of total roots.

Wave Resistance for the Distributed Circular Load
Now we consider the distributed force of action. The loading function P ext is represented as where P 0 represents the magnitude of loading,R an effective radius of the loading, and f (r) = e −πr 2 the shape of the loading as a Gaussian distribution, normalized as 2π ∞ 0 f (r)rdr = 1 (Yeung and Kim 2000) [28]. The inversed point loadP which acts as an Gaussian distribution isP where r represents the distance from the center of the load to the origin, and J 0 (·) the zerothorder Bessel function of the first kind. A substitution of Equation (30) into Equation (15) yields the equation of the wave resistance: where I(k) is given in Equation (20). By the residue theorem, we give the solution directly: where where P(k) and Q(k) are given in Equations (26) and (27), respectively.

Results
With the physical parameters in Squire et al. (1996) [40], i.e., E = 5 Gpa, H = 350 m, ρ = 1024 kg/m 3 , g = 9.8 m/s 2 , and ρ e = 917 kg/m 3 , numerical analysis is applied to obtain the solutions of the equation P 2 + Q 2 = 0. Changes of non-dimensional wave resistance with different parameters are presented in the following figures. In general, for the case of a moving source with a fixed speed, when the viscoelastic parameter τ is smaller, the other parameters will have more drastic effects. Figure 3 shows that the wave resistance declines monotonically as the moving speed grows with V > c min . There is no steady wave resistance when V < c min . When V = c min , with the increase in τ, the wave resistance will reach its peak value. When the value of τ drops, the wave resistance will be considerably decreased. For a larger value of τ, the attenuation rate of the wave resistance is bigger than in other situations. When the velocity is growing, the wave resistance tends to more likely remain a constant.
As shown in Figure 4, adding more plate thickness causes the wave resistance to decrease. The maximum wave resistance will be produced by small values of the parameter when d = 1 m, which takes the viscoelasticity of the plate at a small value into account. With larger value of τ, the wave resistance is more likely to be very small. However, when τ = 0.001 s, the line of wave resistance converges quickly with a rise of d. For other conditions, the curves always have a relatively steady decrease. With the addition of d, the three curves all tend towards a minuscule value that is close to 0.   Figure 5 demonstrates that as τ is small, the wave resistance decreases dramatically with the increase in τ. The smaller the velocity of the moving load, the larger the wave resistance. With a bigger value of V, the larger the slope of the curve. However, when V is smaller, the wave resistance R descends relatively slowly.
As indicated in Figure 6, there are significant differences among the curves of the wave resistance. When Poisson's ratio ν tends to be a tiny value, the wave resistance will go up with the decrease in τ. For example, when τ = 0.001 s, the wave resistance is much larger than those in other conditions. Moreover, there is an obvious decrease with the increase in ν. When τ = 0.01 s, the curve tends to be gentle and approaches a value finally. When τ = 0.1 s, the wave resistance will be a much smaller value, i.e., near 0, and there is almost no fluctuation.
For the distributed circle load, Figure 7 demonstrates that when a larger radiusR of distributed load is inflicted, the added wave resistance will be much greater. Figure 8 indicates that the small parameter of viscoelastic τ will also play an important role in this kind of situation. If τ is small enough, the wave resistance will be great. In comparison with the numerical results of Yeung and Kim (2000) [28], we find a similar trend for the wave resistance.    In conclusion, if the velocity is equal to the minimum of the group velocity, the wave resistance is not effected by the viscoelastic parameter τ. When the velocity of the moving source is bigger than c min , it has a significant difference from the two-dimensional model investigated by Li and Lu (2017) [36]. In the three-dimensional model considered here, the wave resistance monotonically changes as the viscoelastic parameter varies. The greater value of the parameter which represents the viscoelasticity, the smaller the wave resistance.

Conclusions
The steady-state wave resistance on a viscoelastic plate at the surface of a fluid caused by the motion of a point load in a monolayer fluid has been investigated. The analytical solution and numerical results demonstrate the impact of physical quantities on the wave resistance, namely, the thickness of the plate, the Poisson ratio, and the viscoelastic parameter.
To sum up, it is found that for the case of a uniform moving source, there is no steady wave resistance when V < c min . It is clear that the wave resistance is monotonically diminishing when the velocity of the moving source exceeds the minimum of the phase velocity. When the viscoelastic parameter is lower, the wave resistance decreases fast. After a specific value, it will decline steadily. It is found that an increasing Poisson's ratio (and an increasing thickness of the plate, too) will lead to a decreasing wave resistance. However, when a greater viscoelastic parameter exits, the effect of Poisson's ratio on the wave resistance could be nearly ignored. Based on the analytical solution presented here, we can recognize the influences of various physical parameters on hydroelastic wave resistance. Apart from the change in physical quantities, the form of the loads also undergoes a change. The wave resistance, which varies directly with the distribution of the loads, decreases monotonically as the radius of the load varies. An increase in the radius of the distributed load leads to a corresponding increase in wave resistance as well.
Perfect fluid and small-amplitude wave hypothesis are made for analytical solutions in the present work. The sea ice model is deemed to be a homogeneous, continuous, infinitely long thin plate, where actually, the sea ice model is much more complicated. The ice sheet could also be inhomogeneous and polyporous. The effects of lateral stress on the ice sheet should also be considered. For analytical solutions of complex sea ice models, further study is expected.
Author Contributions: Investigation, methodology, calculation, visualization, validation, formal analysis and writing-original draft preparation, Z.W.; Conceptualization, methodology, resources, supervision, project administration, funding acquisition and writing-review and editing, D.L. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.