Introduction. Issues regarding the dynamic analysis of the behavior of structures and buildings featuring complex foundations on multilayer bases, as well as the analysis of multilayer systems with local inclusions and inhomogeneities [1], necessitate the further development of dynamic analysis methods and the creation of efficient algorithms for calculating wave fields in such media [2]. Widely used finite element methods encounter a number of difficulties when applied to semi-infinite layered media, primarily due to the presence of reflecting boundaries and the requirement to satisfy wave radiation conditions at infinity. Enhancing the efficiency of numerical calculation methods for systems with a large number of layers, specifically through the use of various parallelization schemes, is a topic of significant interest. This paper advances an approach that combines finite element methods (FEM) and boundary element methods (BEM) for the analysis of multilayer media in contact with complex structures.

Hybrid methods combining the finite element method and the boundary element method have been extensively studied in the scientific literature for solving problems in acoustics, fluid-structure interaction, and wave scattering. In [3], an iterative approach is presented for analyzing the interaction between saturated porous media and sealed elastic media using a combined FEM-BEM technique. This method accounts for nonlinearities within the FEM subdomains. The authors [4] develop an explicit direct FEM-BEM coupling procedure for nonlinear dynamics, allowing subdomains to be analyzed independently with different time steps, thereby enhancing the flexibility and efficiency of coupled analyses. In [5], an FEM-BEM coupling method for time-domain fluid-structure interaction analysis is introduced. This method employs the Poincaré-Steklov operator and enables the derivation of a priori error estimates for model problems. Paper [6] analyzes an overlapping continuous model for the Helmholtz wave propagation problem in unbounded domains. The method accurately preserves the radiation condition and demonstrates efficiency for heterogeneous media. Paper [7] develops a FEM-BEM coupling strategy for time-harmonic acoustic scattering in media with variable sound speed. The method employs an interface variable at the coupling boundary and demonstrates quasi-optimal convergence. The author [8] proposes a hybrid numerical method for solving phase-change problems with moving boundaries. The method combines radial basis functions with a grid-based approach, enabling the solution to solidification and dissolution problems. Paper [9] investigates a posteriori error estimates for the virtual element method (VEM) and FEM-BEM-based methods. The method allows for calculating the distance between the finite element gradient and its post-processed version, as well as identifying points of superconvergence. The author [10] proposes a regularizer based on On-Surface Radiation Conditions (OSRC) to stabilize the coupling between the FEM and the BEM. This method improves the convergence of iterative solvers, specifically at high frequencies, and reduces computational costs. Paper [11] develops a deep learning framework for predicting transfer functions in structural-acoustic models, enabling a significant reduction in computation time compared to traditional FEM-based calculations. The authors [12] propose a framework for acoustic analysis and shape optimization of 3D doubly periodic multilayer structures. The method employs the isogeometric boundary element method (IG-BEM) to optimize the shape of acoustic structures. Paper [13] presents a semi-analytical finite element method for modeling 2D guided waves in composite laminates, providing the analysis of wave scattering at delamination and the determination of reflection and transmission coefficients. In [14], high-order and plane-wave basis functions for solving Helmholtz problems are compared. It demonstrates that the NURBS basis outperforms the Lagrange basis and analyzes the errors and numerical dispersion associated with various bases. Paper [15] proposes a hybrid FEM-BEM method for solving acoustic-structural problems. The method employs a discrete shear gap to construct a locking-free FEM formulation and utilizes BEM collocation based on the Burton-Miller formulation. In [16], a symmetric FEM-BEM coupling scheme is developed for the scattering of transient acoustic waves by inhomogeneous anisotropic obstacles. The method enables the discretization of transmitted and scattered waves and demonstrates stability in Hilbert space.

