Free Surface Boundary Conditions Research Articles

Finite-difference method for modeling the surface wave propagation with surface topography in anisotropic-viscoelastic media

The growth of urbanization makes it increasingly important to accurately investigate underground spaces. Surface-wave-based inversion methods are becoming popular due to their cost-effectiveness and efficiency in reconstructing underground space. However, these methods may produce inaccurate results in the presence of media anisotropy or/and viscoelasticity. To investigate the influence of media anisotropy or viscoelasticity on surface waves, it is crucial to develop an efficient and accurate forward modeling method in anisotropic-viscoelastic (AV) media. The finite-difference (FD) method is a widely used forward modeling method in surface-wave-based inversion methods. However, implementing the free-surface boundary condition and representing non-flat topography can be challenging. Considering the challenges described above, this study proposes a simple and efficient FD method for modeling surface-wave propagation. The proposed method employs a parameter-modified strategy to implement the stress-free condition by modifying the model parameters near the (non-) flat free-surface boundary. A staircase discretization strategy is then used to represent the non-flat free surface. To suppress the staircase diffractions caused by the staircase discretization strategy, an independent wavefield superposition is adopted with modeling results of different parameter configurations. Compared with conventional FD methods, the proposed method can reduce computational costs while accurately simulating surface waves. Numerical experiments have demonstrated the accuracy and feasibility of the proposed method. The proposed method provides a theoretical basis for surface-wave-based inversion methods and contributes to the development of accurate underground space investigation.

Effect of free surface boundary condition on reflected waves in a half-space piezoelectric semiconductor medium

Effect of free surface boundary condition on reflected waves in a half-space piezoelectric semiconductor medium

3D Sensitivity Kernels With Full Attenuation Computed by a Combination of the Strong Stability Preserving Runge‐Kutta Method and the Scattering‐Integral Method

AbstractTo date developments of seismic attenuation models of the Earth have lagged those of velocity models. This is partly due to difficulties in isolating waveform perturbations caused by attenuation and velocity heterogeneities and partly due to different theories (e.g., ray theory, finite frequency theory) used in most previous and current studies to approximate sensitivity kernels and invert for the attenuation structure. We present in this paper a new method for computing the 3D waveform Fréchet kernels that account for full physical‐dispersion and dissipation attenuation. For solving the 3D isotropic anelastic wave equation described by the generalized Maxwell body model, we extend our previously proposed full waveform modeling method in Cartesian coordinates to that in spherical coordinates, which can provide stable numerical solutions even in the presence of strong attenuation. Then, we apply the scattering‐integral method for calculating 3D travel time and amplitude sensitivity kernels with respect to velocity and attenuation structures. We demonstrate the accuracy of our forward method and the effectiveness of the implementation of absorbing and free surface boundary conditions through numerical tests. Moreover, by choosing the Northwestern United States region as a realistic example, we verify the accuracy of the computed 3D sensitivity kernels through comparing the waveform measurements with predictions from the kernels. Finally, we discuss the importance of calculating full anelastic sensitivity kernels including both effects of physical dispersion and dissipation, where we specially explore the effect of scattering due to random velocity and attenuation heterogeneities on waveform measurements.

PolyStokes: A Polynomial Model Reduction Method for Viscous Fluid Simulation

Standard liquid simulators apply operator splitting to independently solve for pressure and viscous stresses, a decoupling that induces incorrect free surface boundary conditions. Such methods are unable to simulate fluid phenomena reliant on the balance of pressure and viscous stresses, such as the liquid rope coil instability exhibited by honey. By contrast, unsteady Stokes solvers retain coupling between pressure and viscosity, thus resolving these phenomena, albeit using a much larger and thus more computationally expensive linear system compared to the decoupled approach. To accelerate solving the unsteady Stokes problem, we propose a reduced fluid model wherein interior regions are represented with incompressible polynomial vector fields. Sets of standard grid cells are consolidated into super-cells, each of which are modelled using a quadratic field of 26 degrees of freedom. We demonstrate that the reduced field must necessarily be at least quadratic, with the affine model being unable to correctly capture viscous forces. We reproduce the liquid rope coiling instability, as well as other simulated examples, to show that our reduced model is able to reproduce the same fluid phenomena at a smaller computational cost. Futhermore, we performed a crowdsourced user survey to verify that our method produces imperceptible differences compared to the full unsteady Stokes method.

