# On Unsteady Effects in WIG Craft Aerodynamics

Authors: Nikolai Kornev
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

The paper presents the analysis of unsteady forces and their influence on the aerodynamics and motion of a wing-in-ground (WIG) effect craft. Two-dimensional and three-dimensional aerodynamic models based on the potential flow are coupled with time domain simulations in the longitudinal plane. A special attention is paid to the explanation of the dynamic ground effect on both the sink and pitching motions. The influence of unsteady and quasi-steady forces on the dynamic ground effects and the craft motion is analyzed for different heights of flight.


---

## Body

## 1. Introduction

WIG (wing-in-ground) crafts are still considered as a very promising means of transport since they possess amphibian properties and can attain speeds which are unattainable for waterborne vessels. A very strong advantage of a WIG craft is the high lift to drag (L/D) ratio which theoretically can attain values over 30 and even more. For comparison, the L/D ratio of modern airplanes is less than twenty and the possibilities to improve it have been sufficiently exhausted. Unfortunately, to the knowledge of the authors, the L/D ratios over 20 have never been attained for real WIG configurations due to the following reasons. First, a very high L/D ratio is documented at extremely low clearance between the craft and the ground when the ratio of the clearance under the trailing edge h referred to as the mean aerodynamic chord h¯ is less than 0.05. In reality, the minimum height of flight near the free surface is limited by waves and safety considerations. Therefore, only a large craft can fly at small relative heights of flight and completely utilize the ground effect attaining a high L/D ratio. All crafts known to the authors both designed and tested have the cruise height of flight h¯ over 0.15. Second, since the fuselage of the WIG craft is designed to reduce the dynamic loadings during the takeoff and landing, its form is far from the aerodynamically optimal one. For instance, one utilizes very often the stepped bottom to reduce the drag in the takeoff mode. Third, to solve the problem of the longitudinal stability, most of the WIG crafts have big tail units located relatively high and far from the central main wing (see wing arrangements 1, 2, and 3 in Figure 1). The tail unit area of the best craft is varied between 22 and 30 percent of the main wing area. Such an arrangement requires the big fuselage which reduces the L/D ratio. Fourth, the majority of the WIG crafts built so far are small-sized ones and have relatively low speed U. To achieve the required lift, the cruise angle of attack is relatively large (from four to six degrees) and does not correspond to the optimal angle of attack at which the L/D ratio is maximal. Fifth, the main wing has usually small aspect ratio AR. It is well known from aerodynamics that the increase of the aspect ratio is the most efficient way to raise the L/D ratio. There is an opinion that the small aspect ratio is sufficient for WIG because the L/D ratio is high enough just due to the ground effect. Additional reasons speaking for small AR are the possibility to perform sharp turning maneuvers at big roll angles without water surface contact and increase of the efficiency of the power augmentation regime (PAR) (the most recent study, see, for instance, [1]). However, the moderate L/D ratio attained for the existing craft motivated designers to come back to the well-tried recipe to improve the L/D ratio by increase of AR. The last WIG Russian designs like SKB and marine passenger ekranoplan MPE designed by D. Synitsin (see the configuration 2 in Figure 1) in the Bureau of Alexeyev or SKP designed by Volkov and Ganin at KSRI have side wings of high aspect ratio located above the main wing plane by using V-angle either for the whole side wing (SKP) or for its part (SKB and MPE). Due to this arrangement change, the overall AR was increased from the usual values of two or three to five and higher. A strong limitation for AR increase is the raise of the derivative dh/dU which reduces the “binding” between the craft and ground [2]. This can cause jumps of the craft in a vertical direction due to wind gusts or thrust increase. With the growth of craft size many of these problems can be overcome and L/D can be sufficiently increased for ultra large WIG craft which are getting probable with the appearance of modern light materials.Figure 1
Typical wing arrangements of the wing in the ground effect craft [3]: (a) airplane type, (b) composite wing, (c) Lippish type, (d) tandem, (e) Volga type.
(a)
(b)
(c)
(d)
(e)At present, there is a large body of literature on the ground effect aerodynamics although some of the key works were published in various reports and are still being not available to readers. For instance, according to the author’s opinion, the most comprehensive and informative experimental study of the PAR regime was performed by the group of Prof. Epstein from the Central Aerohydrodynamic Institute (TsAGI). This group studied the effect of the wing planform, aspect ratio, distance between the engine nozzles and the leading edge of the wing, the blowing angles in horizontal and vertical directions, jet gas temperature, and the free surface (Froude number (Fn) effects) on the PAR efficiency. It was shown particularly that the optimal PAR efficiency determined as the ratio of the lift to thrust is attained at a small aspect ratio of the main wing around 0.7 (see Figure 48.14 in [3]). This explains the concept of the composite wing developed by many designers in Russia and Germany (see, for instance, configurations 2 and 3 in Figure 1 and in [4]). The wing is subdivided into two sections along the span, one of which is used for PAR whereas the second one is to generate lift in cruise regime. The free surface effect has a dramatical influence at Fn~0.7 similarly to common air cushion ships. At cruise Fn numbers, no influence of the free surface deformation on WIG aerodynamics was documented in the PAR regime.Most of physical phenomena in ground effect aerodynamics have been already described in various papers and monographs [5–7]. In this introduction, we would like to mention just three facts which are often not sufficiently interpreted in the literature. First, it concerns the behavior of the drag when the height h decreases. Very often, it is said that the drag also decreases due to reduction of the induced resistance. In reality, this effect depends on the wing arrangement. The drag can increase whereas the lift to drag ratio increases almost in all cases when the lift gets higher. Second, as correctly mentioned in many papers, the lift decreases when the Venturi effect comes into play. A less known fact is the lift reduction at high lift conditions. This effect which is documented, for instance, in the experimental study by Gadezky [8], can be explained from the analysis of the pressure distribution. The overpressure is strongly increased on the pressure side at small h whereas the under pressure slightly decreases on the suction side. At small angles of attack without flaps, the first effect is proved to be much stronger resulting in the substantial lift increase. For strong flap deflections or at high pitch angles and small clearances, the pressure under the wing cannot become higher than the stagnation one whereas the under pressure on the upper side keeps the tendency for reduction. As a result, the lift decreases when h goes to zero for high lift conditions in extreme ground effect. Third, there is a popular opinion that reversed delta wing arrangement 3 (Figure 1) proposed by Lippisch is the only one that has proven to be inherently stable in ground effect [9]. Stability of the WIG craft is determined by the distance between the aerodynamic centers in height Xh and in pitch Xϑ, position of the center of gravity, damping properties, and binding criterion.Numerical analysis shows that there exist a few other wing arrangements with common trapeze-like planform and tail unit which have similar stability properties as the reversed delta wing. The reader can easily prove it using the software Autowing [10] (http://www.lemos.uni-rostock.de/cfd-software/).The problems mentioned above are not the main barrier for a sustainable development of WIG technology. The amphibian property of WIG is the sufficient factor for their use in the transport system despite of the moderateL/D ratio. The main problem of WIG is the stability and the flight safety. Most detailed description of the dynamics, stability, and control of the WIG craft is given by Zhukov [11] and Diomidov [12] who made serious contribution to the design of all famous Soviet ekranoplans. An overview of some results can be found in the book of Rozhdestvensky [6] and in the SNAME handbook [3]. In fact, WIG crafts designed by Russian and German specialists were thoroughly proved using the linear stability theory and nonlinear motion simulations. Although all famous crafts satisfy the criteria of both static and dynamic stabilities, the trial tests and operation of WIG crafts have been accompanied by stability problems. A possible reason is the fact that the classical stability theory is valid for the case of small perturbations and it is based on the linear representation of forces which are strongly nonlinear near the ground. Avvakumov [13] introduced the term of the stability “in small” meaning the stability predicted by the linear stability analysis to differentiate it from the stability “on the whole.” The realistic perturbations cannot be considered as small ones, and the conclusions drawn from the linear stability theory are not sufficient to guarantee safe flight.As shown in our recent paper [14], the motion of the WIG craft which is optimal from the point of the classical stability theory including Irodov’s criterion [15], position of the center of gravity, and binding criterion can become unstable under relatively weak perturbations like wind gust. Therefore, the nonlinear simulation of the motion under perturbations is the inevitable part of the flight safety analysis to study the stability “on the whole.”Full 3DOF time domain simulation of the WIG motion coupled with the force determination is the task which can be solved using the modern CFD (computational fluid dynamics) codes. However, it requires big computational resources and for study of WIG dynamics with variation of craft and perturbations parameters in reasonable time is invisible. A more practical way is the decoupled consideration of dynamics and aerodynamics. Aerodynamic calculations are performed as preprocessing providing the coefficients which are then used for integration of the motion equations. Such a procedure is sufficiently based on the way of force representation in the form of Taylor or Fourier series on kinematic parameters. The set of kinematic parameters representing the WIG motion is still being a discussion subject. This problem is addressed in Section2.1 of the present paper. Representation (6) was widely used in the Russian WIG community. It is based on the assumption of linear representation of forces caused by unsteady effects. The coefficients of unsteady forces are usually determined under the assumption that the craft performs harmonic oscillations with small amplitude and small frequency. This raises the following questions:
(i)
Is this representation accurate enough to simulate nonlinear unsteady WIG dynamics in the ground effect
(ii)
How is calculation of the coefficients in this representation done
(iii)
Is the quasi-steady approach suitable for motion simulationsThese three questions are in the focus of the present paper. Before proceeding further, we need to give the definition of the unsteady forces. Under unsteady force, we imply the fraction of the total aerodynamic force arising due to change of kinematic parameters and caused by the following vortex wake effects:(i)
Creation of unsteady vortex wake in the form of vortices with dominating transversal orientation
(ii)
Change of the tip vortex circulation along its length
(iii)
Oscillations of the tip vortex in horizontal and vertical directionsThe second and the third effects which are typical only for three dimensional flows have influence on the upwind velocities and forces arising on the tail unit.Contrary to unsteady forces, the quasi-steady time-dependent forces are calculated using the assumption that the aerodynamics is fully determined by instantaneous values of kinematic parameters neglecting unsteady vortex wake effects.The motivation for this paper is the so-called dynamic ground effect mentioned in some recent papers (see, for instance, [16–18]). Surprisingly, this effect has never been explicitly discussed in the community involved in the development of WIG crafts. To understand this fact, this phenomenon is analyzed in Section 3.2 of this paper. The analysis is performed on the basis of the potential flow models from the following reasons. First, separate analysis of the quasi-steady and unsteady force components, necessary to interpret the dynamic ground effect, is only possible within the potential theory. Second, at moderate pitch angles and very high Reynolds numbers, the flow around the WIG is free of flow separations and influence of the viscosity is negligible. Third, as our viscous simulation experience shows [19], deformation of grids and their topology during the motion close to the ground makes the interpretations of the result difficult. It is necessary to separate the grid change effects from the real aerodynamic ones.


## 2. Mathematical Models

### 2.1. Force Representation in WIG Aerodynamics

The aerodynamic model is based on the representation of forces in a form of truncated Taylor series. This idea taken from aviation has sufficiently been modified by Prof. Treshkov in the late sixties. The original force coefficientC representation used in aviation is formulated as
(1)Cα,t=C0α0+∂C∂αα0⋅αt−α0+∂C∂α̇¯α0⋅α̇bU+⋯.Here,α is calculated in radians and α̇¯=α̇b/U. In the proximity to the ground, the force coefficient depends additionally on the height of flight h. As Treshkov noted, representation (1) is not convenient for WIG crafts because of ambiguous definition of the angle of attack. In aviation, the angle of attack can be produced either by the craft pitching ϑ or by the vertical motion h. Thus, it can be stated that
(2)α=ϑorα=−ḣU.Therefore, there is no difference in the derivatives between two ways in the change of the angle of attack:(3)∂C∂ϑϑ0=−∂C∂ḣ¯α0,where ḣ¯=ḣt/U is the nondimensional vertical velocity.The motion is steady if(4)ϑ̇=0,h¨=0.In proximity to the ground, the vertical motion is principally unsteady and it is not similar to the pitching rotation. Thus, relation (3) is not valid and representation (1) has to be extended for motion in the proximity to the ground. For aerodynamic coefficients, a linear representation can be used:
(5)Cϑ,h,t=C0ϑ0,h0+∂C∂ϑϑ0,h0⋅ϑt−ϑ0+∂C∂h¯ϑ0,h0⋅ht−h0b+∂C∂ϑ̇¯ϑ0,h0⋅ϑ̇tbU+∂C∂ḣ¯ϑ0,h0⋅ḣtU+⋯,where the angle of attack was excluded and replaced by the pitch angle to escape the ambiguous definition. ϑ̇¯ is the nondimensional angular velocity ϑ̇tb/U. Since the steady force fraction described by the second and fourth terms in (5) is strongly nonlinear, it is used further in the nonlinear form:
(6)Cϑ,h,t=C0ϑ,h+∂C∂ϑ̇¯ϑ,h⋅ϑ̇tbU+∂C∂ḣ¯ϑ,h⋅ḣtU+⋯.The last two terms on the r.h.s. in (6) contain both quasi-steady and unsteady fractions of the force.For the study of unsteady force effects, we use a hierarchy of two models. In the next subsection, we present the coupled two-dimensional aerodynamic and motion models with two degrees of freedom based on the potential formulation. Since the resistance and induced drag are not captured, the change of the motion speed is neglected. To take the three-dimensional effects into account, we use the nonlinear vortex lattice method to calculate the coefficients in representation (6) which are then used for 3DOF motion simulation.In what follows, for the sake of brevity, we use designationsCl,mϑ, Cl,mh, Cl,mḣ, Cl,mϑ̇, and Cl,mh¨ meaning the derivatives on the pitch angle in radians and nondimensional parameters h¯, ḣ¯, ϑ̇¯, and h¨¯. Determination of these derivatives discussed in this paper is of a big importance because they are used as input parameters for the design of WIG control systems (see details in [11, 12]).


### 2.2. Two-Dimensional Coupled Aerodynamic and Motion Models with Two Degrees of Freedom (2D+2DOF Model)

Simulation of WIG motion coupled with the simultaneous determination of aerodynamic forces is a laborious problem requiring big computational resources. If the primary aim of the study is the qualitative analysis, it is worth to work first with simplified fast models retaining main typical aerodynamic features and containing minimum free parameters typical for CFD (e.g., grid and turbulence model). In this work, we studied influence of unsteady aerodynamic effects on the WIG dynamics using the two-dimensional potential flow model incorporated into two degrees of freedom motion simulation. The lifting system consists of the main airfoil with the chord length ofb=1.0 and tail unit with the chord of 0.25 (see Figure 2). Both lifting elements have the same setup angle of attack. The distance between the trailing edge of the airfoil and the leading edge of the tail unit in the horizontal direction is 0.65 and 0.5 in the vertical direction. The flow is assumed to be inviscid and has potential. The no-penetration condition is enforced on airfoils. The ground effect is modeled using the mirror image. The condition of zero circulation around contours C1 and C2 (see Figure 2) is satisfied according to the Thomson theorem. The airfoils and vortex wake are modeled using the method of discrete vortices [20]. The vortex sheet is represented by point vortices both on airfoils and in the vortex wake. The no penetration condition is satisfied at collocation points located between point vortices on the profile. The first free point vortex in the wake lies at the line tangential to the airfoil. The trailing edge is exactly in the middle between the first free wake vortex and the last point vortex closest to the trailing edge on the profile. Such an arrangement of point vortices and collocation points secures fulfillment of the Kutta condition. Motion of free vortices in the wake is determined from the equations of fluid particle trajectories.Figure 2
Sketch of the lifting system, the coordinate system, and contoursC1 and C2.The pressure is calculated from the Bernoulli equation written for unsteady potential flow in the reference system moving with the speed of airfoilU∗t:
(7)∂Φ∂t+Ua22−U→aU→∗+pρ=paρ,where U→a=U→r+U→∗ is the absolute velocity, U→r is the relative velocity, and Φ is the potential. Introducing the nondimensional quantities φ=Φ/Ub, ua=Ua/U, u∗=U∗/U, τ=tU/b, and Cp=p−pa/ρU2/2, one obtains after some transformations the difference of the pressure coefficient between the pressure and the suction sides
(8)ΔCp=Cp−−Cp+=−2∂∂τφ−−φ+−ur−2−ur+2.The relative velocities on the suction and pressure sides are(9)ur−=u00+γ2,ur+=u00−γ2,where γ is the strength of the vorticity sheet and u00 is the direct value of the velocity on the airfoil (see Figure 2). Introducing circulation Γ=φ+−φ−, finally, we have
(10)ΔCp=Cp−−Cp+=−2γu00−∂Γ∂τ.IntegratingΔCp and x/bΔCp along the profile, one obtains the force and moment coefficients. The two-dimensional code was successfully verified by comparison with analytical data for a single airfoil at h=∞.The motion is calculated in the semiconnected coordinate systemxoy (see Figure 2). The horizontal speed is assumed to be constant. Two degrees of freedom motion equations read
(11)my¨=FNcosϑ−mg,Jzzϑ¨=Mz,where FN is the normal force, Mz is the pitching moment, and y is the distance between the center of gravity and the ground. In the nondimensional form, the system reads
(12)y¨¯=d2y¯dτ2=μCNcosϑ−Fn−2,ϑ¨¯=d2ϑdτ2=izCm,where μ=ρb2L/2m, Fn−2=gb/U2, iz=ρb4L/2Jzz, CN=FN/ρU2/2bL, and Cm=Mz/ρU2/2b2L. Here, L=1 is the wing span size.This model is nonlinear with respect to the vortex wake influence and is applicable for arbitrary nonlinear dynamics with the only restriction of the constant horizontal speed.


### 2.3. Simulation of Motion on the Basis of Aerodynamic Coefficients Obtained from the Vortex Lattice Method (3D+3DOF Model)

Within this model, the aerodynamic coefficients are precomputed using the vortex lattice method [10] and then saved in the form of lookup tables. The data from the tables approximated with B-splines are used then for motion simulations. The vortex lattice method assumes the potential inviscid flow. The ground effect is modeled by the mirror image method. The lifting surfaces are modeled by vortex sheets represented by a set of M discrete horseshoe vortices with unknown circulations Γk determined from the no-penetration condition transformed to a system of linear equations with respect to Γk:
(13)∑kMΓkwjknj=Fjnj,j=1,M,where wjk is the velocity induced by kth horseshoe vortex with unit circulation at the jth collocation point and nj and Fj are the normal vector and the vector of incident flow speed, respectively, at the jth point. All time-dependent quantities are supposed to be harmonically changing with small frequencies and small amplitudes. Γkt, wjkt, and Fjt can be represented in the form of truncated Fourier series:
(14)Γkt=Γk0+∑iNΓkqiqi+Γkq̇iq̇i+Oωi2,qi∗2,ωiqi∗2,where qi=qi∗cosωit is the ith kinematic parameter changing with the small amplitude qi∗ and frequency ωi,q̇i=dqi/dt. The resting small terms of orders higher than one are neglected. In the limit of small frequencies and amplitudes, system (13) is reduced to three systems of linear equations solved sequentially:
(15)∑kMΓk0wjk0nj=Fj0nj,(16)∑kMΓkqwjk0nj=Fjqnj−∑kMΓk0wjkqnj,(17)∑kMΓkq̇wjk0nj=Fjq̇nj−∑kMΓk0wjkq̇nj−∑kMΓkqwjkunsteadynj.The first system (15) describes the steady problem and the second one (16) provides the derivatives on kinematic parameters qi whereas unsteady derivatives on q̇i are calculated from the last system (17). System (16) is independent of (17) only in the limit of small frequencies. Function w is a pure geometrical function depending on coordinates of horseshoe vortices. For surge, heave, and pitch motions, the latter does not depend on q̇i, i.e., wjkq̇=0. Each discrete vortex element produces unsteady vortices to satisfy the theorem of Thomson which is fulfilled for each element separately. As a result, the vortex sheet shed from each horseshoe element moves along the lifting surface and propagates further downstream in the wake. The intensity of such a vortex sheet is proportional to Γkq. The last term in (17) describes the contribution of the unsteady vortex sheet to the normal component of velocity. Within this work, the position of all free vortices is prescribed. They oscillate along with the lifting surfaces as rigid constructions connected to lifting surfaces. The longitudinal vortices are straight lines being tangential to local chords. The transversal vortices of the vortex sheet are parallel to elements which they shed from. In summary, the aerodynamic model is linearized with respect to the angle of attack, wake form influence, and airfoil thickness effect. However, the following very important factor of the nonlinearity in the ground effect is considered, namely, the no penetration condition is satisfied at the actual position of the lifting surface at each time moment not being projected to any plane parallel to the ground as it is the case in fully linear wing theories. Consideration of instantaneous orientation of lifting surfaces with respect to the ground is the main source of nonlinearity of aerodynamic coefficients on the height of flight and pitch angle. The importance of this factor is discussed in [10]. Once Γ and its derivatives are known, the forces and moments are calculated from the Joukowsky theorem. Further details can be found in [10].Although representation (6) is usually used in the frequency domain simulation formulation, it is utilized here within the nonlinear three DOF time domain simulation based on the system of three ordinary differential equations:
(18)mU̇xcg=T−ρUa2S2Cd,mU̇ycg=ρUa2S2Cl−mg,Jzzω̇z=ρUa2Sb2Cm,where ḣcg=Uycg, ϑ̇=ωz, Uy=Uycg−wyt, and Ua=Uxcg−wxt. Index cg denotes the center of gravity. The possible lift component and pitching moment due to the propulsion are neglected. Aerodynamic coefficients Cl, Cd, and Cm are represented in the form of (6).The Vortex lattice method was implemented in the software Autowing [10] (http://www.lemos.uni-rostock.de/cfd-software/). Description of the method and its thorough validation are given by Kornev and Matveev [2], Kornev and Treshkov [10], and Benedict et al. [21].The3D+3DOF model described above was used to simulate the dynamics of a WIG configuration of the Lippisch type (Figure 3). The fuselage and the side rudders were neglected. The configuration consists of the central wing which models the aerodynamic influence of the fuselage, main wing, tail unit, endplates, and winglets (Table 1). All sizes are referred to the root chord of the central wing b=3.85m. All lifting elements have the profile Clark Y. The endplate is modeled as the thin plate. The Venturi effect was not documented at pitch angles larger than minus three degrees at the minimal height of flight of 0.1 m. The mass of the craft is 1.7 tons when the craft flies with a zero pitch angle at height h=0.5m. The inertia moment is 4694 kg·m2 for these conditions. The craft is faultless from the point of the classical stability analysis. It possesses both static stability and dynamic stability up to the height of flight of 2.5 meters at a zero pitch angle. At the height of 0.5 m, Irodov’s criterion referred to b is about 0.13. The speed of the craft is U=41.7m/s if not explicitly mentioned further. The binding criterion (see [2, 11])
(19)dhdU=2UCl∂Cl/∂h¯XϑXh−Xϑ,calculated at this speed and a height of flight of 0.5 m, is 0.05 and can be considered as quite satisfactory. In contrast to CFD viscous methods, the VLM approach possesses a fast convergence reaching the final result at a few thousands of panels at heights of flight larger than h/b≥0.1. The total number of vortex elements utilized for the calculations of the hypothetical WIG craft was around two thousand using square panels.Figure 3
Hypothetical WIG craft used for the study.Table 1
Main geometrical characteristics of the WIG craft used in the calculations.
Part
Chord
Aspect ratio
Taper ratio
Sweep angle
Pitch angle
V-angle
Central wing
1.0
0.3
1.0
0°
6°
0°
Main wing
1.0
0.5
0.3
−20°
6°
−12°
Tail unit
0.2
4.0
1.0
0°
2°
0°
Winglet
0.2
1.4
0.7
25°
5°
35°


## 2.1. Force Representation in WIG Aerodynamics

The aerodynamic model is based on the representation of forces in a form of truncated Taylor series. This idea taken from aviation has sufficiently been modified by Prof. Treshkov in the late sixties. The original force coefficientC representation used in aviation is formulated as
(1)Cα,t=C0α0+∂C∂αα0⋅αt−α0+∂C∂α̇¯α0⋅α̇bU+⋯.Here,α is calculated in radians and α̇¯=α̇b/U. In the proximity to the ground, the force coefficient depends additionally on the height of flight h. As Treshkov noted, representation (1) is not convenient for WIG crafts because of ambiguous definition of the angle of attack. In aviation, the angle of attack can be produced either by the craft pitching ϑ or by the vertical motion h. Thus, it can be stated that
(2)α=ϑorα=−ḣU.Therefore, there is no difference in the derivatives between two ways in the change of the angle of attack:(3)∂C∂ϑϑ0=−∂C∂ḣ¯α0,where ḣ¯=ḣt/U is the nondimensional vertical velocity.The motion is steady if(4)ϑ̇=0,h¨=0.In proximity to the ground, the vertical motion is principally unsteady and it is not similar to the pitching rotation. Thus, relation (3) is not valid and representation (1) has to be extended for motion in the proximity to the ground. For aerodynamic coefficients, a linear representation can be used:
(5)Cϑ,h,t=C0ϑ0,h0+∂C∂ϑϑ0,h0⋅ϑt−ϑ0+∂C∂h¯ϑ0,h0⋅ht−h0b+∂C∂ϑ̇¯ϑ0,h0⋅ϑ̇tbU+∂C∂ḣ¯ϑ0,h0⋅ḣtU+⋯,where the angle of attack was excluded and replaced by the pitch angle to escape the ambiguous definition. ϑ̇¯ is the nondimensional angular velocity ϑ̇tb/U. Since the steady force fraction described by the second and fourth terms in (5) is strongly nonlinear, it is used further in the nonlinear form:
(6)Cϑ,h,t=C0ϑ,h+∂C∂ϑ̇¯ϑ,h⋅ϑ̇tbU+∂C∂ḣ¯ϑ,h⋅ḣtU+⋯.The last two terms on the r.h.s. in (6) contain both quasi-steady and unsteady fractions of the force.For the study of unsteady force effects, we use a hierarchy of two models. In the next subsection, we present the coupled two-dimensional aerodynamic and motion models with two degrees of freedom based on the potential formulation. Since the resistance and induced drag are not captured, the change of the motion speed is neglected. To take the three-dimensional effects into account, we use the nonlinear vortex lattice method to calculate the coefficients in representation (6) which are then used for 3DOF motion simulation.In what follows, for the sake of brevity, we use designationsCl,mϑ, Cl,mh, Cl,mḣ, Cl,mϑ̇, and Cl,mh¨ meaning the derivatives on the pitch angle in radians and nondimensional parameters h¯, ḣ¯, ϑ̇¯, and h¨¯. Determination of these derivatives discussed in this paper is of a big importance because they are used as input parameters for the design of WIG control systems (see details in [11, 12]).


## 2.2. Two-Dimensional Coupled Aerodynamic and Motion Models with Two Degrees of Freedom (2D+2DOF Model)

Simulation of WIG motion coupled with the simultaneous determination of aerodynamic forces is a laborious problem requiring big computational resources. If the primary aim of the study is the qualitative analysis, it is worth to work first with simplified fast models retaining main typical aerodynamic features and containing minimum free parameters typical for CFD (e.g., grid and turbulence model). In this work, we studied influence of unsteady aerodynamic effects on the WIG dynamics using the two-dimensional potential flow model incorporated into two degrees of freedom motion simulation. The lifting system consists of the main airfoil with the chord length ofb=1.0 and tail unit with the chord of 0.25 (see Figure 2). Both lifting elements have the same setup angle of attack. The distance between the trailing edge of the airfoil and the leading edge of the tail unit in the horizontal direction is 0.65 and 0.5 in the vertical direction. The flow is assumed to be inviscid and has potential. The no-penetration condition is enforced on airfoils. The ground effect is modeled using the mirror image. The condition of zero circulation around contours C1 and C2 (see Figure 2) is satisfied according to the Thomson theorem. The airfoils and vortex wake are modeled using the method of discrete vortices [20]. The vortex sheet is represented by point vortices both on airfoils and in the vortex wake. The no penetration condition is satisfied at collocation points located between point vortices on the profile. The first free point vortex in the wake lies at the line tangential to the airfoil. The trailing edge is exactly in the middle between the first free wake vortex and the last point vortex closest to the trailing edge on the profile. Such an arrangement of point vortices and collocation points secures fulfillment of the Kutta condition. Motion of free vortices in the wake is determined from the equations of fluid particle trajectories.Figure 2
Sketch of the lifting system, the coordinate system, and contoursC1 and C2.The pressure is calculated from the Bernoulli equation written for unsteady potential flow in the reference system moving with the speed of airfoilU∗t:
(7)∂Φ∂t+Ua22−U→aU→∗+pρ=paρ,where U→a=U→r+U→∗ is the absolute velocity, U→r is the relative velocity, and Φ is the potential. Introducing the nondimensional quantities φ=Φ/Ub, ua=Ua/U, u∗=U∗/U, τ=tU/b, and Cp=p−pa/ρU2/2, one obtains after some transformations the difference of the pressure coefficient between the pressure and the suction sides
(8)ΔCp=Cp−−Cp+=−2∂∂τφ−−φ+−ur−2−ur+2.The relative velocities on the suction and pressure sides are(9)ur−=u00+γ2,ur+=u00−γ2,where γ is the strength of the vorticity sheet and u00 is the direct value of the velocity on the airfoil (see Figure 2). Introducing circulation Γ=φ+−φ−, finally, we have
(10)ΔCp=Cp−−Cp+=−2γu00−∂Γ∂τ.IntegratingΔCp and x/bΔCp along the profile, one obtains the force and moment coefficients. The two-dimensional code was successfully verified by comparison with analytical data for a single airfoil at h=∞.The motion is calculated in the semiconnected coordinate systemxoy (see Figure 2). The horizontal speed is assumed to be constant. Two degrees of freedom motion equations read
(11)my¨=FNcosϑ−mg,Jzzϑ¨=Mz,where FN is the normal force, Mz is the pitching moment, and y is the distance between the center of gravity and the ground. In the nondimensional form, the system reads
(12)y¨¯=d2y¯dτ2=μCNcosϑ−Fn−2,ϑ¨¯=d2ϑdτ2=izCm,where μ=ρb2L/2m, Fn−2=gb/U2, iz=ρb4L/2Jzz, CN=FN/ρU2/2bL, and Cm=Mz/ρU2/2b2L. Here, L=1 is the wing span size.This model is nonlinear with respect to the vortex wake influence and is applicable for arbitrary nonlinear dynamics with the only restriction of the constant horizontal speed.


## 2.3. Simulation of Motion on the Basis of Aerodynamic Coefficients Obtained from the Vortex Lattice Method (3D+3DOF Model)

Within this model, the aerodynamic coefficients are precomputed using the vortex lattice method [10] and then saved in the form of lookup tables. The data from the tables approximated with B-splines are used then for motion simulations. The vortex lattice method assumes the potential inviscid flow. The ground effect is modeled by the mirror image method. The lifting surfaces are modeled by vortex sheets represented by a set of M discrete horseshoe vortices with unknown circulations Γk determined from the no-penetration condition transformed to a system of linear equations with respect to Γk:
(13)∑kMΓkwjknj=Fjnj,j=1,M,where wjk is the velocity induced by kth horseshoe vortex with unit circulation at the jth collocation point and nj and Fj are the normal vector and the vector of incident flow speed, respectively, at the jth point. All time-dependent quantities are supposed to be harmonically changing with small frequencies and small amplitudes. Γkt, wjkt, and Fjt can be represented in the form of truncated Fourier series:
(14)Γkt=Γk0+∑iNΓkqiqi+Γkq̇iq̇i+Oωi2,qi∗2,ωiqi∗2,where qi=qi∗cosωit is the ith kinematic parameter changing with the small amplitude qi∗ and frequency ωi,q̇i=dqi/dt. The resting small terms of orders higher than one are neglected. In the limit of small frequencies and amplitudes, system (13) is reduced to three systems of linear equations solved sequentially:
(15)∑kMΓk0wjk0nj=Fj0nj,(16)∑kMΓkqwjk0nj=Fjqnj−∑kMΓk0wjkqnj,(17)∑kMΓkq̇wjk0nj=Fjq̇nj−∑kMΓk0wjkq̇nj−∑kMΓkqwjkunsteadynj.The first system (15) describes the steady problem and the second one (16) provides the derivatives on kinematic parameters qi whereas unsteady derivatives on q̇i are calculated from the last system (17). System (16) is independent of (17) only in the limit of small frequencies. Function w is a pure geometrical function depending on coordinates of horseshoe vortices. For surge, heave, and pitch motions, the latter does not depend on q̇i, i.e., wjkq̇=0. Each discrete vortex element produces unsteady vortices to satisfy the theorem of Thomson which is fulfilled for each element separately. As a result, the vortex sheet shed from each horseshoe element moves along the lifting surface and propagates further downstream in the wake. The intensity of such a vortex sheet is proportional to Γkq. The last term in (17) describes the contribution of the unsteady vortex sheet to the normal component of velocity. Within this work, the position of all free vortices is prescribed. They oscillate along with the lifting surfaces as rigid constructions connected to lifting surfaces. The longitudinal vortices are straight lines being tangential to local chords. The transversal vortices of the vortex sheet are parallel to elements which they shed from. In summary, the aerodynamic model is linearized with respect to the angle of attack, wake form influence, and airfoil thickness effect. However, the following very important factor of the nonlinearity in the ground effect is considered, namely, the no penetration condition is satisfied at the actual position of the lifting surface at each time moment not being projected to any plane parallel to the ground as it is the case in fully linear wing theories. Consideration of instantaneous orientation of lifting surfaces with respect to the ground is the main source of nonlinearity of aerodynamic coefficients on the height of flight and pitch angle. The importance of this factor is discussed in [10]. Once Γ and its derivatives are known, the forces and moments are calculated from the Joukowsky theorem. Further details can be found in [10].Although representation (6) is usually used in the frequency domain simulation formulation, it is utilized here within the nonlinear three DOF time domain simulation based on the system of three ordinary differential equations:
(18)mU̇xcg=T−ρUa2S2Cd,mU̇ycg=ρUa2S2Cl−mg,Jzzω̇z=ρUa2Sb2Cm,where ḣcg=Uycg, ϑ̇=ωz, Uy=Uycg−wyt, and Ua=Uxcg−wxt. Index cg denotes the center of gravity. The possible lift component and pitching moment due to the propulsion are neglected. Aerodynamic coefficients Cl, Cd, and Cm are represented in the form of (6).The Vortex lattice method was implemented in the software Autowing [10] (http://www.lemos.uni-rostock.de/cfd-software/). Description of the method and its thorough validation are given by Kornev and Matveev [2], Kornev and Treshkov [10], and Benedict et al. [21].The3D+3DOF model described above was used to simulate the dynamics of a WIG configuration of the Lippisch type (Figure 3). The fuselage and the side rudders were neglected. The configuration consists of the central wing which models the aerodynamic influence of the fuselage, main wing, tail unit, endplates, and winglets (Table 1). All sizes are referred to the root chord of the central wing b=3.85m. All lifting elements have the profile Clark Y. The endplate is modeled as the thin plate. The Venturi effect was not documented at pitch angles larger than minus three degrees at the minimal height of flight of 0.1 m. The mass of the craft is 1.7 tons when the craft flies with a zero pitch angle at height h=0.5m. The inertia moment is 4694 kg·m2 for these conditions. The craft is faultless from the point of the classical stability analysis. It possesses both static stability and dynamic stability up to the height of flight of 2.5 meters at a zero pitch angle. At the height of 0.5 m, Irodov’s criterion referred to b is about 0.13. The speed of the craft is U=41.7m/s if not explicitly mentioned further. The binding criterion (see [2, 11])
(19)dhdU=2UCl∂Cl/∂h¯XϑXh−Xϑ,calculated at this speed and a height of flight of 0.5 m, is 0.05 and can be considered as quite satisfactory. In contrast to CFD viscous methods, the VLM approach possesses a fast convergence reaching the final result at a few thousands of panels at heights of flight larger than h/b≥0.1. The total number of vortex elements utilized for the calculations of the hypothetical WIG craft was around two thousand using square panels.Figure 3
Hypothetical WIG craft used for the study.Table 1
Main geometrical characteristics of the WIG craft used in the calculations.
Part
Chord
Aspect ratio
Taper ratio
Sweep angle
Pitch angle
V-angle
Central wing
1.0
0.3
1.0
0°
6°
0°
Main wing
1.0
0.5
0.3
−20°
6°
−12°
Tail unit
0.2
4.0
1.0
0°
2°
0°
Winglet
0.2
1.4
0.7
25°
5°
35°


## 3. Results

### 3.1. Validation Tests

To the author’s knowledge, there are no well-tried benchmark tests proven in different independent experiments for validation of numerical simulations of wing arrangements in the close proximity to the ground. Very often, the available data raise the doubts at small heights of flight because of inaccurate ground modeling especially if a simple plate without any suction of the plate boundary layer is utilized. This concerns both steady and unsteady aerodynamic parameters.In this section, we evaluate the accuracy of the potential Vortex lattice method (VLM) by comparison of Autowing simulation data with available measurements and CFD calculations using the OpenFOAM code. Figure4 presents the lift coefficient of a single wing with the aspect ratio of 1.5 and the profile NACA 6409 for two nondimensional heights of flight h/b=0.05 and 0.3. The OpenFOAM simulations (see [19]) were performed using the limitedLinear scheme with the mix 90% CDS +10%UDS for the convection term. The Autowing computations show a very strong convergence when the number of vortex elements grows. The results were obtained with 60 panels along the span and 40 along the chord, i.e., with square panels. The ground is modeled by the mirror image method in VLM and the moving wall in OpenFOAM. As seen, the results of simulations are in very good agreement with measurements [22] for the moderate height h/b=0.3. For very small h/b=0.05, CFD shows a good accuracy whereas VLM results deviate sufficiently from experiments especially at moderate pitch angles. Higher lift in VLM can be explained by insufficient modeling of the Venturi effect. Indeed, the thickness influence in VLM is partly taken into account by the sources located on the mean wing surface. This approach models properly the potential displacement effect at moderate heights. At extremely low heights, the profile form should be modeled more accurately. The second reason of the overestimation of the lift coefficient is the neglect of viscosity in the potential theory. As mentioned in Introduction, the heights of h¯=0.05 are unrealistic for WIG crafts because of wave limitations. Our experience shows that at more realistic heights larger than 0.15, such a thickness model is appropriate.Figure 4
Lift coefficient of the wing with the aspect ratio of 1.5 and profile NACA 6409 versus the angle of attack (or the pitch angle). Comparison between VLM and OpenFOAM simulations and experiment [22].OpenFOAM calculations (dashed lines) are used to validate the VLM simulations (solid lines) for the configuration of the SeaFalcon WIG [23] with a schematic fuselage model (Figure 5). The results shown for different pitch angles ϑ are taken from [19]. It should be noted that the setup angle of the main wing is 5.5 deg, so that the resulting angle of the wing is 5.5 degrees larger. The VLM results are in good agreement with CFD for the lift coefficient (Figure 5) and for the Irodov criterion (Figure 6). The latter result is a convincing indication of the good VLM accuracy because the Irodov criterion depends on derivatives X=Cmh/Clh−Cmϑ/Clϑ which proper computation is a more challenging problem than determination of coefficients. The pitching moment coefficient is small (Figure 6) because the reference point used for its calculation was located close to the pressure center. Therefore, a small difference between VLM and CFD for the moment coefficient is not important.Figure 5
Half of the schematic configuration of the SeaFalcon WIG (a, b) and lift coefficient (c) versus the nondimensional height.
(a)
(b)
(c)Figure 6
Pitching moment coefficient (a) and the Irodov criterion (b) versus the nondimensional height.
(a)
(b)Validation results for the WIG craft Hoverwing 20 [24] presented in Figure 7 demonstrate an excellent agreement between VLM and measurements for the pitch angles of the craft of 0 and 2 degrees. The setup angle of the main wing is 6 degrees. The disagreement of approximately ten percent appears at ϑ=4 deg which corresponds to the pitch angle of 10 deg for the main wing. Additional measurements are necessary to clarify the reason for this growing discrepancy. Among others, it could be an improper modeling of the ground in experiments using plates without boundary layer suction. This results in a higher lift. Generally, in the past, any detailed validations of WIG craft simulations were seriously complicated by the lack of extensive experimental wind tunnel studies which were not available to small-design companies because of their high costs.Figure 7
Lift (a) and pitching moment coefficient (b) for the WIG craft Hoverwing 20 [24] versus height of flight.
(a)
(b)Validation for the unsteady forces is performed for a single wing with the aspect ratio of 3.0 and the profile NACA 0015 moving with prescribed harmonic vertical oscillationsh¯=h0+0.045sinωt with the dimensionless frequency p=ωc/V=0.25, where V=12 m/s is the wing speed. VLM results are compared with unsteady OpenFOAM simulations obtained by Andreas Groß (unpublished material obtained by the research group lead by the author) using the morphing grid technology and k−ωSST turbulence model. Forces in VLM unsteady simulations are calculated in form (6). The ground is modeled by the mirror image method in VLM and the moving wall in OpenFOAM.Figure8 illustrates a good agreement between two approaches at a moderate height of flight h0=0.3. The VLM provides the averaged lift coefficient of Cl0=0.768 whereas OpenFOAM Cl0=0.756. The time-dependent force part obtained from the VLM reads ΔClt=Clhht+Clḣḣ/V=−0.28ht−4.39ḣ/V. Closer to the ground at h0=0.15, the discrepancy between two methods increases. While the averaged values are in good agreement Cl0=0.824 for the VLM versus Cl0=0.838 for the OpenFOAM, there is a certain deviation between unsteady parts: ΔClt=Clhht+Clḣḣ/V=−0.44ht−4.91ḣ/V for the VLM and ΔClt=−0.88ht−5.13ḣt/V=Aht+Bḣt/V for the OpenFOAM. Since at small frequencies, B≈Clḣ and both B are close to each other (−5.13 versus −4.91), one can conclude that the accuracy of the linear representation of unsteady parts in (6) for wing oscillations with relatively big oscillation amplitude and moderate frequencies is acceptable. The increase of derivative Clḣ with decreasing height (−4.39 at h0=0.3 versus −4.91 at h0=0.15) is due to the dynamic ground effect discussed in the next chapter. The difference between unsteady parts is obviously due to the difference between A coefficients which include steady Clh and unsteady Clh¨ parts, i.e., A=Clh+ωc/V2Clh¨. At small frequencies, the second term ω2Clh¨ can be neglected. The analysis shows that the discrepancy between CFD and VLM arises due to a strong nonlinear change of the steady coefficient Cl resulting in a large derivative Clh. At the small height of flight h0=0.15 and amplitude a/b=0.045, the instantaneous height of flight becomes around 0.1. In this case, the accuracy of the VLM computations is degraded due to insufficient modeling of the profile thickness effect and neglect of viscosity as mentioned above for the steady single-wing test case. The data fluctuations in the CFD curve for h0=0.15 are of a pure numerical character, caused by the worsening of the mesh quality when the morphing grid technology is applied in close proximity to the ground.Figure 8
The unsteady time-dependent force partΔClt=Clhht+Clḣḣ/BV at h0=0.3 (a) and h0=0.15 (b). Dashed line: change of the flight height asinωt; solid line: VLM calculations using Autowing; filled squares: OpenFOAM CFD computations.
(a)
(b)Other validation samples for single-wing and WIG craft arrangements are presented by Kornev and Matveev [2], Kornev and Treshkov [10], and Benedict et al. [21].Concluding, the validation results presented above can be considered as satisfactory and the VLM model can be utilized for the study in which the primary aim is the determination of the qualitative effects of the ground effect aerodynamics.


### 3.2. Phenomenon of the Dynamic Ground Effect

One of the ground effect phenomena attracting the attention of experts in the latest time is the so-called dynamic ground effect. For the review of works on the dynamic ground effect, the reader is referred to papers [17, 18]. The essence of this effect is a substantial increase of the total lift force when WIG sinks near the ground. Outside of the ground effect, the lift increase is purely due to an additional angle of attack caused by the sink rate of the craft. In the ground effect, the lift is proved to be much stronger. Within force representation (6), the additional dynamic force caused by the sink rate is described by the term
(20)∂Cl∂ḣϑ,h⋅ḣU=ClḣḣU.The fact that derivativeClḣ gets much larger in the ground effect is not new. It is well known for a long time and can be found in some papers widely available to readers. For instance, in the monograph of Rozhdestvensky [6] (see Figures 3.14 and 3.15 in [6]), we find that this derivative on ḣ multiplied by the height of flight hClḣ is approximately 1.4 whereas the derivative on the vertical acceleration hClh¨ is 0.5 in the limit of a zero Strouhal number. Therefore, according to the linear theory based on the matched asymptotic expansion method, both derivatives become infinite when the craft is approaching the ground. In reality, it is restricted as shown in our paper [10] (see Figure 11 in [10]).Calculations using the VLM method for the Lippisch configuration presented in Figure9(a) illustrate three interesting facts. First, the force and moment derivatives on ḣ at large and moderate heights of flight are close to these on pitch angle ϑ (see Figure 9(a)). This means that the effect of the sink rate is almost the same as the effect of the pitch angle. Change of the force is almost due to additional angle of attack caused by the sink rate. At small heights of flight, there is a big discrepancy between two derivatives since the orientation of the craft with respect to the ground plays a significant role. This result is in accordance with studies [17, 18] in which the height range is subdivided into three regions. At high and intermediate heights, the effect of the sink rate is almost the same as that of the pitch angle. At small heights, the authors make compression work effects responsible for the increase of the unsteady force. The results in [17, 18] were obtained for a two-dimensional flow around an airfoil. The dynamic ground effect in the three-dimensional formulation was studied in [16].Figure 9
Influence of the ground effect on aerodynamic derivatives for the Lippisch configuration craft. Red lines: full unsteady calculations; black lines with circles: quasi-steady calculations.
(a)
(b)Second, the derivatives onḣ determined using the quasi-steady approach, i.e., the aerodynamic force is determined by the instantaneous value of ḣ and the unsteady wake contribution is zero wjkunsteady=0, are almost equal to these obtained using the unsteady vortex lattice method, i.e., with account for unsteady vortices. For this case, the equation for Γkḣ can be reduced to that similar to steady problem system (15):
(21)∑kMΓkḣwjk0nj=Fjḣnj.Therefore,Γkḣ for the flat wing can be expressed through Γk0 as Γkḣ=Γk0Fkḣnk/Fk0nk=Γk0tan−1ϑ. The obvious similarity between (15) and (21) points out that the reason for the so-called dynamic ground effect consisting in the increase of the lift at the sink rate is just the same as that in the increase of the lift in a steady flight, i.e., due to stagnation of the flow between the ground and the lifting surface. Unsteady aerodynamic effects connected to change of the lift and creation of free vortices play a minor role.Third, the dynamic ground effect is strongly pronounced in the pitching motion much more than in the sink one. For this case, the unsteady aerodynamic effects have a strong impact on the derivatives on the pitching velocityϑ (see Figure 9(b)).Concluding, the dynamic ground effect has been known in the ground effect aerodynamics for a long time although not discussed explicitly. For the simple sink rate motion, it is caused by the stagnation of the flow between the ground and the lifting surface.


### 3.3. Influence of Unsteady Aerodynamic Effects on Motion of a Lifting System Near the Ground

The ground effect has a substantial influence on mean forces and moments as well on their derivatives onh, ϑ, ḣ, and ϑ̇. Moreover, the unsteady aerodynamic effects increase when h→0 especially for derivatives on ϑ̇. The growth of unsteady derivatives in ground effect suggests an idea that their influence on motion and stability analysis gets stronger and the accuracy of their prediction cannot be satisfactory when using either the methods based on the assumption of harmonic change of parameters with small amplitudes and frequencies or a quasi-steady approach. Indeed, the unsteady part of forces in (6) depends on ḣ and ϑ̇ linearly although the ground effect is an essentially nonlinear aerodynamic phenomenon. To clarify these questions, we perform a study of the sensibility of motion parameters to unsteady forces obtained from different models at large and small heights of flight.The first calculations were performed using the2D+2DOF model which is applicable to arbitrary motion and does not use the assumption of small harmonic oscillations. The mass m of the schematic craft with b=1m shown in Figure 2 and the position of the center of gravity were chosen from the force and moment balance equation for each height. The inertia moment is taken as Jzz=m⋅10kg⋅m2. Simulations were performed at a constant speed U=40m/s and middle pitch angle ϑ0=5deg. The numbers of discrete vortices on the main airfoil are N1=100 and N2=25 on the tail unit. To get the stable solution with vortex sheet deformation, the time step should be proportional to the distance between discrete vortices on the airfoil, i.e., Δτa=b/N1=0.25b/N2=10−2 (see [20]). Unfortunately, this time step is too coarse to get the stable motion computation. Increase of N1 was impossible due to numerical instability problems originating from the fact that the evolution of thin vortex sheets belongs to the ill-posed mathematical problems. To secure the stability of the overall simulation, two time steps are utilized. The determination of aerodynamic forces and motion parameters is carried out with the step Δτ=10−3. Within the time step Δτa>Δτ, the circulation of the first wake free vortices closest to the trailing edges is changed whereas their positions remain constant with respect to the lifting system. The shedding of free vortices from the trailing edges into wake and motion of free vortices occurs with Δτa=10−2. This two-level time step procedure allows one to hold the distance between free vortices proportional to N1−1, to avoid the Kelvin Helmholtz instability and to reduce the number of free vortices in the wake. The lifting system moves under step-like wind perturbations
(22)wxU=−0.1,if20<τ≤30,0,else.Results are presented in Figure10. The red lines correspond to full nonlinear unsteady simulation, whereas the black lines with circles present the simulation with quasi-steady determination of forces, i.e., without the last term in equation (10).
(23)ΔCp=Cp−−Cp+=−2γu00.Figure 10
Trajectory of the schematic WIG craft under step-like wind perturbations. Black lines with circles correspond to the simulation with quasi-steady determination of forces.The most important conclusion from this analysis is that the difference between the unsteady aerodynamics formulation and quasi-steady one decreases drastically when the lifting system approaches the ground. This conclusion was reproduced in all simulations regardless ofm, b, U, N1, N2, and Δτ and the method of motion equation integration.The next motion simulation was carried out with the use of the3D+3DOF model. The perturbed motion was studied for two initial heights of flight 0.5 m and 3 m. The corresponding nondimensional values, referred to the mean aerodynamic chord are, respectively, 0.13 and 0.78. Motion equation (18) was integrated using the Runge-Kutta method of the 4th order with the time step of 10-4 s. The mass of the WIG craft and position of the center of gravity were chosen from the force and moment balance equation for each height. The trajectory of the WIG craft (see Figure 3) under sinusoidal wind perturbations wx=wx0sin2πt/T, where wx0=−2m/s and T=10s is shown in Figure 11. Figure 12 shows the trajectory of the WIG craft under step-like wind perturbations
(24)wx=−2m/s,if0<t≤4s,0,ift>4s.Figure 11
Trajectory of the WIG craft under sinusoidal wind perturbationswx=wx0sin2πt/T, where wx0=−2m/s and T=10s. Black lines with circles correspond to the simulation with quasi-steady determination of forces.
(a)
(b)Figure 12
Trajectory of the WIG craft under step-like wind perturbations. Black lines with circles correspond to the simulation with quasi-steady determination of forces.
(a)
(b)We come to the same conclusion as this with2D+2DOF simulation. For calculation of motion, the importance of consideration of unsteady effects in the calculation of aerodynamic forces goes down drastically when the height of flight decreases. Moreover, the quasi-steady approach and the force representation are proved to be sufficient for application near the ground. This sufficiently simplifies the nonlinear motion simulation used to estimate the stability of the WIG craft “on the whole” (see definition in Introduction) and to design the WIG control systems.The explanation of this rather surprising result is the increase of restoring forces close to the ground. The increase of unsteady effects, i.e., effects caused by the formation of free vortices and their dynamics, whenh→0, has no remarkable effect on WIG dynamics because of big restoring forces which also grow near the ground. The motion is mostly determined by instantaneous values of kinematic parameters. Far from the ground, in the area of weak stability, small differences in forces result in large deviation in trajectories. Note that the quasistatic approach ḣ=0,ϑ̇=0 is fully nonapplicable regardless of height of flight.The influence of unsteady effects on the position of the stability roots (see [4, 11]) is illustrated in Figure 13. Locus of the first three roots λi, i=1,3, was obtained for heights of 0.5 and 2.0 meters using quasi-steady and unsteady aerodynamic forces. The first root corresponds to the aperiodic motion Imλ1=0 whereas two others Reλ2=Reλ3, Imλ2=−Imλ3 correspond to the long-period oscillatory mode. The roots corresponding to the short-period fast-decaying mode are not shown. The difference between unsteady and quasi-steady approaches for the height of 2.0 m is negligible, whereas two roots corresponding to the long-period oscillatory mode for the height of 0.5 m are slightly different. Since −Reλ2,3 is relatively large for h=0.5, the oscillatory mode decays relatively quickly and the difference between unsteady and quasi-steady approaches is not remarkable in the direct time domain simulation of motion.Figure 13
Locus of three first roots of the WIG craft at two heights of flight obtained under unsteady and quasi-steady aerodynamic forces.


## 3.1. Validation Tests

To the author’s knowledge, there are no well-tried benchmark tests proven in different independent experiments for validation of numerical simulations of wing arrangements in the close proximity to the ground. Very often, the available data raise the doubts at small heights of flight because of inaccurate ground modeling especially if a simple plate without any suction of the plate boundary layer is utilized. This concerns both steady and unsteady aerodynamic parameters.In this section, we evaluate the accuracy of the potential Vortex lattice method (VLM) by comparison of Autowing simulation data with available measurements and CFD calculations using the OpenFOAM code. Figure4 presents the lift coefficient of a single wing with the aspect ratio of 1.5 and the profile NACA 6409 for two nondimensional heights of flight h/b=0.05 and 0.3. The OpenFOAM simulations (see [19]) were performed using the limitedLinear scheme with the mix 90% CDS +10%UDS for the convection term. The Autowing computations show a very strong convergence when the number of vortex elements grows. The results were obtained with 60 panels along the span and 40 along the chord, i.e., with square panels. The ground is modeled by the mirror image method in VLM and the moving wall in OpenFOAM. As seen, the results of simulations are in very good agreement with measurements [22] for the moderate height h/b=0.3. For very small h/b=0.05, CFD shows a good accuracy whereas VLM results deviate sufficiently from experiments especially at moderate pitch angles. Higher lift in VLM can be explained by insufficient modeling of the Venturi effect. Indeed, the thickness influence in VLM is partly taken into account by the sources located on the mean wing surface. This approach models properly the potential displacement effect at moderate heights. At extremely low heights, the profile form should be modeled more accurately. The second reason of the overestimation of the lift coefficient is the neglect of viscosity in the potential theory. As mentioned in Introduction, the heights of h¯=0.05 are unrealistic for WIG crafts because of wave limitations. Our experience shows that at more realistic heights larger than 0.15, such a thickness model is appropriate.Figure 4
Lift coefficient of the wing with the aspect ratio of 1.5 and profile NACA 6409 versus the angle of attack (or the pitch angle). Comparison between VLM and OpenFOAM simulations and experiment [22].OpenFOAM calculations (dashed lines) are used to validate the VLM simulations (solid lines) for the configuration of the SeaFalcon WIG [23] with a schematic fuselage model (Figure 5). The results shown for different pitch angles ϑ are taken from [19]. It should be noted that the setup angle of the main wing is 5.5 deg, so that the resulting angle of the wing is 5.5 degrees larger. The VLM results are in good agreement with CFD for the lift coefficient (Figure 5) and for the Irodov criterion (Figure 6). The latter result is a convincing indication of the good VLM accuracy because the Irodov criterion depends on derivatives X=Cmh/Clh−Cmϑ/Clϑ which proper computation is a more challenging problem than determination of coefficients. The pitching moment coefficient is small (Figure 6) because the reference point used for its calculation was located close to the pressure center. Therefore, a small difference between VLM and CFD for the moment coefficient is not important.Figure 5
Half of the schematic configuration of the SeaFalcon WIG (a, b) and lift coefficient (c) versus the nondimensional height.
(a)
(b)
(c)Figure 6
Pitching moment coefficient (a) and the Irodov criterion (b) versus the nondimensional height.
(a)
(b)Validation results for the WIG craft Hoverwing 20 [24] presented in Figure 7 demonstrate an excellent agreement between VLM and measurements for the pitch angles of the craft of 0 and 2 degrees. The setup angle of the main wing is 6 degrees. The disagreement of approximately ten percent appears at ϑ=4 deg which corresponds to the pitch angle of 10 deg for the main wing. Additional measurements are necessary to clarify the reason for this growing discrepancy. Among others, it could be an improper modeling of the ground in experiments using plates without boundary layer suction. This results in a higher lift. Generally, in the past, any detailed validations of WIG craft simulations were seriously complicated by the lack of extensive experimental wind tunnel studies which were not available to small-design companies because of their high costs.Figure 7
Lift (a) and pitching moment coefficient (b) for the WIG craft Hoverwing 20 [24] versus height of flight.
(a)
(b)Validation for the unsteady forces is performed for a single wing with the aspect ratio of 3.0 and the profile NACA 0015 moving with prescribed harmonic vertical oscillationsh¯=h0+0.045sinωt with the dimensionless frequency p=ωc/V=0.25, where V=12 m/s is the wing speed. VLM results are compared with unsteady OpenFOAM simulations obtained by Andreas Groß (unpublished material obtained by the research group lead by the author) using the morphing grid technology and k−ωSST turbulence model. Forces in VLM unsteady simulations are calculated in form (6). The ground is modeled by the mirror image method in VLM and the moving wall in OpenFOAM.Figure8 illustrates a good agreement between two approaches at a moderate height of flight h0=0.3. The VLM provides the averaged lift coefficient of Cl0=0.768 whereas OpenFOAM Cl0=0.756. The time-dependent force part obtained from the VLM reads ΔClt=Clhht+Clḣḣ/V=−0.28ht−4.39ḣ/V. Closer to the ground at h0=0.15, the discrepancy between two methods increases. While the averaged values are in good agreement Cl0=0.824 for the VLM versus Cl0=0.838 for the OpenFOAM, there is a certain deviation between unsteady parts: ΔClt=Clhht+Clḣḣ/V=−0.44ht−4.91ḣ/V for the VLM and ΔClt=−0.88ht−5.13ḣt/V=Aht+Bḣt/V for the OpenFOAM. Since at small frequencies, B≈Clḣ and both B are close to each other (−5.13 versus −4.91), one can conclude that the accuracy of the linear representation of unsteady parts in (6) for wing oscillations with relatively big oscillation amplitude and moderate frequencies is acceptable. The increase of derivative Clḣ with decreasing height (−4.39 at h0=0.3 versus −4.91 at h0=0.15) is due to the dynamic ground effect discussed in the next chapter. The difference between unsteady parts is obviously due to the difference between A coefficients which include steady Clh and unsteady Clh¨ parts, i.e., A=Clh+ωc/V2Clh¨. At small frequencies, the second term ω2Clh¨ can be neglected. The analysis shows that the discrepancy between CFD and VLM arises due to a strong nonlinear change of the steady coefficient Cl resulting in a large derivative Clh. At the small height of flight h0=0.15 and amplitude a/b=0.045, the instantaneous height of flight becomes around 0.1. In this case, the accuracy of the VLM computations is degraded due to insufficient modeling of the profile thickness effect and neglect of viscosity as mentioned above for the steady single-wing test case. The data fluctuations in the CFD curve for h0=0.15 are of a pure numerical character, caused by the worsening of the mesh quality when the morphing grid technology is applied in close proximity to the ground.Figure 8
The unsteady time-dependent force partΔClt=Clhht+Clḣḣ/BV at h0=0.3 (a) and h0=0.15 (b). Dashed line: change of the flight height asinωt; solid line: VLM calculations using Autowing; filled squares: OpenFOAM CFD computations.
(a)
(b)Other validation samples for single-wing and WIG craft arrangements are presented by Kornev and Matveev [2], Kornev and Treshkov [10], and Benedict et al. [21].Concluding, the validation results presented above can be considered as satisfactory and the VLM model can be utilized for the study in which the primary aim is the determination of the qualitative effects of the ground effect aerodynamics.


## 3.2. Phenomenon of the Dynamic Ground Effect

One of the ground effect phenomena attracting the attention of experts in the latest time is the so-called dynamic ground effect. For the review of works on the dynamic ground effect, the reader is referred to papers [17, 18]. The essence of this effect is a substantial increase of the total lift force when WIG sinks near the ground. Outside of the ground effect, the lift increase is purely due to an additional angle of attack caused by the sink rate of the craft. In the ground effect, the lift is proved to be much stronger. Within force representation (6), the additional dynamic force caused by the sink rate is described by the term
(20)∂Cl∂ḣϑ,h⋅ḣU=ClḣḣU.The fact that derivativeClḣ gets much larger in the ground effect is not new. It is well known for a long time and can be found in some papers widely available to readers. For instance, in the monograph of Rozhdestvensky [6] (see Figures 3.14 and 3.15 in [6]), we find that this derivative on ḣ multiplied by the height of flight hClḣ is approximately 1.4 whereas the derivative on the vertical acceleration hClh¨ is 0.5 in the limit of a zero Strouhal number. Therefore, according to the linear theory based on the matched asymptotic expansion method, both derivatives become infinite when the craft is approaching the ground. In reality, it is restricted as shown in our paper [10] (see Figure 11 in [10]).Calculations using the VLM method for the Lippisch configuration presented in Figure9(a) illustrate three interesting facts. First, the force and moment derivatives on ḣ at large and moderate heights of flight are close to these on pitch angle ϑ (see Figure 9(a)). This means that the effect of the sink rate is almost the same as the effect of the pitch angle. Change of the force is almost due to additional angle of attack caused by the sink rate. At small heights of flight, there is a big discrepancy between two derivatives since the orientation of the craft with respect to the ground plays a significant role. This result is in accordance with studies [17, 18] in which the height range is subdivided into three regions. At high and intermediate heights, the effect of the sink rate is almost the same as that of the pitch angle. At small heights, the authors make compression work effects responsible for the increase of the unsteady force. The results in [17, 18] were obtained for a two-dimensional flow around an airfoil. The dynamic ground effect in the three-dimensional formulation was studied in [16].Figure 9
Influence of the ground effect on aerodynamic derivatives for the Lippisch configuration craft. Red lines: full unsteady calculations; black lines with circles: quasi-steady calculations.
(a)
(b)Second, the derivatives onḣ determined using the quasi-steady approach, i.e., the aerodynamic force is determined by the instantaneous value of ḣ and the unsteady wake contribution is zero wjkunsteady=0, are almost equal to these obtained using the unsteady vortex lattice method, i.e., with account for unsteady vortices. For this case, the equation for Γkḣ can be reduced to that similar to steady problem system (15):
(21)∑kMΓkḣwjk0nj=Fjḣnj.Therefore,Γkḣ for the flat wing can be expressed through Γk0 as Γkḣ=Γk0Fkḣnk/Fk0nk=Γk0tan−1ϑ. The obvious similarity between (15) and (21) points out that the reason for the so-called dynamic ground effect consisting in the increase of the lift at the sink rate is just the same as that in the increase of the lift in a steady flight, i.e., due to stagnation of the flow between the ground and the lifting surface. Unsteady aerodynamic effects connected to change of the lift and creation of free vortices play a minor role.Third, the dynamic ground effect is strongly pronounced in the pitching motion much more than in the sink one. For this case, the unsteady aerodynamic effects have a strong impact on the derivatives on the pitching velocityϑ (see Figure 9(b)).Concluding, the dynamic ground effect has been known in the ground effect aerodynamics for a long time although not discussed explicitly. For the simple sink rate motion, it is caused by the stagnation of the flow between the ground and the lifting surface.


## 3.3. Influence of Unsteady Aerodynamic Effects on Motion of a Lifting System Near the Ground

The ground effect has a substantial influence on mean forces and moments as well on their derivatives onh, ϑ, ḣ, and ϑ̇. Moreover, the unsteady aerodynamic effects increase when h→0 especially for derivatives on ϑ̇. The growth of unsteady derivatives in ground effect suggests an idea that their influence on motion and stability analysis gets stronger and the accuracy of their prediction cannot be satisfactory when using either the methods based on the assumption of harmonic change of parameters with small amplitudes and frequencies or a quasi-steady approach. Indeed, the unsteady part of forces in (6) depends on ḣ and ϑ̇ linearly although the ground effect is an essentially nonlinear aerodynamic phenomenon. To clarify these questions, we perform a study of the sensibility of motion parameters to unsteady forces obtained from different models at large and small heights of flight.The first calculations were performed using the2D+2DOF model which is applicable to arbitrary motion and does not use the assumption of small harmonic oscillations. The mass m of the schematic craft with b=1m shown in Figure 2 and the position of the center of gravity were chosen from the force and moment balance equation for each height. The inertia moment is taken as Jzz=m⋅10kg⋅m2. Simulations were performed at a constant speed U=40m/s and middle pitch angle ϑ0=5deg. The numbers of discrete vortices on the main airfoil are N1=100 and N2=25 on the tail unit. To get the stable solution with vortex sheet deformation, the time step should be proportional to the distance between discrete vortices on the airfoil, i.e., Δτa=b/N1=0.25b/N2=10−2 (see [20]). Unfortunately, this time step is too coarse to get the stable motion computation. Increase of N1 was impossible due to numerical instability problems originating from the fact that the evolution of thin vortex sheets belongs to the ill-posed mathematical problems. To secure the stability of the overall simulation, two time steps are utilized. The determination of aerodynamic forces and motion parameters is carried out with the step Δτ=10−3. Within the time step Δτa>Δτ, the circulation of the first wake free vortices closest to the trailing edges is changed whereas their positions remain constant with respect to the lifting system. The shedding of free vortices from the trailing edges into wake and motion of free vortices occurs with Δτa=10−2. This two-level time step procedure allows one to hold the distance between free vortices proportional to N1−1, to avoid the Kelvin Helmholtz instability and to reduce the number of free vortices in the wake. The lifting system moves under step-like wind perturbations
(22)wxU=−0.1,if20<τ≤30,0,else.Results are presented in Figure10. The red lines correspond to full nonlinear unsteady simulation, whereas the black lines with circles present the simulation with quasi-steady determination of forces, i.e., without the last term in equation (10).
(23)ΔCp=Cp−−Cp+=−2γu00.Figure 10
Trajectory of the schematic WIG craft under step-like wind perturbations. Black lines with circles correspond to the simulation with quasi-steady determination of forces.The most important conclusion from this analysis is that the difference between the unsteady aerodynamics formulation and quasi-steady one decreases drastically when the lifting system approaches the ground. This conclusion was reproduced in all simulations regardless ofm, b, U, N1, N2, and Δτ and the method of motion equation integration.The next motion simulation was carried out with the use of the3D+3DOF model. The perturbed motion was studied for two initial heights of flight 0.5 m and 3 m. The corresponding nondimensional values, referred to the mean aerodynamic chord are, respectively, 0.13 and 0.78. Motion equation (18) was integrated using the Runge-Kutta method of the 4th order with the time step of 10-4 s. The mass of the WIG craft and position of the center of gravity were chosen from the force and moment balance equation for each height. The trajectory of the WIG craft (see Figure 3) under sinusoidal wind perturbations wx=wx0sin2πt/T, where wx0=−2m/s and T=10s is shown in Figure 11. Figure 12 shows the trajectory of the WIG craft under step-like wind perturbations
(24)wx=−2m/s,if0<t≤4s,0,ift>4s.Figure 11
Trajectory of the WIG craft under sinusoidal wind perturbationswx=wx0sin2πt/T, where wx0=−2m/s and T=10s. Black lines with circles correspond to the simulation with quasi-steady determination of forces.
(a)
(b)Figure 12
Trajectory of the WIG craft under step-like wind perturbations. Black lines with circles correspond to the simulation with quasi-steady determination of forces.
(a)
(b)We come to the same conclusion as this with2D+2DOF simulation. For calculation of motion, the importance of consideration of unsteady effects in the calculation of aerodynamic forces goes down drastically when the height of flight decreases. Moreover, the quasi-steady approach and the force representation are proved to be sufficient for application near the ground. This sufficiently simplifies the nonlinear motion simulation used to estimate the stability of the WIG craft “on the whole” (see definition in Introduction) and to design the WIG control systems.The explanation of this rather surprising result is the increase of restoring forces close to the ground. The increase of unsteady effects, i.e., effects caused by the formation of free vortices and their dynamics, whenh→0, has no remarkable effect on WIG dynamics because of big restoring forces which also grow near the ground. The motion is mostly determined by instantaneous values of kinematic parameters. Far from the ground, in the area of weak stability, small differences in forces result in large deviation in trajectories. Note that the quasistatic approach ḣ=0,ϑ̇=0 is fully nonapplicable regardless of height of flight.The influence of unsteady effects on the position of the stability roots (see [4, 11]) is illustrated in Figure 13. Locus of the first three roots λi, i=1,3, was obtained for heights of 0.5 and 2.0 meters using quasi-steady and unsteady aerodynamic forces. The first root corresponds to the aperiodic motion Imλ1=0 whereas two others Reλ2=Reλ3, Imλ2=−Imλ3 correspond to the long-period oscillatory mode. The roots corresponding to the short-period fast-decaying mode are not shown. The difference between unsteady and quasi-steady approaches for the height of 2.0 m is negligible, whereas two roots corresponding to the long-period oscillatory mode for the height of 0.5 m are slightly different. Since −Reλ2,3 is relatively large for h=0.5, the oscillatory mode decays relatively quickly and the difference between unsteady and quasi-steady approaches is not remarkable in the direct time domain simulation of motion.Figure 13
Locus of three first roots of the WIG craft at two heights of flight obtained under unsteady and quasi-steady aerodynamic forces.


## 4. Conclusion

Unsteady forces arising due to change of kinematic parameters and formation of the vortex wake have strong impact on the aerodynamical characteristics of lifting systems near the ground. They are responsible for the so-called dynamic ground effect consisting in a substantial increase of the total lift force when a wing sinks near the ground. Here, the lift ath=∞ is compared with that at any small h and the same sink rate. Strictly speaking, this effect has been known in the ground effect aerodynamics for a long time although not discussed explicitly. At high and intermediate heights, the effect of the sink rate is almost the same as that of the pitch angle change. At small heights, the orientation of the craft with respect to the ground plays a significant role and the unsteady lift due to the sink rate is much larger than this of the pitch angle change. The present analysis based on the potential flow model shows clearly that the reason for the dynamic ground effect due to the sink is just the same as the increase of the lift in a steady flight, i.e., due to stagnation of the flow between the ground and the lifting surface. Unsteady aerodynamic effects connected to the change of the lift and creation of free vortices play a minor role. The dynamic ground effect is strongly pronounced in the pitching motion much more than that in the sink one.Although the fraction of forces caused by unsteady effects increases near the ground, the influence of unsteady aerodynamic effects on the motion of a lifting system gets smaller when the height of flight decreases. The motion simulation can be predicted quite properly using the quasi-steady approach. The explanation of this rather surprising result is the increase of restoring forces close to the ground. The increase of unsteady effects, i.e., effects caused by the formation of free vortices and their dynamics, whenh→0, has no remarkable effect on WIG dynamics because of big restoring forces which also grow near the ground. The motion is mostly determined by instantaneous values of kinematic parameters. Far from the ground, in the area of weak stability, small differences in forces result in large deviation in trajectories.



---

*Source: 8351293-2019-05-06.xml*