General Solutions to the Navier-Stokes Equations for Incompressible Flow

Article information

J. Ocean Eng. Technol. 2024;38(5):315-324

Publication date (electronic) : 2024 October 18

doi : https://doi.org/10.26748/KSOE.2024.051

JangRyong Shin^1,

¹Engineer, Offshore structure design department, Hanwha Ocean Co., LTD, Geoje, Korea

Corresponding author JangRyong Shin: +82-55-680-5117, jrshin9797@hanwha.com

Received 2024 April 24; Revised 2024 July 19; Accepted 2024 September 25.

Abstract

1. Introduction

Water waves are considered one of the most fascinating natural phenomena in the marine environment. Waves in the open ocean are generally generated by winds. The transmission of energy from the wind to the water is facilitated by friction between these two fluid media. This energy is transmitted as a wave across the water surface. Consequently, friction plays a pivotal role in the circulation of fluids and is a crucial factor in wave generation. The energy transferred from the wind to waves diminishes as the wind intensity decreases, leading to the transformation of waves into swells that eventually dissipate. Therefore, friction also contributes to the dissipation of waves. Nevertheless, the mathematical description of water waves has been conducted under the assumption of irrotational, inviscid flow (Chaplin, 1980; Chappeleakimr, 1961; De, 1955; Dean, 1965; Fenton, 1988; Kishida and Sobey, 1988; Korteweg and de Vries, 1895; Rienecker and Fenton, 1981; Shin, 2016, 2019; Stokes, 1847; Stoke, 1880; Skjelbreia and Hendrickson, 1960; Vanden-Broeck and Schwartz., 1979) or under the assumption of rotational, inviscid flow (Rankine, 1863; Chen and Zou, 2019; Constantin, 2005; Dalrymple, 1974; Henry; 2008; Nwogu, 2009; Shin, 2022, 2023). Hence, these studies cannot account for the natural dissipation of waves. In a few exceptional cases, viscosity is considered in the linearized Navier-Stokes equations, e.g., water waves over a viscous mud bottom and seiche phenomena in coastal regions. The equation (Dean and Dalrymple, 1984) for long waves with bottom friction is analogous to that of a damped harmonic oscillator. The amplitude decreases exponentially with time, and the angular frequency also changes, providing an important clue for solving the Navier–Stokes (NS) equations. NS equations are a set of nonlinear partial differential equations that govern the flow of Newtonian fluids. These equations are valuable because of their ability to elucidate numerous physical phenomena. The structural design process commences by evaluating the external appearance of the design object and the loads on the surface generated by fluid flow. The NS equations are used for these determinations. Despite the wide range of practical applications, the existence of smooth solutions in three dimensions has not been proven definitively. Theoretical solutions remain limited for most cases, with only a few exceptions. Analytical approaches are deemed almost impossible for solving these problems. Thus, only numerical and experimental methods are allowed. This study applied a theoretical approach with tensor calculus to address these limitations. General solutions were introduced without considering the boundary conditions and initial conditions.

NS equations are the Euler equations with the addition of the gradient of the viscous stress tensor. The viscous stress tensor is proportional to the strain rate, making the gradient of viscous stress a linear function of the velocity field. The convective acceleration is the only nonlinear term included in the Euler equation. Shin (2022, 2023) reported the complete solutions to the Euler equations. They transformed the Euler equations into Bernoulli’s principle when the velocity field was described by the solutions of the Helmholtz equation. Using the methodology used in the studies conducted by Shin (2022, 2023) with the solution for a damped harmonic oscillator, the findings were extended to cover general three-dimensional viscous flow in the current study.

When the velocity field is described by solutions of the Helmholtz equation, the gradient of the viscous stress tensor for incompressible flow is proportional to the velocity field. The product of a solution to the Helmholtz equation and any time function is also a solution to the Helmholtz equation. By substituting the product into the NS equations, these equations can be decomposed into the Euler equations and a linear ordinary differential equation for the time function. Based on Shin’s (2022) research, this study solved the Euler equations for general flow. Some exact solutions to NS equations exist. The Taylor–Green vortex (Taylor and Green, 1937) and Antonio’s solution (Antuono, 2020) are two interesting examples. The Taylor–Green vortex was an unsteady flow of a decaying vortex in Cartesian coordinates. Antonio’s solution is considered a three-dimensional generalization of the classic two-dimensional Taylor–Green vortex. For verification, the Taylor–Green vortex and Antonio’s solutions are special cases of the solutions presented in this study. Substituting the results into the vorticity equations confirmed that the vorticity equations were satisfied.