An overset mesh-free finite-difference method for seismic modeling including surface topography

An overset mesh-free finite-difference method (FDM) is developed for 2D time-domain seismic wave modeling using second-order elastic displacement wave equations in a generally isotropic medium that exhibits complex surface topography and subsurface structures. The complex surface topography is discretized using an overset mesh-free node generation method, thus naturally arranging surface nodes on the surface topography and outperforming the Cartesian grid in terms of geometric flexibility. A Cartesian grid is used to discretize the entire computational region in terms of computational efficiency. The QR-decomposition radial-basis-function FDM (QR-RBF-FDM) and FDM are adopted to solve the equation independently in the mesh-free node and Cartesian grid. The data transfers between them are conducted through QR-RBF interpolation to facilitate mesh-free calculation, thus eliminating the requirement for node-to-grid matching. The improved free-surface boundary condition is implemented in the QR-RBF-FDM. Furthermore, a more appropriate formula for the outward normal direction is developed. The numerical solutions of our method are compared with those of the curvilinear grid FDM and standard FDM to verify the effectiveness of our method. The results indicate that our method avoids numerical scattering induced by staircase approximations of Cartesian grids while saving considerable computational costs. Therefore, this method is a better candidate than the mesh-free FDM for seismic numerical simulation with free surface topography.

Numerical modeling on global-scale mantle water cycle and its impact on the sea-level change

Using numerical modeling of the global-scale mantle water cycle, I have clarified the mechanisms governing the deep Earth water content variations and evaluated their effects on sea-level changes. The following points are found: 1. The temporal variations of water content in the deep Earth are mainly determined by a balance between the ingassing flux from the ocean and the dehydration flux in the mantle wedge. 2. The primary control of ingassing flux comes from the water content of oceanic crust, rather than the style of plate subduction promoted by the free surface boundary condition. 3. The water content of oceanic crust ranges from 0.1 to 1 wt. percent to match the ingassing flux with the inferred one from the geological analysis.Based on the best-fit scenario for the ingassing flux from the ocean to the deep interior, I obtain an average rate of sea-level change of −1 to −1.5 m/Myrs over 1 billion year, which is nearly consistent with the long-term trend of sea-level change estimated from geological and paleontological data (∼−1 m/Myr). Note, however, that the detailed profile of sea-level change in the numerical modeling is difficult to explain the observation one. This underlines the importance of dehydration process in the mantle wedge in understanding the sea-level change resulting from the water budget in the Earth's deep interior over the Earth's history.

A new elastic instability in gravity-driven viscoelastic film flow