A significant number of research studies are focused on developing numerical methods for analyzing layered, anisotropic, and functionally graded media. Paper [17] proposes a modified scaled boundary finite element method (SBFEM) for analyzing the dynamic response of anisotropic layered media. This method provides the solution to problems involving cross-anisotropic and generally anisotropic media, while also simplifying calculations for cross-anisotropic layered media. In [18], the SBFEM is applied to determine the dynamic response of an infinite three-dimensional soil medium. A continuous subdivision algorithm allows for the derivation of the dynamic stiffness matrix of layered soil in the frequency domain and facilitates the analysis of embedded foundations. In [19], a 3D approach is developed that couples FEM and BEM in the frequency domain to analyze the dynamic response of structures on layered transversally isotropic half-spaces. The method uses Green's functions and allows taking into account the singular behavior of stresses at interfaces. In [20], a new methodology is proposed for the analysis of longitudinally invariant problems of soil-structure interaction in elastodynamics. The method uses 2.5D FEM-BEM coupling and improves computational efficiency by using the fundamental solution method to model wave propagation. In [21], a method is developed for calculating 2.5D elastodynamic Green's functions for layered half-spaces in Cartesian coordinates. The method provides obtaining accurate solutions for homogeneous and layered media. Paper [22] proposes a 2.5D FEM-BEM methodology to analyze the dynamic response of segmented tunnels in saturated soil. The method allows taking into account the effect of tunnel connections on dynamic stresses and pore water pressure, which is underestimated in homogeneous models. In [23], the boundary-based finite element method (BFEM) for stress analysis in anisotropic functionally graded materials is presented. The method uses piecewise homogeneous layers and takes into account the effects of anisotropy and gradient without constructing a mesh of the area.

Numerous modern research studies focus on the optimization and stress analysis of heterogeneous, composite, and functionally graded materials. A stabilized iRBF method for analyzing elastoplastic problems is proposed in [24]. This method minimizes computational costs while providing high accuracy. An approach for the topological optimization of damping layers subjected to harmonic excitation is developed in [25]. The method employs a piecewise-constant level-set function and sensitivity analysis.

While research in this field enables the analysis of a wide range of complex systems, the dynamics of a multilayered foundation system containing a partially embedded object remain insufficiently studied. The objectives of this work are to analyze the stress-strain state of a complex structure resting on a multilayered foundation, utilizing a coupled FEM-BEM approach, and to conduct a comparative analysis of the efficiency of various parallel computing technologies for optimizing the calculation of these characteristics. Achieving these goals requires examining a modeling methodology for multilayered foundations that combines these methods, as well as addressing the modeling of the complex structure using a custom-developed software tool that implements this methodology and employs parallel computing technologies based on Fortran.

Materials and Methods. The Lamé equation provides a fundamental description of the motion of an isotropic linear-elastic medium in terms of displacements linking the Lamé elastic constants to the displacement vector field and body forces. Applying the Fourier transform to this equation enables a transition from a differential formulation in the spatial domain to algebraic relations in the frequency domain, where fundamental solutions take an explicit form. This serves as the basis for constructing the boundary element method, in which integral representations involving the Fourier transform are used to formulate systems of boundary integral equations that reduce the three-dimensional elasticity problem to a two-dimensional discretization of the boundary.

To illustrate this transition, let us consider the behavior of some object located on the surface or partially embedded in a multilayer half-space in the mode of steady-state harmonic oscillations with a frequency. As an example, consider a construction site with a pile foundation. The solution to the problem is constructed separately: for part of the layered soil and the object with a foundation using the finite element method, which allows taking into account all the design features of the object and foundation, and the remaining part of the foundation based on the method of boundary integral equations (boundary elements) (Fig. 1).

The behavior of all components of the structure and the N-layer base is considered linear elastic. In this case, the properties of the medium in its j-th component are described by density ρj and Lamé coefficients λj, μj (j = 0, 1, …, N – 1). The contact of all media is assumed to be rigid, determining the equality of the normal and tangential elements of the stress and displacement vectors at the interface boundaries.