This study adjusted Shin’s results (2022, 2023) to comply with NS equations. The process through which waves evolve into swells was elucidated. The modified wave theory was verified by comparing the decay modulus with the experimental data of Shiau and Rumer (1974) and the theoretical solution presented by Dean and Dalrymple (1984). The comparison results showed a high level of agreement.

2. Navier-Stokes Equations for Incompressible Flow

Tensor notation was used for simplicity, and the following rule was adopted when using tensor notation. Whenever the subscript appeared twice as a letter in a grouping, the expressions of this grouping were summed, with the repeated index taking on successively the values 1, 2, and 3. The non-repeated letter subscripts or numerical subscripts remain fixed during the process (Dym and Shames, 1973). The partial derivatives of a tensor are denoted using commas and indices as ∂()∂xi=(),i where x_i is the position vector; the conservation of mass is presented as follows:

A time derivative will be denoted with the ‘over-dot’ notation. Therefore, v̇_i is the local acceleration. ρ is the density of the fluid, which is a constant; p = p(x_i, t) is pressure field; μ is kinematic viscosity, which is a constant; f_i is body force per unit mass, which is a conservative force. The NS equations for an incompressible fluid are composed of Eqs. (1) and (2). The vorticity is defined as follows:

Taking the inner product of the velocity v_p and the deformation tensor v_q,p, the convective acceleration in Eq. (2) is expressed as

Replacing p = j, q = i, and i = k, Eq. (5) can be presented as follows because v_pΩ_i∊_ipq = v_pΩ_i∊_qip =−v_pΩ_i∊_qpi :

Substituting Eq. (6) into Eq. (2), the alternative form of the momentum equation is represented as follows:

Because f_i is a conservative force, it has potential energy. The first term on the right-hand side represents the outer product of the velocity and vorticity. Therefore, the term is zero when the vorticity and the velocity are parallel to each other.

3. Decomposition of the Navier-Stokes Equations

A definitive clue to solving the NS equations can be found in the solution of a single-degree-of-freedom damped harmonic oscillator or long waves with bottom friction (Dean and Dalrymple, 1984). In the oscillator, the damping force is proportional to the velocity ż of the object because it is due to fluid viscosity: −cż, where c is the damping coefficient. Newton’s second law for a damped harmonic oscillator is presented as follows:

Theorem 1. When the velocity fields are represented as a product of a time function and solutions of the Helmholtz equation, the NS equations can be decomposed into the Euler equations and an ordinary equation for the time function.

Proof. Substituting Eq. (10) into Eq. (1), the conservation of mass is presented as follows:

The gradient of viscous stress is presented as μv_i,jj = μbv_i = μbF (t)u_i, and the local acceleration is presented as follows:

Substituting Eq. (12) and μv_i,jj = μbF (t)u_i into Eq. (7), Eq. (7) is converted to the Euler momentum equation, provided that the following condition is satisfied:

Because u_i ≠ 0, the solution to Eq. (13) can be presented as follows:

The constant, b, must be negative because the flow must eventually dissipate. Eq. (14) shows similarities to the exponential decay of the amplitude, as expressed in Eq. (9). Eq. (7) is presented as follows:

Theorem 2. When the velocity field is a solution of the Helmholtz equation, then the vorticity is also a solution of the Helmholtz equation.

Proof. Ω_i,jj =(∊_ipqv_q,p )_,jj =(∊_ipqv_q,jj)_,p =(∊_ipqbv_q )_,p = b(∊_ipqv_q,p )= bΩ_i

As a result, the velocity field and the vorticity field were determined using the Helmholtz equation. The next step was to determine the pressure field using Eq. (15). “ 0 < F (t) < 1”, and it has no dimensions. Hence, it is considered a non-dimensional perturbation parameter. Using the perturbation method, Eq. (15) can be solved in the following chapters for two cases: u_i is steady flow and u_i is unsteady flow.

4. Pressure Field for Steady Flow

When the velocity field, u_i is independent of time, Eq. (10) is presented as follows:

Therefore, the local acceleration in Eq. (15) is zero. The pressure field is presented as the sum of the dynamic pressure and static pressure. The pressure field must be presented by a quadratic function with respect to F(t) because the dynamic pressure is proportional to the square of the velocity, and static pressure is proportional to the body force, as follows.

The first and second terms mean the static and dynamic pressures, respectively. Substituting Eqs. (16)–(17) into Eq. (15), the result is presented as