We examine the linear stability of the gravity-driven flow of a viscoelastic fluid film down an inclined plane. The viscoelastic fluid is modeled using the Oldroyd-B constitutive equation and, therefore, exhibits a constant shear viscosity and a positive first normal stress difference in viscometric shearing flows; the latter class of flows includes the aforesaid film-flow configuration. We show that the film-flow configuration is susceptible to two distinct purely elastic instabilities in the inertialess limit. The first instability owes its origin entirely to the existence of a free surface and has been examined earlier [Shaqfeh et al., “The stability of gravity driven viscoelastic film-flow at low to moderate Reynolds number,” J. Non-Newtonian Fluid Mech. 31, 87–113 (1989)]. The second one is the analog of the centermode instability recently discovered in plane Poiseuille flow [Khalid et al., “Continuous pathway between the elasto-inertial and elastic turbulent states in viscoelastic channel flow,” Phys. Rev. Lett. 127, 134502 (2021)] and owes its origin to the base-state shear; it is an example of a purely elastic instability of shearing flows with rectilinear streamlines. One may draw an analogy of the aforesaid pair of unstable elastic modes with the inertial free-surface and shear-driven instabilities known for the analogous flow configuration of a Newtonian fluid. While surface tension has the expected stabilizing effect on the Newtonian and elastic free-surface modes, its effect on the corresponding shear modes is, surprisingly, more complicated. For both the Newtonian shear mode and the elastic centermode, surface tension plays a dual role, with there being parameter regimes where it acts as a stabilizing and destabilizing influence. While the Newtonian shear mode remains unstable in the limit of vanishing surface tension, the elastic centermode becomes unstable only when the appropriate non-dimensional surface tension parameter exceeds a threshold. In the limit of surface tension being infinitely dominant, the free-surface boundary conditions for the film-flow configuration reduce to the centerline symmetry conditions satisfied by the elastic centermode in plane Poiseuille flow. As a result, the regime of instability of the film-flow centermode becomes identical to that of the original channel-flow centermode. At intermediate values of the surface tension parameter, however, there exist regimes where the film-flow centermode is unstable even when its channel-flow counterpart is stable. We end with a discussion of the added role of inertia on the aforementioned elastic instabilities.

Frequency-domain elastic wave modeling for vertical transversely isotropic media with an average-derivative optimal method

Based on an average-derivative method (ADM), it is feasible to develop an optimal elastic wave finite-difference nine-point scheme for vertical transversely isotropic (VTI) media in the frequency domain. The optimization coefficients of ADM can be determined by minimizing phase-velocity dispersion error. Nevertheless, for the complex VTI media, the point-by-point optimization, which only considers a single anisotropic parameter pair, will require high computational cost. To address this issue, we take into account the other two low-cost schemes: the overall simultaneous optimization in terms of various anisotropy parameter pairs and the interpolation based on the completely computed optimization coefficients. Applying the 58 measured anisotropy parameter pairs of sedimentary rocks, we test the performances of point-by-point optimization, the overall optimization, interpolation, and classical nine-point method (CM). The results indicate that the number of grid points per qS wavelength of the interpolation scheme is closest to that of the point-by-point optimization scheme and much lower than that of the overall optimization scheme and CM. In addition, the surface of solid earth is a strong impedance-discontinuity surface (free surface), resulting in the more complicated wavefields. To realize the full-waveform simulation for elastic VTI media, it is necessary to implement the free-surface boundary condition in detail. Based on the equivalence between the stress-image and adaptive parameter-modified schemes for implementing the free-surface boundary condition in the time domain, we give a new frequency-domain free-surface expression for the VTI media. We still call this expression the stress-image scheme. Numerical experiments demonstrate that (1) the optimal ADM can produce more accurate results than CM at almost the same computational cost, (2) the interpolation scheme has almost identical accuracy with the point-by-point optimization scheme, and (3) the stress-image scheme is feasible for the implementation of free-surface boundary condition in a VTI medium.

Finite-difference modeling with topography using 3D viscoelastic parameter-modified free-surface condition

The free-surface boundary condition is a crucial aspect in the numerical modeling of (visco)elastic wave equations, especially when using a finite-difference (FD) method in the presence of surface topography. The parameter-modified method is a widely used approach to solve this problem. In this regard, the adaptive Poisson’s ratio parameter-modified method has proven to be effective in accurately simulating seismic surface waves within the FD discretization framework. Based on the equivalent medium theory, vacuum approximation, and mathematical limit, we develop a viscoelastic parameter-modified (VPM) method for the implementation of the free-surface boundary condition in the 3D viscoelastic wave equation. Our approach modifies the viscoelastic constitutive relation and density at the free surface and provides a formulation in terms of displacement and stress. We determine that our VPM method is more general than the viscoelastic stress-image method as it includes the latter as a limit case when the Poisson’s ratio equals zero. The presented free-surface method is represented implicitly within the FD grid, and we provide implementation details when simulating surface waves with topography in the standard staggered-grid FD scheme. We support our theory’s feasibility and accuracy through numerical examples.