The motion of the medium is determined through solving the system of Lamé equations regarding the amplitude functions of the displacements of points of the medium [26]:

(1)

Here, the coordinate system (x, y, z) is related to the surface of the multilayer half-space, as shown in Figure 1.

Without limiting the generality of the problem statement, we believe that oscillations are excited by forces distributed in limited area Ω on the flat boundary of the half-space at x = 0:

At this, linear parameter α in formula (1) corresponds to the characteristic size of this region. The problem formulation is completed by the requirement that the solutions satisfy the limiting absorption principle [27].

To describe the dynamic behavior of the lower part, represented by a multilayer half-space with a cavity, given a formally known stress vector τ on its surface S and specified surface stresses Rj, we employ the boundary integral equation method based on the dynamic reciprocity theorem and a specific construction of the fundamental solution matrices.

The form of the fundamental solution matrices is based on the solution to an auxiliary problem concerning the vibrations of a homogeneous half-space (x > 0) characterized by the elastic parameters of the j-th layer subject to a prescribed stress vector X1(y, z) on its surface x = 0 in a local coordinate system.

To this end, we express the displacement field in a matrix form subjected to a Fourier transform:

(2)

A stress vector can be represented in a similar way:

Vector C(α, β) is determined from the boundary condition on the surface x = 0 in the form:

The resulting analytical matrix representations facilitate the solution to the problem for a layer and a multilayer half-space, as well as the efficient application of computational parallelization methods.

To construct the fundamental solution matrices, we assume that a concentrated force acts at point r0 = (x0, y0, z0) of the j-th elastic layer 0 < x < hj:

Then the displacements of a point with coordinates (x, y, z), satisfying equation (1) with the addition of a term for the body force qδ(r – r0) on the left-hand side are given by an integral of the form:

(3)

By analogy, the fundamental solutions for the stress state of a layer with a source in the form of a double Fourier transform are determined:

(4)

To construct fundamental solution matrices in the Fourier-transformed form, the principle of superposition of the solutions to three problems is used:

(5)

where , corresponding to formula (2), represents the field of waves radiated from and reflected by the layer boundary; x = 0; — field of waves radiated and reflected from the layer boundary x = hj; — field of an oscillation source in unbounded space; easily obtained from by switching to a new coordinate system: x → h – x, change of signs of the unit vectors ex, ez and the introduction of a new stress vector instead of .

The matrix components can be obtained by applying the re-expansion formulas for spherical waves in terms of a set of plane waves, expressed as a double Fourier integral [28], which results in:

,

,

,

,

,

,

,

.

The stress components admit a similar representation. If the source of oscillations is absent from the selected j-th layer, the third term in representation (5) is omitted.

To definitively determine the fundamental solutions, it is required to eliminate the vector functions for each component of the multilayered half-space. This is implemented by applying continuity conditions for the displacement and stress vectors u(j+1) = u(j), t(j+1) = t(j) expressed in the form of Fourier transforms at the interfaces between the x = xj, j = 0, 1 ,…, N – 2 of the j-th and j + 1-th components of the layered medium, as well as the free-boundary conditions at x = 0. The resulting system of linear algebraic equations and its determinant characterize the dispersion properties of the multilayer half-space.

Using the dynamic reciprocity theorem, we subsequently derive boundary integral equations on the surface S for the nodal values of the displacement and stress vectors, which take the following form [29]:

(6)

Here, up, τk — components of the stress and displacement vectors, respectively, — components of the constructed matrices of fundamental displacements and stresses at point from the action of the source at point r0, Rk — components of the external force vector on the boundary of the half-space, nk — components of the outer normal vector to S.

By discretizing system (6) using the linear boundary element method, a system of linear algebraic equations SLAE is obtained for the unknown nodal displacements uS and forces FS:

(7)

The problem under consideration reduces to the simultaneous solution to the BEM SLAE and the FEM [30]:

(8)