Eq. (18) is decomposed with two equations because Eq. (18) is valid for all t, and F (t) ≠ 0, using the perturbation method, as follows. The zero-order perturbation equation is presented as follows:

The second-order perturbation equation is presented as

Integrating Eq. (19), the potential energy of the body force is determined as follows:

The static pressure field is determined using Eq. (21). u_i is a solution of the Helmholtz equation. Hence, u_j,pp = bu_j. Therefore, the last term of Eq. (20) is presented as follows:

Using Eq. (4),

Using Eq. (11), u_p,pj = 0, and uj,pp=Ω¯r,p∊rpj. Substituting it into Eq. (22),

Where δ_ij is the Kronecker delta. Since divergence of curl is zero, we have Ω¯k,k=0. Therefore, Ω¯kΩ¯i,k=(Ω¯kΩ¯i),k. Let’s introduce a quantity, E,idef=(Ω¯kΩ¯i),k. Therefore, E=∫(Ω¯kΩ¯i),k dx_i is a type of kinetic energy. Because Ω¯k is also a solution of the Helmholtz equation and is determined by the velocity field, the energy, E can be easily determined. Substituting Eq. (24) into Eq. (20), the result is presented as follows.

Bernoulli’s principle can be expressed as Eq. (26) by integrating Eq. (25) with respect to x_i :

The dynamic pressure field is determined using Eq. (26). The second term represents the kinetic energy associated with irrotational motion. In contrast, the third term represents the kinetic energy related to the vortex. The total kinetic energy is the linear sum of the two terms because the constant, b is negative, which is consistent with the law of conservation of energy. The ratio, r, is defined as

Ω¯kΩ¯k≥2E≥0 because the kinetic energy is a positive-definite function and the constant, b, is negative. Therefore 0≤r≤1.

The energy E is referred to as an anti-vortex energy because it reduces the kinetic energy due to the vortex like frictional energy. Two cases of anti-vortex energy are considered as follows:

5. Pressure Field for Unsteady Flow

When the velocity field, u_i is unsteady flow, the velocity field can be expressed as Eq. (10). A typical example is progressive water waves (Shin, 2022, 2023). Using the moving coordinate system proposed by Dean (1965), an unsteady flow can be converted into a steady flow, as presented in the previous Chapter. The moving coordinate system is defined as

G_i (t) has a length unit. In a damped harmonic oscillator, angular frequency is modified with a damping ratio, “ ζ,” as shown in Eq. (9). Therefore, the function G_i (t) shows that viscosity also affects the angular frequency of water waves. The first term in Eq. (15) must be a quadratic term of F (t), it is assumed that the function G_i (t) is defined as

C_i,j = 0 because C_i is a constant. Hence,

Using Eqs. (29) and (30),

Using the coordinate system, the velocity field in Eq. (10) can be expressed as

The pressure field can be presented by a quadratic function with respect to F (t) as follows:

The convective acceleration in Eq. (15) is presented as

Differentiating Eq. (35), the gradient of the pressure is presented as

Substituting Eqs. (36)–(38) into Eq. (15), the result is expressed as follows:

The zero-order perturbation equation is presented as the same in Eq. (19), and the second-order perturbation equation is presented as

Because dX_i = dx_i, integrating the above equation with respect to dx_i gives Bernoulli’s principle as follows:

The dynamic pressure field was determined using Eq. (41). Substituting Eq. (14) into Eq. (30) and integrating the result results in

The integral constant is determined as

6. Results and Verification

Some exact solutions to the NS equations exist for specific conditions (Antuono, 2020; Ethier and Steinman, 1994; Taylor and Green, 1937; Wang, 1991). For verification, the Taylor–Green vortex and Antonio’s solution were proven to be specific cases of the general solutions presented in this study.

6.5 Two-Dimensional Progressive Water Waves

Using the conventional coordinate system considered in Shin (2022, 2023), the phase is defined as β = kx −ωt where k is the wave number. Therefore, kG corresponds to ωt = k Ct for 2-D inviscid flow. Using Eqs. (42) and (44),

(55) kG(t)=ωμb(eμbt−1)

Using power series expansion, kG = ωt + O((μbt)²), where O stands for Big-O notation. Therefore, for inviscid flow, Eq. (55) gives ωt as defined in Shin (2022, 2023). “ μb” has a unit of angular frequency. Thus, a dimensionless quantity called the “damping ratio” is defined as ξdef=−μb/ω. Using the damping ratio, the phase for viscous flow is represented as