Numerical study on mean wave drift forces and moment of ships using Rankine panel method in frequency domain

A three-dimensional frequency-domain seakeeping model based on the Rankine panel method is developed to compute the mean wave drift forces and moment acting on ships with forward speeds. Two different linearization models (the double-body and Neumann–Kelvin models) are applied to evaluate the mean wave drift forces and moment using a near-field method. In these models, the Rayleigh’s artificial friction is adopted in the free-surface boundary condition to suppress longer wave components in computational regions apart from the ship. A novel upwind difference scheme is introduced to calculate the first-order and second-order derivatives of the velocity potential on the free surface. Wigley I, Wigley III and S175 containership are adopted in numerical simulations to validate the effectiveness of the developed model. In the present study, the motion response amplitude operators in head, beam and oblique seas obtained using the two models are compared, and their contributions to the mean wave drift forces and moment are analyzed. Results show that in head seas, the pitch motions have significant effects on the surge-directional wave drift force Fx. In beam seas, the heave motions have the most effect on the mean wave drift forces and moment (Fx, Fy and Mz).

Pressure analysis in nonlinear waves by revisiting the breaking wave onset

Breaking waves are numerically simulated by solving the two dimensional nonlinear free surface boundary conditions in potential theory. Onset of breaking wave is determined for two configurations: A focused wave in intermediate depth and a dambreaking flow leading to a soliton. Two criteria are identified: first, a kinematic criterion follows from the ratio of the horizontal velocity at the crest to the phase velocity; second, a dynamic criterion is formulated in terms of the vertical pressure gradient. Thresholds of these criteria are proposed in the literature. For the two configurations studied here, it is shown that these thresholds are not universal but they are reasonable estimates of the onset of breaking wave. To facilitate this study, the Hessian matrix of the pressure is analyzed. It is shown that there exists a continuous curve that links the crest of the wave to the sea bottom. This line follows from the colinearity of the pressure gradient and one of the two eigenvectors of the Hessian matrix. The analysis of the pressure, its gradient and its second derivatives along this line gives more insights into the definition of the thresholds associated to the onset of breaking wave.

A unified steady and unsteady formulation for hydrodynamic potential flow simulations with fully nonlinear free surface boundary conditions

This work discusses the correct modeling of the fully nonlinear free surface boundary conditions to be prescribed in water waves flow simulations based on potential flow theory. The main goal of such a discussion is that of identifying a mathematical formulation and a numerical treatment that can be used both to carry out transient simulations, and to compute steady solutions — for any flow admitting them. In the literature on numerical towing tank in fact, steady and unsteady fully nonlinear potential flow solvers are characterized by different mathematical formulations. The kinematic and dynamic fully nonlinear free surface boundary conditions are discussed, and in particular it is proven that the kinematic free surface boundary condition, written in semi-Lagrangian form, can be manipulated to derive an alternative non penetration boundary condition by all means identical to the one used on the surface of floating bodies or on the basin bottom. The simplified mathematical problem obtained is discretized over space and time via Boundary Element Method (BEM) and Implicit Backward Difference Formula (BDF) scheme, respectively. The results confirm that the solver implemented is able to solve steady potential flow problems just by eliminating null time derivatives in the unsteady formulation. Numerical results obtained confirm that the solver implemented is able to accurately reproduce results of classical steady flow solvers available in the literature.

Modeling evaporation with a meshfree collocation approach