Here, M — mass matrix, C — resistance matrix, K — structural stiffness matrix, u — amplitude functions of the nodal displacement vector, — nodal load vector, ω — frequency of steady-state oscillations. Note that vectors u and F contain variables for nodes lying on the common boundary S, which are unknown quantities.

After jointly solving the SLAE (7, 8), the resulting solution is substituted into the known representations of wave fields (2, 3, 4), which makes it possible to calculate the dynamic characteristics of the stress-strain state of a multilayer foundation and consider its interaction with the associated partially embedded object.

Consider the problem of determining normal stresses along the axis X at a specific point within a three-layer spatial foundation subjected to dynamic surface loading in the 15–25 Hz frequency range. The bottom layer of the foundation is an elastic half-space, while a structure comprising a thin-walled steel cube with a partially embedded concrete foundation is situated on the surface (Fig. 2).

Several problem formulations are considered — with and without a building structure, and with varying stiffness configurations of the foundation media. The physical properties of the foundation layers and structural elements are presented in Tables 1 and 2, respectively, where ρ — density, E — Young's modulus, ν — Poisson's ratio.

Load q = 1 kPa is applied to the surface at point (x0; y0; z0) = (0; 4.7; 0), one of the bottom corners of the cube is located at point (x1; y1; z1) = (0; –1; –1), and the observation point coordinates are (x; y; z) = (0; 1; 0).

To solve the problem, a combination of a custom software product, implementing the authors' BEM for determining the stress-strain state of a multilayered foundation, and the Ansys Mechanical software suite, using the FEM to analyze the structure and surrounding soil, is used. The numerical experiment was conducted on a PC equipped with an Intel Core i5-13600K processor and 32 GB of RAM (6000 MT/s).

Research Results. Data on the amplitude-frequency characteristics of the observation point are obtained (Fig. 3).

For cases a and b (stiff surface layers and layers of equal stiffness), the presence of the structure reduces the displacement of the point under consideration, and its behavior becomes more monotonic. In the latter case c, the low-frequency finite resonances of the structure are more pronounced, and amplitudes increase.

Data are also obtained on the displacement of points within the base at the layer interfaces beneath the structure. The results are presented in Figure 4. Plots a–c correspond to the interface between the upper and lower layers, while plots d–f correspond to the interface between the lower layer and the half-space.

A comparison of the plots leads to the conclusion that vertical displacements depend significantly, first and foremost, on the structure of the foundation soil. The presence of a structure introduces additional displacement maxima, particularly in the case of soft surface layers.

The calculations performed proved to be highly time-consuming, which necessitated consideration of parallel computing. Using the proposed method, parallelization issues in constructing the fundamental solutions for the BEM were examined.

During the computation of fundamental solutions, the stress vector was calculated 3,096,676 times; debugging tools recorded a total execution time of 56.1 seconds for the software module, with a processor time of 55 seconds.

We examine the implementation of various parallel computing techniques in a Fortran program, utilizing three technologies: Fortran Coarray, Message-Passing Interface (MPI), and OpenMP.

Coarray is a technology in which a Fortran program is replicated across multiple processors (images). This technology was introduced in Fortran 2008 to support parallel programming. It provides a simple and intuitive syntax for performing parallel computations without explicit message passing. In Fortran, Coarray is a variable declared with one or more additional codimensions. Each image maintains its own local copy of the Coarray data [31].

MPI (Message Passing Interface) is a standardized and portable message-passing system designed to facilitate communication and coordination between parallel processes in a distributed-memory parallel computing environment. MPI is commonly used in high-performance computing (HPC) applications to enable communication between multiple processes running on different cluster nodes. In Fortran, MPI offers a set of subroutine calls for direct communication among parallel processes [32].