(56) β=kx−ωt{1−e−ξωtξωt}

Multiplying e^{− ξωt} by the velocity field in Shin (2022, 2023) gave a solution to the NS equation. Multiplying e^{− 2ξωt} to the dynamic pressure in Shin (2022, 2023) gives the pressure field for viscous flow. In the conventional coordinate system, the kinematic free surface boundary condition (KFSBC) is presented as

(57) v2=∂η∂t+v1∂η∂x1

where v₁ and v₂ are horizontal and vertical velocities, respectively, for viscous flow, and η is the wave profile. The KFSBC stipulates that there must be no flow crossing any fluid interface. The free surface, which is a fluid interface, remains unaffected by the viscosity. Consequently, the kinematic boundary condition should remain unaltered by viscosity. Using Eq. (10), the KFSBC is presented as

(58) e−ξωtu2=∂η∂t+e−ξωtu2∂η∂x1

where u₁ and u₂ are the horizontal and vertical velocities, respectively, for inviscid flow. The dimensionless elevation of the wave profile is defined as γ = k(η +h) in Shin (2022, 2023), where h is the water depth. Using the phase defined in Eq. (56),

(59) ∂η∂t=∂η∂β∂β∂t=∂β∂tddβ(γk−h)=−ωke−ξωtdγdβ

and

(60) ∂η∂x1=∂η∂β∂β∂x1=∂β∂x1ddβ(γk−h)=dγdβ

Substituting Eqs. (59) and (60) into Eq. (58), results in Shin (2022, 2023) because e^{− ξωt} ≠0. The KFSBC in Shin (2022, 2023) is still valid for viscous flow. Applying the dynamic boundary condition on the free surface to Eq. (35), the wave height for viscous flow is expressed as

(61) H=Hoe−2ξωt

where H_o is the wave height at t = 0. Waves emanating from a remote source are referred to as swells, whereas waves produced in the nearby vicinity are termed seas. Specifically, waves formed due to the wind-induced forced oscillation of water are generally identified as seas. In contrast, those originating from the natural oscillation of water post the cessation of wind are labeled as swells. Eqs. (56) and (61) indicate that the height of swell waves decreases when wind energy no longer influences wave formation, leading to a longer wave period and an increase in wavelength over time. The quantity “b” represents the strength of vorticity, σ defined in Shin (2022). b =− (kσ)² because the length is normalized by the wavelength in Shin (2022, 2023). Using ξω =− μb, the damping ratio is determined as follows:

(62) ξ=μ(kσ)2ω

A logarithmic decrease is used to determine the damping ratio. Using Eq. (61), the damping ratio can be expressed as

(63) ξ=14πln(H1H2)

where H₁, and H₂ represent the wave height of any two successive peaks. Although the strength of vorticity, σ, is determined using the average water particle velocities on the free surface and the seabed in Shin (2022, 2023), it can also be determined using Eqs. (62) and (63), which indicate that the currents considered in Shin (2022, 2023) are generated by the same source as the waves, but they differ from the Stokes drift. This study has proven that viscosity is the source of circulation and current in fluids.

Dean and Dalrymple (1984) provide the decay modulus for long waves with frictional damping as follows:

(64) α=πP1/4

where P is the Proudman number, Pdef=μ2/gk2h5. Shiau and Rumer (1974) carried out a series of experiments to examine the decay of seiches in a basin. The experiments were conducted in very shallow water (0.15 < h < 8.5 cm). The theoretical values of α in Eq. (64) were comparable to the experimental data in Fig. 1.

Fig. 1.

Decay modulus versus Proudman number [Experiment data from Shiau and Rumer (1974)].

The velocity fields for viscid flows, v_i, were checked. The velocity fields for inviscid flows, u_i, were compared with experimental data by Le Méhauté et al. (1968) in Shin (2022, 2023). Therefore, according to Eq. (10), the function, F(t) was checked instead of the velocity field, v_i, as follows. Eq. (61) is valid for long waves because Shin (2022, 2023) is valid even for ultra-shallow water waves. α = 2ξωT = 4πξ from Eq. (61) because the decay modulus is defined by the relative reduction in amplitude over one wave period. Using Eq. (62), the decay modulus is expressed as

(65) α=4π(kh)2σ2P

Note that ω² = gk²h for shallow water waves and μ=Pgk2h5. Using Eq. (11) in Shin (2023), Eq. (65) is expressed as