In this paper, a new model for the below-boiling point evaporation process with a meshfree collocation method is developed. In order to capture the phase change process, two different approaches are proposed: multi-phase and single-phase. First, a multi-phase approach is considered, where a novel mass transfer model assumes that the diffusion driven by the vapor concentration gradient in the air phase near the interface is the primary driving force for the mass transfer between phases as both the liquid water and air/vapor phases are simulated. Then, a water-only single-phase approach is also proposed, in which only the liquid water phase is simulated. For this, appropriate free surface boundary conditions are developed based on the convective mass transfer theory to model evaporation and incorporate airflow effects without explicitly simulating the air phase. In order to validate the proposed models, a series of experiments with varying air temperature, relative humidity, and airflow rate is conducted. The numerical results show a good agreement with the evaporation rate measured in the experiments. The multi-phase simulations agree better with the experiments, while the single-phase simulations also produce good results with a much lower computational effort.

Forces of fully nonlinear interfacial periodic waves on a cylindrical pile in a two-layer fluid with free-surface boundary conditions

In the frame of fully nonlinear potential flow theory, series solutions of interfacial periodic gravity waves in a two-layer fluid with free surface are obtained by the homotopy analysis method (HAM), and the related wave forces on a vertical cylinder are analyzed. The solution procedure of the HAM for the interfacial wave model with rigid upper surface is further developed to consider the free surface boundary. And forces of nonlinear interfacial periodic waves are estimated by both the classical and modified Morison equations. It is found that the estimated wave forces by the classical Morison equation are more conservative than those by the modified Morison’s formula, and the relative error between the total inertial forces calculated by these two kinds of Morison’s formulae remains over 25% for most cases unless the upper and lower layer depths are both large enough. It demonstrates that the convective acceleration neglected in the classical Morison equation is rather important for inertial force exerted by not only internal solitary waves but also interfacial periodic waves. All of these should further deepen our understanding of internal periodic wave forces on a vertical marine riser.

Fractional-step finite difference schemes for incompressible elasticity on overset grids

Efficient finite-difference schemes for the numerical solution of the time-dependent equations of incompressible linear elasticity on complex geometry are presented. The schemes are second-order accurate in space and time, and solve the equations in displacement-pressure form. The algorithms use a fractional-step approach in which the time-step of the displacements is performed separately from the solution of a Poisson problem to update the pressure. Complex geometry is treated with curvilinear overset grids. Compatibility boundary conditions and an upwind dissipation for wave equations in second-order form are included to ensure stability for overset grids, and for the difficult case of traction (free-surface) boundary conditions. A divergence damping term is added to keep the dilatation small. The stability of the schemes is studied with GKS mode analysis. Exact eigen-mode solutions are obtained for several problems involving rectangular, cylindrical and spherical configurations, and these are valuable as benchmark problems. The accuracy and stability of the approach in two and three space dimensions is illustrated by comparing the numerical solutions with the exact solutions of the benchmark problems.

A regularity result for the incompressible inviscid magnetohydrodynamics equations in the Arbitrary Lagrangian-Eulerian coordinates

Consider herein the incompressible inviscid magnetohydrodynamics equations in a moving domain with free surface boundary conditions. Without smallness assumption in the fluid volume, a priori estimates for solutions of this model are established when the initial data in H2.5+δ (δ>0) upon the Arbitrary Lagrangian-Eulerian coordinates under the boundary condition p+12|B|2=0. Indeed, this is achieved by reformulating the system into a new formulation with the Arbitrary Lagrangian-Eulerian variables, presenting the uniform estimates for the pressure, the tangential estimates for the system, as well as the curl and divergence estimates.

The Effect of a Linear Free Surface Boundary Condition on the Steady-State Wave-Making of Shallowly Submerged Underwater Vehicles