OpenMP (Open Multi-Processing) is an API that supports cross-platform shared-memory parallel programming in C, C++, and Fortran. OpenMP simplifies the development of parallel programs by providing a set of compiler directives and library routines for expressing parallelism in code [33]. This technology is particularly well-suited for loop-level and region-level parallelization. OpenMP manages thread safety, including the handling of private and shared variables. Drawbacks of OpenMP include the need to account for potential data race conditions and the requirement for manual synchronization.

Modifications were made to the original software solution to enable parallel computing. During the experiment, the workload was distributed across 2, 4, 6, and 8 threads (images). The results are presented in Tables 3 and 4 and in Figure 5.

Discussion. The results obtained confirm theoretical concepts regarding the interaction between a structure, including its foundation, and the surrounding soil mass. An analysis of how vertical displacements at a point within a multilayer foundation depend on the layer configuration and the presence of a structure with a partially embedded foundation nearby yielded the expected results. In a configuration where layers of stiffer materials overlie a highly compressible half-space, the introduction of a structure with a partially embedded foundation, which possesses even greater stiffness, reduces vibrations at adjacent points caused by surface loading. A similar result is observed in a low-compressibility medium of uniform stiffness. Conversely, in a configuration where a rigid half-space is overlain by layers of elastic soil, the presence of a structure can intensify resonance phenomena, particularly at low frequencies.

Analyzing the displacement of points along the segment within the domain under consideration made it possible to evaluate the impact of the examined factors, including spatial considerations, on the stress-strain state of the medium. As expected, the greatest displacements were observed in the vicinity of the load application point and gradually decreased with increasing distance. The frequencies selected for analysis — 18 Hz and 22 Hz — represented the natural frequencies of the structure and were determined through modal analysis.

Given that solving this problem is a resource-intensive and time-consuming process, employing parallel computing technologies to optimize calculations represents a natural stage in the software product evolution. An analysis of the time required to compute fundamental solutions using Coarray, MPI, and OpenMP technologies reveals a high degree of similarity in performance gains across all three methods, demonstrating the maturity of these technologies within the Fortran language ecosystem. Transitioning from a sequential implementation to a two-thread setup results in a reduction of calculation time by approximately 30%. Increasing the number of threads by two results in a similar relative performance gain. This trend continues at higher thread counts. Thus, while increasing the number of threads, and consequently the number of active processor cores, reduces computation time. In absolute terms the time difference for a large number of threads will be small. This is evidenced by the growth in processor time consumed as the number of threads increases for both MPI and OpenMP, which is associated with higher synchronization overhead for parallel processes. Differences in the runtime behavior of the Coarray-based software module may stem from the increased time required to initialize parallel images or, potentially, from the specific execution characteristics of the software on the given hardware. Previous studies indicate that further increasing the number of Coarray images will result in processor time costs comparable to those of other technologies.

Conclusion. Based on the results obtained, it can be concluded that the matrix-based technique for constructing fundamental solutions for multilayer media containing surface and embedded objects is an efficient tool. It enables the accurate accounting of the medium dispersion properties and energy radiation conditions at infinity during the interaction between the object and the foundation. It also allows for the natural parallelization of computational processes through the combined use of boundary integral equation and finite element methods. The application of Coarray, MPI, and OpenMP technologies for parallel computing makes it possible to significantly reduce the software module execution time without incurring a substantial increase in resource consumption. The use of various technologies yields a similar increase in the software module performance, while differences in CPU time consumption stem from the specific implementation details of the parallelization technologies. Selecting the optimal technology depends on the specific implementation of the application and the constraints of the selected platform. Thus, when the developed module is compiled as a DLL library and integrated into a comprehensive solution, it should be noted that the Coarray technology requires the base application to be implemented in Fortran as well, while MPI applications are initialized via the mpiexec utility. Alternatives to the technologies discussed include OpenCL, OpenACC, and CUDA, which enable efficient computation using graphics accelerators. Investigating methods for applying these technologies to the analysis of complex structural systems and multilayered foundations represents a priority area for future research.