(66) α=4π{cos−1(1−u¯s1−u¯b)}2P

where u_s is the average speed in the horizontal direction over a period on still water and u_b is the average speed in the horizontal direction over a period on the seabed. ū_s = u_s/C and ū_b = u_b/C, where C is the celerity. Upon comparing Eq. (64) and Eq. (66), two distinct differences emerge between these equations. Eq. (64) is based on long-wave theory. Owing to the relatively small water depth in relation to the wavelength, variations in velocity in the vertical direction were disregarded, and only the average velocity was considered in Eq. (64). Consequently, the decay modulus in Eq. (64) is not dependent on the velocity distribution. On the other hand, the decay modulus in Eq. (66) is closely associated with the velocity distribution. Dean and Dalrymple (1965) used a quadratic friction law where the shear stress is directly proportional to the square of the velocity. The decay modulus in Eq. (64) was calculated assuming that bottom shear stress is presented by τ12=ρωμu1, where u₁ is the horizontal velocity. Therefore, the shear stress is directly proportional to μ. The present study used a Newtonian fluid model in which shear stress is linearly related to the strain rate, expressed as τ₁₂ = ρμ(u_1,2 + u_2,1), where x₁ is the horizontal direction and x₂ is the vertical direction. Therefore, the shear stress is directly proportional to μ . Consequently, this leads to a change in the exponent of P from 1/4 to 1/2. When the same constitutive equation and the average concept considered in the long wave theory are used, the average shear stress gradient is as follows.

(67) τ12,2=ρωμ{cos−1(1−u¯s1−u¯b)}u1h

Substituting Eq. (67) into the right side of Eq. (13), the decay modulus is presented as

(68) α=4π{cos−1(1−u¯s1−u¯b)}P1/4

For the velocity ratio, 1-U¯s1-U¯b=cos(14)=0.97, Eqs. (64) and (68) are equivalent to one another, as shown in Fig.1. Eqs. (64) and (68) are compared in Fig. 1. The comparison results indicated a significant level of concordance between Eq. (64) and Eq. (68). In contrast, this indicates a distinction between Eq. (64) and Eq. (66). The distinction elucidates the rationale behind Dean and Dalrymple’s (1965) choice to use a quadratic friction law instead of a Newtonian fluid model. Eqs. (67) and (68) express the applicability of this study even for non-Newtonian fluids, provided that the linearization of shear stresses is achieved using the technique described by Dean and Dalrymple (1984).

7. Conclusions

NS equations play a crucial role in elucidating various scientific and engineering phenomena. They are used in modeling ocean currents, water waves, and other water flows. They help with the design of the structure. Despite the wide range of practical applications, the equations are mainly solved using numerical or experimental methods, while there are some exact solutions to specific problems.

General solutions to NS equations for incompressible flow are presented in this study. The solutions were calculated by combining the solutions of the Helmholtz equation and the solution of a single-degree-of-freedom damped harmonic oscillator. They are analytical exact solutions to NS equations.

The velocity fields and the vorticity fields were presented with solutions of the Helmholtz equations. When the velocity fields are expressed as the product of the solution of the Helmholtz equation and a univariate function of time, the NS momentum equation is decomposed into the Euler momentum equations and an ordinary differential equation for the univariate function of time. The Euler momentum equations were converted to Bernoulli’s principle to determine the pressure field. The time function was determined by solving the ordinary differential equation. Consequently, the general solutions were derived through the linear combinations of the solutions to the Helmholtz equation. The solutions presented in this study are functions of class C^∞ because solutions of the Helmholtz equations are smooth (Lin, 1988).

The Taylor–Green vortex solution and the solutions proposed by Antuono (2020) are particular cases of these solutions, and these solutions satisfy the vorticity equations. The solutions presented by Shin (2022, 2023) were enhanced to integrate viscid water waves with shear currents.

As time elapses or as a swell moves away from its generating area, its height decreases while the period and wavelength increase. Eventually, the swell dissipates entirely because of the removal of the energy source. This phenomenon was substantiated mathematically in the present study. The decay modulus was compared with the experimental data from Shiau and Rumer (1974) and the solution proposed by Dean and Dalrymple (1984). The agreement was excellent. Unlike Dean and Dalrymple (1984), the study showed that the decay modulus is influenced by the velocity gradient in the depth direction.

The technique presented in this study can also be applied to non-Newtonian fluid, provided that the linearization of shear stresses is achieved through the technique described by Dean and Dalrymple (1984).

Article information Continued