Near-surface simulation methods for shallowly submerged underwater vehicles are necessary for the population of a variety of free-surface-affected, coefficient-based maneuvering and seakeeping models. Simulations vary in complexity and computational costs, often sacrificing accuracy for simplicity and speed. One particular simplifying assumption, the linearization of the free surface boundary conditions, is explored in this study by comparing the steady-state wave-making characteristics of a shallowly submerged prolate spheroid using two different simulation methods at several submergence depths and forward speeds. Hydrodynamic responses are compared between a time-domain boundary element method that makes use of a linearized free surface boundary condition and an inviscid, volume of fluid Reynolds-Averaged Navier–Stokes computational fluid dynamics code that imposes no explicit free surface boundary condition. Differences of up to 22.6%, 32.5%, and 33.3% are found in the prediction of steady state surge force, heave force, and pitch moment, respectively. The largest differences between the two simulation methods arise for motions occurring at small submergences and large wave-making velocities where linear free-surface assumptions become less valid. Nonlinearities that occur in such cases are revealed through physical artifacts such as wave steepening, wave breaking, and high-energy waves. A further examination of near-surface viscous forces reveals that the viscous drag on the vessel is depth dependent due to the changing velocity profile around the body.

A Simplified Numerical Method for Simulating the Generation of Linear Waves by a Moving Bottom

In this study, simplified linear numerical method that can simulate wave generation and transformation by a moving bottom is introduced. Numerical analysis is conducted in wave number domain after continuity equation, linear dynamic and kinematic free surface boundary conditions and linear kinematic bottom boundary condition are Fourier transformed, and the results are expressed in space domain by an inverse Fourier transform. In the wavenumber domain, the dynamic free water surface boundary condition and the kinematic free water surface boundary condition are numerically calculated, and the velocity potential in the mean water level (<i>z</i> = 0) satisfies the continuity equation and the kinematic bottom boundary condition. Wave generation and transformation are investigated when the triangular and rectangular shape of bottoms move periodically. The results of the simplified numerical method are compared with the results of previous analytical solutions and agree well with them. Stability of numerical results according to the calculation time interval (△<i>t</i>) and the calculation wave number interval (△<i>k</i>) was also investigated. It was found that the numerical results were appropriate when △<i>t</i>≤<i>T</i>(period)/1000 and △<i>k</i>≤π/100.

Generalized finite difference method applied to solve seismic wave propagation problems. Examples of 3D simulations

The simulation of seismic wave propagation generally requires dealing with complex tridimensional geometries that are irregular in shape and have non‐uniform properties, features that make the application of the generalized finite difference method in this field interesting. This work continues the extensive developments by the research team focused on the simulation of seismic wave propagation in two‐dimensional domains. In this new contribution, the general formulation and the treatment of free surface boundary conditions are extended for the three‐dimensional case, and the results obtained from different examples are analyzed.

Elastic full-wave field simulation in 3D tunnel space with the variable staggered-grid finite-difference method in cylindrical coordinates

Tunnel advance detection technology is an important method for determining the structure of complex geological bodies in front of the tunnel face. Among the tunnel advance detection technologies, the seismic method is one of the most accurate methods with long detection distances. In seismic tunnel advance detection, the cylindrical configuration aggravates the complexity of the wave field in the 3D tunnel space and significantly influences the accuracy of the detection results. Thus, it is crucial to simulate an accurate seismic full-wave field of the tunnel space and to understand the propagation and wave-field characteristics of individual seismic waves for seismic tunnel advance detection. Usually, in 3D Cartesian coordinates, the tunnel wall is approximated with a staircase boundary, but it is not sufficiently accurate in shape and generates numerical dispersion in the simulation, especially in the presence of surface waves. Therefore, we developed a variable staggered-grid finite-difference method in cylindrical coordinates to simulate the elastic full-wave field in a 3D tunnel space. Properly setting free-surface boundary conditions solves the propagation of surface waves on the tunnel wall and face. The interference of the instability and discontinuity of the pole axis in the seismic wave field simulation was eliminated using our method. According to practical engineering situations, models containing three types of geological bodies in front of the tunnel face were designed. Their seismic wave propagation and wave field characteristics were analyzed. Compared with the simulation results in Cartesian coordinates, the results in cylindrical coordinates show that numerical dispersion is negligible and conclude that a higher signal-to-noise ratio and more accurate seismic wave field can be simulated with cylindrical coordinates in the 3D tunnel space. Our simulation method provides theoretical and practical guidance for analyzing and interpreting seismic wave fields in tunnel advance detection.