Open AccessArticle

A Study on Adaptive Implicit–Explicit and Explicit–Explicit Time Integration Procedures for Wave Propagation Analyses

Structural Engineering Department, Federal University of Juiz de Fora, Juiz de Fora CEP 36036-330, MG, Brazil

LAMEMO, PEC-COPPE, Federal University of Rio de Janeiro, Rio de Janeiro CEP 21941-611, RJ, Brazil

Author to whom correspondence should be addressed.

Acoustics 2024, 6(3), 651-680; https://doi.org/10.3390/acoustics6030036

Submission received: 2 May 2024 / Revised: 12 July 2024 / Accepted: 16 July 2024 / Published: 23 July 2024

Browse Figures

Versions Notes

Abstract

This study delves into the effectiveness of two time integration techniques, namely the adaptive implicit–explicit (imp–exp) and explicit–explicit (exp–exp) methods, which stand as efficient formulations for tackling intricate systems characterized by multiple time scales. The imp–exp technique combines implicit and explicit procedures by employing implicit formulations for faster components and explicit calculations for slower ones, achieving high accuracy and computational efficiency. Conversely, the exp–exp method, a variation of explicit methods with sub-cycling, excels in handling locally stiff systems by employing smaller sub-steps to resolve rapid changes while maintaining stability. For both these approaches, numerical damping may be activated by adaptive time integration parameters, allowing numerical dissipation to be locally applied, if necessary, as a function of the considered discrete model and its computed responses, enabling a highly effective numerical dissipative algorithm. Furthermore, both these techniques stand as very simple and straightforward formulations as they rely solely on single-step displacement–velocity relations, describing truly self-starting procedures, and they stand as entirely automated methodologies, requiring no effort nor expertise from the user. This work provides comparative studies of the adaptive imp–exp and exp–exp approaches to assess their accuracy and efficiency across a wide range of scenarios, with emphasis on geophysical applications characterized by multiscale problems, aiming to establish under which circumstances one approach should be preferred over the other.

Keywords:

implicit–explicit; explicit–explicit; hybrid analysis; adaptive parameters; sub-cycling; geophysical applications

1. Introduction

Time integration is a crucial procedure in the numerical simulation of dynamic systems, allowing for the resolution of ordinary and partial differential equations that describe the temporal evolution of several different models. Among the various time integration strategies available, two that have garnered significant attention are implicit–explicit and explicit–explicit methodologies [1,2,3,4,5,6,7,8,9,10]. These combined techniques have been shown to be effective in handling complex systems with multiple time scales, where traditional methods may struggle. Furthermore, these procedures have been applied in fields such as fluid dynamics, where they can accurately simulate interactions between the fluid and its surroundings; weather forecasting, where they can handle complex systems with various atmospheric phenomena such as air pressure, temperature, and humidity; and geophysics, where they can model the Earth’s magnetic field or the behavior of seismic waves in media with diverse physical properties, among others.

The literature contains many well-known explicit and implicit algorithms for time-marching analysis [11,12]. Explicit procedures [13,14,15,16] are often favored due to their lower computational effort but they have limitations in terms of stability conditions. On the other hand, implicit approaches [17,18,19,20] may offer unconditional stability but at the cost of needing increased computational resources. The initial studies on implicit–explicit time-marching techniques focused on combining different time integration procedures, such as the central difference (CD) method and the trapezoidal rule (TR), yet their implementation was somewhat complex and cumbersome [1]. Over the past few decades, more advanced integration procedures have been proposed, simplifying the handling of implicit and explicit subdomains and allowing for more effective combined routines [2,3,4]. Gravouil, Combescure, and other researchers have made substantial contributions to this field through their works on different time-marching integration schemes [21,22,23]. Indeed, ongoing research in this area has led to the development of various hybrid time-marching techniques for transient analyses. However, although adaptive approaches now usually consider spatial discretization techniques [24] as well as many other applications, this is not the case regarding hybrid time integration procedures. A recent hybrid methodology for time integration, which is one of the focuses of the present study, was developed by Soares [5] and it includes adaptive time integrators that aim to combine the strengths of both explicit and implicit methods, characterizing an adaptive implicit–explicit (imp–exp) scheme. Later on, by considering only the explicit counterpart of this hybrid methodology, an adaptive explicit–explicit (exp–exp) scheme was developed [9], combining this formulation with adaptive time-steps/sub-cycling splitting procedures.

The imp–exp technique merges the advantages of both implicit and explicit methods to solve differential equations. It is particularly useful when dealing with systems that exhibit multiple time scales, as it allows for appropriate treatment of both slow and fast-evolving parts of the system. By using an implicit method for the fast parts and an explicit method for the slow parts, the imp–exp technique can achieve high computational efficiency while maintaining accuracy and stability. On the other hand, the explicit technique with sub-cycling is a variation of the explicit method that uses smaller time sub-steps to accurately resolve parts of the system that exhibit rapid evolution. This approach can be highly effective in dealing with locally stiff systems, where traditional explicit methods may struggle due to stability constraints. When using sub-cycling, the explicit technique can maintain stability while still achieving high accuracy and efficiency.

In other words, under the imp–exp approach, the method incorporates an implicit operator in the stiffer regions where highly restricted time steps (computed per element) akin to a fully explicit technique would be necessary. Conversely, an explicit operator is employed in the more pliable regions, allowing for larger critical time steps. However, it is worth noting that Soares’ methodology [5] does not specifically emphasize the “coupling” of explicit and implicit routines. Instead, its main objective is to achieve a hybrid explicit and implicit formulation using a single time-marching algorithm. To accomplish this goal, Soares proposes a model/solution-adaptive imp–exp time integration procedure. This approach allows for the adaptation of the time integration based on the specific characteristics of the model and its solution responses. By dynamically and locally adjusting the time integration procedure, the methodology offers a more flexible and efficient analysis of the system, ensuring that the explicit and implicit components of the hybrid formulation are automatically distributed throughout the model and optimizing the overall performance of the algorithm.

In the exp–exp approach [9], while the rigid regions entail small time steps due to stability criteria, sub-cycling enables the more flexible regions to be modeled with larger time steps, treating the entire system in a genuinely explicit manner. This approach subdivides the model into multiple explicit subdomains with distinct time steps, with the “stiffer” subdomains having smaller time steps and the less “stiff” ones having larger time steps, always respecting stability criteria. By following this configuration, the time steps and time integrators of the time-marching technique are locally computed, taking into account the adopted spatial discretization, the properties of the model, and even the evolution of the computed responses. Both the adaptive imp–exp and exp–exp techniques offer powerful tools for accurately simulating complex models and their use has been shown to be highly effective in a wide range of applications. However, establishing which approach may be more appropriate for a given model configuration is still not clear; this paper aims to enlighten this discussion.

By focusing on wave propagation problems, significant research efforts have been directed toward the development of time integration algorithms capable of introducing controlled numerical dissipation to eliminate spurious high-frequency contributions in the response [13,14,15,16,17,18,19,20]. These spurious modes generate non-physical oscillations, which, if left unattended, can significantly impact the entire analysis, leading to deviations in results from the true physical response of the problem. The challenge resides, however, in introducing high-frequency dissipation without excessive damping in crucial low-frequency modes. One way for estimating numerical dissipation is by assessing the spectral radius of the adopted methodology, which corresponds to the magnitude of the largest eigenvalue of the amplification matrix of the referred time integration algorithm. In a stable analysis, the spectral radius ranges from 0 to 1, where a value of 1 signifies a non-dissipative algorithm. When the spectral radius is less than 1, numerical dissipation occurs and it is maximized as the spectral radius approaches zero. The current literature provides a variety of time-marching procedures with dissipative characteristics tailored specifically for hyperbolic problems [13,14,15,16,17,18,19,20,21]; however, most of these dissipative approaches are highly dependent of the adopted discretizations to become effective, a drawback that the here-discussed adaptive formulations aim to surpass.

The discussed imp–exp methodology incorporates two time integration parameters, namely α and γ. The calculation of the γ parameter enhances accuracy and ensures stability, defining the explicit and implicit components of the model. On the other hand, the discussed exp–exp approach just considers a unique time integration parameter on its formulation, the above-referred α. The evaluation of the α parameter, in both imp–exp and exp–exp time integration procedures, focuses on enabling an effective numerical dissipative algorithm. This parameter aims to eliminate the influence of spurious modes and reduce amplitude decay errors, defining the dissipative and non-dissipative elements of the model, which are updated at each time step of the analysis. The referred adaptive strategy is non-iterative, allowing the related parameters to be directly calculated considering the physical and geometrical properties of the elements of the adopted spatial discretization, the applied time-step value, and results computed from previous time steps. Furthermore, these techniques rely solely on single-step displacement–velocity relations, standing as simple truly self-starting approaches.

The present manuscript is structured as follows: firstly, the governing semi-discrete equation of the model is presented, along with the studied adaptive solution methodologies, highlighting the main aspects of both imp–exp and exp–exp techniques. In this sense, an overview of how these methodologies handle the referred semi-discrete equation and adapt themselves based on the model’s characteristics and responses is reported. Next, numerical applications are considered, taking into account theoretical benchmarks and highly complex applied models. These applications aim to discuss and compare the use of the two discussed approaches in different scenarios, analyzing their performance, accuracy, and flexibility. The considered numerical examples showcase the capabilities and limitations of each methodology, providing valuable insights into their practical applications. Finally, the paper concludes with a summary of the properties of the discussed formulations, presenting the key findings and positive aspects of both approaches. The conclusions highlight the strengths and potential applications of the methodologies, giving an overview of their suitability for various wave propagation problems in engineering and geophysics.

2. Methods

The semi-discrete equation for a hyperbolic transient system may be written as:

M \ddot{U} (t) + C \dot{U} (t) + K U (t) = F (t)

(1)

where

M

C

, and

K

represent the mass, damping, and stiffness matrices of the model, respectively,

F (t)

represents the applied force vector, and

\ddot{U} (t)

\dot{U} (t)

, and

U (t)

describe the vectors representing acceleration, velocity, and displacement, respectively. The initial conditions of this model are defined by

U^{0} = U (0)

as the initial displacement vector, and

{\dot{U}}^{0} = \dot{U} (0)

as the initial velocity vector.

The standard Finite Element Method (FEM) is considered herein for spatial discretization. However, the time integration formulations that are studied in this paper are not limited to be exclusively used in conjunction with this specific method. In fact, the discussed time integration procedures are expected to be effective considering any spatial discretization technique that is based on local formulations. Leveraging the conventional FEM, the domain of the problem is subdivided into elements, thereby enabling the computation of local matrices and vectors. As a result, through the integration of Equation (1) over time, at the element level (which is indicated by the subscript “

e

”), taking into account a time increment

Δ t

(where

t^{n + 1} = t^{n} + Δ t

), the following expression can be established:

M_{e} \int_{t^{n}}^{t^{n + 1}} {\ddot{U}}_{e} (t) d t + C_{e} \int_{t^{n}}^{t^{n + 1}} {\dot{U}}_{e} (t) d t + K_{e} \int_{t^{n}}^{t^{n + 1}} U_{e} (t) d t = \int_{t^{n}}^{t^{n + 1}} F_{e} (t) d t

(2)

whose integrals may be evaluated as:

\int_{t^{n}}^{t^{n + 1}} {\ddot{U}}_{e} (t) d t = {\dot{U}}_{e}^{n + 1} - {\dot{U}}_{e}^{n}

(3a)

\int_{t^{n}}^{t^{n + 1}} {\dot{U}}_{e} (t) d t = U_{e}^{n + 1} - U_{e}^{n}

(3b)

\int_{t^{n}}^{t^{n + 1}} U_{e} (t) d t = ∆ t U_{e}^{n} + ½ α_{e}^{n} ∆ t^{2} {\dot{U}}_{e}^{n} + ½ γ_{e}^{n} ∆ t^{2} {\dot{U}}_{e}^{n + 1}

(3c)

{\bar{F}}_{e} = \int_{t^{n}}^{t^{n + 1}} F_{e} (t) d t

(3d)

It can be observed that, in Equation (3c), there are two time integration parameters, which are denoted by

α_{e}^{n}

and

γ_{e}^{n}

. The values of these parameters may vary for each element “

e

” within the adopted spatial discretization and over each discrete time step “

n

” of the considered time-marching methodology. Additionally, in Equation (3d),

\bar{F}

represents the time integral of the force term, which may be analytically evaluated or computed by employing any suitable standard numerical procedure.

By considering the definitions specified by Equations (2) and (3a–d), the recurrence relations that are described in Equation (4a,b) can be obtained to solve Equation (1), establishing a time-domain solution algorithm that is based solely on displacement–velocity relationships, as follows:

(M_{e} + ½ Δ t C_{e} + ½ γ_{e}^{n} Δ t^{2} K_{e}) {\dot{U}}_{e}^{n + 1} = {\bar{F}}_{e} + M_{e} {\dot{U}}_{e}^{n} - ½ Δ t C_{e} {\dot{U}}_{e}^{n} - K_{e} (Δ t U_{e}^{n} + ½ α_{e}^{n} Δ t^{2} {\dot{U}}_{e}^{n})

(4a)

U^{n + 1} = U^{n} + ½ Δ t {\dot{U}}^{n} + ½ Δ t {\dot{U}}^{n + 1}

(4b)

In this case, as one may observe, Equation (4a) allows computation of the velocities of the model once assembling is carried out, whereas Equation (4b) allows the calculation of its displacements. As Equation (4a) indicates, if

γ = 0

is considered and if lumped mass and damping matrices are regarded (as is the case in this work), a diagonal effective matrix takes place, turning the governing system of equations explicit, which highly reduces the computational effort of the indicated solution procedure. Thus, when

γ_{e}^{n}

is defined by a non-null value, it specifies a so-called implicit element and when it is defined by a null value, it specifies an explicit element. In this context, null and non-null values may occur for

γ_{e}^{n}

, taking into account the focused imp–exp formulation; just null values are applied for this parameter (and it may then be disregarded) once the exp–exp approach is considered.

As previously highlighted,

α_{e}^{n}

and

γ_{e}^{n}

are adaptively calculated herein, taking into account the local features of the model. In the next subsection, in which the focused imp–exp approach is discussed, the expressions to locally evaluate these parameters are provided. Thus, for the exp–exp approach, which is discussed in the subsequent subsection, the same equation that is below-presented to calculate

α_{e}^{n}

can be employed, considering

γ_{e}^{n} = 0

on its formulation.

2.1. Implicit–Explicit Approach

As mentioned before, the

γ_{e}^{n}

parameter determines whether explicit or implicit subdomains are defined within the model. Explicit elements are created when

γ_{e}^{n} = 0

, while implicit elements are generated when

γ_{e}^{n} \neq 0

. Implicit elements yield higher accuracy when

0 < γ_{e}^{n} < 1 / 2

is considered. Specifically, by setting

α = 1

and

γ = 0

, the proposed technique becomes spectrally equivalent to the central difference (CD) method. Similarly, when

α

= 1/2 and

γ

= 1/2, it emulates the trapezoidal rule (TR). Therefore, for

0 < γ < 1 / 2

and

α = 1 - γ

(representing a non-dissipative formulation), an intermediate methodology between the CD and the TR emerges. Since the CD results in negative period elongation and the TR leads to positive period elongation, utilizing

0 < γ_{e}^{n} < 1 / 2

helps minimize period elongation errors, contributing to a more accurate technique.

The value of

γ_{e}^{n}

is established based on the maximum sampling frequency of element

Ω_{e}^{m a x}

. This sampling frequency is computed as

Ω_{e}^{m a x} = ω_{e}^{m a x} Δ t

, where

ω_{e}^{m a x}

represents the highest natural frequency of the element (i.e., it is the highest square root of the eigenvalues

λ_{e}

that are calculated based on the local matrices

M_{e}

and

K_{e}

, solving the generalized eigenvalue problem

K_{e} θ_{e} = {λ_{e} M}_{e} θ_{e}

). In this context,

γ_{e}^{n}

may be evaluated as follows [5]:

I f Ω_{e}^{m a x} \leq 2, γ_{e}^{n} = 0

(5a)

I f Ω_{e}^{m a x} > 2, γ_{e}^{n} = \frac{1}{2} \tanh (\frac{1}{4} Ω_{e}^{m a x})

(5b)

where, as specified by Equations (5a,b), the critical sampling frequency for the considered explicit formulation is 2 (similar to that of the CD method) and

0 < γ_{e}^{n} < 1 / 2

for the considered implicit formulation. It is important to observe that Equation (5b) is deduced so that the critical sampling frequency of the method always becomes greater than

Ω_{e}^{m a x}

, ensuring stability for the implicit elements of the model. Consequently, the combination of Equation (5a,b) guarantees stability for all of the elements, both explicit and implicit.

One should observe that, in cases where the element properties do not vary throughout the analysis, as in standard linear models, the values of

γ_{e}^{n}

remain unchanged over time. In fact, if the element matrices remain unaltered, the time integration parameter

γ_{e}^{n}

also remains the same along time, resulting in an unaltered effective matrix throughout the analysis (avoiding its updating and computationally demanding correspondent treatment). Moreover, as previously observed, the incorporation of explicit elements into the analysis leads to a portion of the effective matrix consisting solely of diagonal entries, which can be readily treated apart from the global system of equations. This reduction in dimensionality and computational effort simplifies the solving process of the system. It is important to emphasize that this automatic reduction in the dimension of the global system is achieved by considering only the diagonal terms of the effective matrix, devoid of any pre-established subdomains or user-supplied input data. Furthermore, the above-described division into implicit and explicit subdomains transpires seamlessly, requiring solely the computation of

Ω_{e}^{m a x}

, as indicated by Equations (5a,b), thus highlighting the straightforward implementation process of the discussed approach.

As previously remarked upon, the time integration parameter

α_{e}^{n}

governs the dissipative characteristics of the technique; it can be adjusted based on the solution’s evolution. Its purpose is to introduce numerical dissipation when and where it may be necessary in order to mitigate the presence of undesired non-physical oscillations. The local computation of

α_{e}^{n}

aims to optimize the incorporation of numerical damping into the analysis while minimizing its adverse effects. For

α_{e}^{n} = 1 - γ_{e}^{n}

, no numerical dissipation is introduced into the analysis, whereas, for

α_{e}^{n} > 1 - γ_{e}^{n}

, numerical damping is applied. This feature allows for the selective activation of dissipation in specific regions and/or time periods where/when oscillations occur and deactivation in regions and/or time periods where/when they are not observed. By doing so, an overall algorithmic dissipative pattern can be avoided, thereby preventing excessive numerical damping errors.

The activation of dissipation can be locally carried out based on an oscillatory criterion. In this context, the

α_{e}^{n}

parameter of the elements surrounding an oscillating degree of freedom (i.e., a wavering

u_{i}

value) can be adjusted to introduce numerical dissipation. On the other hand, if no oscillatory behavior is observed,

α_{e}^{n}

may be set to

1 - γ_{e}^{n}

. This process involves computing an oscillatory parameter

φ_{e}^{n}

for each element and time step of the analysis. If

φ_{e}^{n} = 0

(indicating no oscillation),

α_{e}^{n}

is assigned the value of

1 - γ_{e}^{n}

. Otherwise, if

φ_{e}^{n}

is non-zero,

α_{e}^{n}

is set to a value greater than

1 - γ_{e}^{n}

, as shown in the next equations [5]:

φ_{e}^{n} = \sum_{i = 1}^{η_{e}} ||u_{i}^{n} - u_{i}^{n - 2}| - |u_{i}^{n} - u_{i}^{n - 1}| - |u_{i}^{n - 1} - u_{i}^{n - 2}||

(6a)

I f φ_{e}^{n} = 0, α_{e}^{n} = 1 - γ_{e}^{n}

(6b)

I f φ_{e}^{n} \neq 0, α_{e}^{n} = 2 {[2 γ_{e}^{n} + (1 + \frac{ς_{e} Δ t}{2 ρ_{e}}) {(\frac{2}{Ω_{e}^{m a x}})}^{2}]}^{1 / 2} - 1 - γ_{e}^{n} - \frac{ς_{e} Δ t}{2 ρ_{e}} {(\frac{2}{Ω_{e}^{m a x}})}^{2}

(6c)

where

ρ_{e}

and

ς_{e}

, which stand for the mass density and viscous damping coefficient of the element, respectively, represent the physical properties of the medium that are associated with matrices

M_{e}

and

C_{e}

, respectively, and

η_{e}

stands for the total amount of degrees of freedom of the element. Equation (6c) is constructed in such a way that it allows for maximum numerical dissipation to be applied at the highest frequency of the element (i.e., at

Ω_{e}^{m a x}

). Thus, this approach becomes highly effective in mitigating the effects caused by spurious high-frequency modes, dissipating their influences to a significant extent without significantly affecting the contribution of the important low-frequency modes. The oscillatory criterion described by Equation (6a) evaluates if an oscillation takes place considering results computed at the last three time steps. In this case, if the solution increment within two consecutive time steps (i.e.,

|u_{i}^{n} - u_{i}^{n - 2}|

) is not equal to the sum of the corresponding increments within one time step (i.e.,

|u_{i}^{n} - u_{i}^{n - 1}| + |u_{i}^{n - 1} - u_{i}^{n - 2}|

), an oscillation occurs.

For the discussed imp–exp analysis, a crucial consideration revolves around the selection of an appropriate time-step value (Δt). In this sense, the decision to increase or not the Δt value not only affects the accuracy of the considered solution procedure but also involves a delicate balance between reducing (or not reducing) the total amount of time steps required for the analysis and introducing (or not introducing) more implicit elements into the effective matrix, which directly impacts the computational efficiency of the hybrid approach. Thus, in order to find an optimal time-step value that provides maximal computational efficiency, the employment of an optimization algorithm becomes essential. For this goal, in the present study, the Particle Swarm Optimization (PSO) algorithm [25] is employed, which automatically determines the optimal Δt value by minimizing the expected total number of operations involved in the referred hybrid solution process. This approach helps to streamline the solution process by optimizing the time-step selection, ultimately improving the overall efficiency of the referred imp–exp analysis. In this context, the imp–exp technique always becomes more efficient than its purely explicit or purely implicit counterparties. Moreover, this approach also provides an entirely automated formulation, in which all the parameters of the discussed time integration procedure are automatically optimally computed, requiring no effort nor expertise from the user.

2.2. Explicit–Explicit Approach

Focusing on the explicit counterpart of the methodology described by Equations (4a,b), for which

γ_{e}^{n} = 0

and no implicit subdomains are created, the entire effective matrix of the method becomes diagonal, which is a distinctive characteristic of explicit approaches that renders very efficient calculations per time step. However, once this procedure is followed and the

ω_{e}^{m a x}

values of the model are computed, the obtained maximal

ω_{e}^{m a x}

value should be employed to establish the (locally evaluated) critical time step of the analysis, so that stability may be granted. This implies that, if a significantly “stiff” region exists in the model, due to its physical properties or mesh refinement level, this critical time step may become exceedingly small, leading to a substantial escalation in computational costs. To address this problem and enable efficient analyses considering generic (and even challenging) models, explicit techniques may consider locally defined time-step values associated with sub-cycling splitting procedures. This concept is adopted here in association with the explicit counterpart of the solution algorithm (4a,b), defining the referred exp–exp scheme, which subdivides the model into different explicit subdomains with distinct time-step values.

This multiple-time-step approach optimizes the computational effort by adjusting the time step according to the stiffness variations within the model, ultimately leading to a more efficient explicit analysis. In this way, “stiffer” subdomains may employ smaller time steps, while “softer” regions may consider larger time-step values, according to the stability limits of the different subdomains of the model. Moreover, through the implementation of adaptive time steps and sub-cycling splitting routines, the exp–exp scheme allows each computed subdomain to be analyzed, considering locally defined time integration parameters that may then be more properly evaluated. Thus, this implementation not only enhances the efficiency of the referred solution procedure but also improves its accuracy. For the exp–exp methodology, the

α_{e}^{n}

parameters of the method are computed as described by Equations (6a–c), considering

γ_{e}^{n} = 0

. In this case,

α_{e}^{n} > 1

is applied wherever and whenever and numerical damping may be necessary;

α_{e}^{n} = 1

is considered otherwise. As discussed before, by following Equation (6c), maximum numerical damping is imposed at the highest sampling frequency of the elements, resulting in an effective dissipative approach for the highest modes of the model, which are then “tracked” by the referred adaptive methodology.

In this study, an automated algorithm is considered to partition the model’s domain into multiple time-step subdomains, requiring no effort nor expertise from the user to accomplish this subdivision. The referred algorithm undertakes this task by computing and allocating a specific time-step value to each node within the model. This procedure encompasses the division of the model into sets of elements that can share the same assigned time-step value, ensuring compliance with their respective stability limits. The subsequent sequence of commands is employed to automatically establish this subdomain division: (i) first, the critical time steps of all elements are calculated, as indicated by

Δ t_{e} = 2 / ω_{e}^{m a x}

, and the smallest

Δ t_{e}

in the model, denoted as

{Δ t}_{e}^{m i n}

(where

{Δ t}_{e}^{m i n} = m i n (Δ t_{e})

), is identified as the basic time step for the controlled subdivision of the domain; (ii) with

{Δ t}_{e}^{m i n}

established, subsequent time-step values are calculated as multiples of powers of 2 of this minimal time-step value (i.e.,

Δ t_{i}

values are calculated, where

Δ t_{i} = 2^{(i - 1)} {Δ t}_{e}^{m i n}

); (iii) each element is then associated with a computed time-step value (

Δ t_{i}

), where

Δ t_{i} \leq Δ t_{e} \leq Δ t_{i + 1}

and “i” indicates the subdomain to which that element belongs; and (iv) finally, each degree of freedom in the model is assigned a time-step value (subdomain) based on the lowest

Δ t_{i}

value among its surrounding elements. By adhering to these steps, the discussed automated algorithm provides a suitable division of the domain, assigning appropriate time-step values to elements and degrees of freedom. This methodology enhances the efficiency of the analysis while upholding accuracy, as it permits a customized time-step distribution based on the local characteristics of the model.

Upon implementing this subdomain division/sub-cycling algorithm, it becomes essential to tackle the matter of interpolating displacement and velocity values near the boundaries of these multiple time-step subdomains. In this study, the following expressions are employed for these interpolation tasks, which are consistent with the adopted time integration procedure, as follows:

U (t) = ({\dot{U}}^{n + 1} - {\dot{U}}^{n}) t^{2} / (2 Δ t) + {\dot{U}}^{n} t + U^{n}

(7a)

\dot{U} (t) = ({\dot{U}}^{n + 1} - {\dot{U}}^{n}) t / Δ t + {\dot{U}}^{n}

(7b)

In Equations (7a,b),

t

stands for the current increment of time (

0 \leq t \leq Δ t

) for the focused subdomain and

∆ t

is the time-step value of the degree of freedom being interpolated, which is related to the neighboring subdomain.

In Figure 1, a schematic representation of this interpolation procedure and sub-cycling computations are depicted. Figure 1a illustrates two subdomains: subdomain 1, characterized by a time-step Δt, and subdomain 2. When calculating displacements within subdomain 2, it becomes necessary to incorporate results from nodes located outside this subdomain that are within elements that share nodes with subdomain 2 (as indicated in the figure). These values are obtained through interpolation, utilizing the computed responses from subdomain 1, as depicted in Figure 1a, and considering Equation (7a,b). By adopting this approach, the explicit methodology enabled by Equation (4a,b) can be seamlessly applied to each subdomain of the model, as indicated in Figure 1b.

As one may observe, in contrast to the imp-exp methodology, the exp-exp approach considers a less straightforward subdomain division and a more complex multi-domain approach, needing to map and interpolate degrees of freedom that are connected to subdomain boundaries etc. However, this approach may allow several different subdomains to take place at once, according to the properties of the discretized model, whereas the reported imp-exp formulation only permits a maximal of two subdomains to occur (i.e., the implicit and the explicit subdomains). Thus, greater flexibility is provided by this approach, which may result in more effective analyses.

2.3. Properties of the Methods

To study the properties of the referred techniques, following standard guidelines, a single-degree-of-freedom problem may be considered and the correspondent amplification matrices that arise for each methodology may be investigated [5]. In Figure 2 and Figure 3, the spectral radii of these amplification matrices (which corresponds to the largest modulus of the eigenvalues of these matrices) are provided, considering the implicit and the explicit counterparts of the discussed solution algorithm described by Equations (4a,b), respectively. As previously remarked and as one may observe in Figure 2a, Equation (5b) is formulated so that the critical sampling frequency of the method always becomes greater than

Ω_{e}^{m a x}

, thus ensuring stability for the implicit elements of the model. As also previously highlighted, Equation (6b) describes a non-dissipative approach and, as so, it provides a unity value for the spectral radius of the method, within its stable domain. By following Equation (6c), on the other hand, the spectral radius may become lower than 1, introducing numerical dissipation into the analysis. In this case, as is illustrated in Figure 2b, Equation (6c) ensures that the bifurcation sampling frequency of the method equals

Ω_{e}^{m a x}

, applying maximal numerical dissipation at the highest frequency of the element, which describes an optimal dissipative approach. Analogous features may be observed considering Figure 3. In this case, as previously highlighted and as Figure 3a illustrates, the discussed methodology becomes spectrally equivalent to the CD once Equations (5a) and (6b) are followed (i.e., once

γ = 0

and

α = 1

For the discussed imp–exp approach, implicit elements only occur when

Ω_{e}^{m a x} > 2

; thus, Figure 2 is depicted only considering

Ω_{e}^{m a x}

values greater than or equal to 2. It is also important to observe that, for the exp–exp approach, since

Δ t_{i} = 2^{(i - 1)} {Δ t}_{e}^{m i n}

is considered when specifying the time-step subdomains of the model and

{Δ t}_{e}^{m i n} = 2 / m a x (ω_{e}^{m a x})

Ω_{e}^{m a x}

will never be less than one in this case;

1 < Ω_{e}^{m a x} \leq 2

stands for all the elements of the model following this solution strategy. As one may observe upon comparing Figure 2b and Figure 3b, the exp–exp approach tends to enable a more intense dissipative approach than the imp–exp methodology, since the discussed explicit formulation allows for lower values for the spectral radius of the method than the implicit formulation. Therefore, for problems in which a great amount of spurious oscillations occurs, the exp–exp technique will probably become more accurate than the imp–exp formulation. However, it is important to observe that the discussed imp–exp approach considers the optimal (in terms of efficiency) distribution of explicit and implicit elements along the model by computing an optimized time-step value, as described at the end of Section 2.1. Thus, for usual wave propagation models, a great number of explicit elements will regularly automatically take place in the analysis and the imp–exp approach will then also be very effective in solving these challenging problems.

In Figure 4 and Figure 5, standard period elongation and amplitude decay errors [5] are depicted, considering the implicit and the explicit counterparts of the discussed solution algorithm, respectively. As one may observe in Figure 4a, the discussed non-dissipative implicit approach is always more accurate than the TR, replicating this technique for

Ω_{e}^{m a x} \to \infty

(i.e., for

γ = 1 / 2

and

α = 1 / 2

). For the dissipative implicit formulation, very low period elongation errors are also observed, as indicated in Figure 4b. It is important to remark that, for the imp–exp approach, considering non-dissipative calculations (which mostly occurs, unless an oscillation is locally perceived), a counterbalance of period elongation errors may take place once the discussed implicit and explicit formulation renders positive (Figure 4a) and negative (Figure 5a) period elongation errors, respectively. Thus, this hybrid approach may become very accurate.

In the next section, wave propagation problems considering different spatially discretized models and/or material distributions are regarded so that the performance of the referred multi-domain adaptive time-marching formulations can be properly discussed. In this way, several different subdomain compositions and time-integration parameter values may occur along the analyses, allowing to properly access the effectiveness of each hybrid formulation and to establish which approach may be more appropriate considering a given set of configurations.

3. Results and Discussions

In this section, five numerical applications are studied, aiming at assessing both the accuracy and efficiency of the discussed imp–exp and exp–exp techniques. The initial focus within the first subsection presented here centers on the examination of two theoretical models for which well-established analytical solutions exist. The first of these models scrutinizes a homogeneous infinite system, whose analytical solution is described by Green’s function of the acoustic problem [26]. Leveraging the availability of this analytical solution, one may meticulously gauge the errors inherent in the numerically computed responses. The second example delves into the investigation of a heterogeneous rod comprising two distinct materials [27], providing, once again, a valuable opportunity to appraise the precision of the discussed solution methodologies and gain insights into the advantages and efficiency stemming from rigid-flexible configurations. Moving on to the latter part of this section, three synthetic models crafted to emulate geological scenarios closely resembling real-world applications are explored. These models are considered to serve as compelling demonstrations of the efficacy of the discussed techniques when confronted with the analysis of extensive geophysical and other wave propagation modeling challenges, particularly those commonly encountered in the oil and gas industry. In this context, the Buzios [28], the 2DEW [29], and the 2004 BP [30] models are studied herein. These benchmark models are distinguished by their intricate structures, replete with multiple layers housing heterogeneous properties, including significant salt regions, representing complex acoustic [28], and elastodynamic [29,30] wave propagation configurations.

To assess the numerical results obtained from the discussed imp–exp and exp–exp approaches, a comparative analysis with well-established standard implicit and explicit methods is engaged. This arrangement encompasses the classic trapezoidal rule (TR) [17], the central difference (CD) method, the implicit generalized α (IG-α) method [18] (with

ρ_{\infty}

= 0.5), the explicit generalized α (EG-α) method [13] (with

ρ_{b}

= 0.3665, as recommended by the authors), the Bathe–Baig implicit composite (IC) method [19], and the Noh–Bathe explicit composite (EC) method [14] (with p = 0.54, as recommended by the authors). Additionally, the purely implicit (imp) and explicit (exp) counterparts of the studied methodology are also included for comparison. In this work, the above referred explicit techniques are applied considering their maximum allowable time-step values for stability (which are determined through element-level evaluations), resulting in the most efficient analysis for each approach. In contrast, the referred implicit methods utilize the same time-step value that is optimally computed herein for the imp–exp solution procedure (as discussed in the last paragraph of Section 2.1), allowing comparisons that consider equivalent discretization levels as well. Through these comparative studies, our objective is to evaluate the performance of the discussed techniques in terms of accuracy and computational efficiency. This comprehensive assessment will furnish valuable insights into the applicability of the imp–exp and exp–exp schemes across various engineering applications, with a particular emphasis on their suitability for addressing large-scale geophysical models.

In the next subsection, in which analytical solutions are available, the errors of the computed responses are evaluated using the following expression:

Error = {[\sum_{n = 1}^{N} {(u^{n} - u_{A} (t^{n}))}^{2} ∕ \sum_{n = 1}^{N} {(u_{A} (t^{n}))}^{2}]}^{1 / 2}

(8)

where the computed field, representing the time history of a given degree of freedom, is denoted by

u

, its analytical counterpart is represented by

u_{A}

, and the total number of time steps in the analysis is indicated by N. The international unit system is followed in the subsections that follow.

3.1. Theoretical Models

In the first example within this subsection, an infinite acoustic model, subjected to an impulsive source, is analyzed. Four FEM meshes, with different refinement levels, are considered to discretize the model; all of them centered at the point of source application. These meshes are designed to have a higher concentration of elements in the vicinity of the source application point, as illustrated in Figure 6. The number of linear triangular elements in each mesh is as follows: (i) discretization 1–50,000 elements; (ii) discretization 2–100,000 elements; (iii) discretization 3–150,000 elements; and (iv) discretization 4–200,000 elements. These four meshes are generated following the same pattern, with dimensions of 50 × 50 m and a 5 m-thick perfectly matched layer (PML) [31] on each side (this setup is designed to prevent wave reflections at the borders of the discretized model). Since the referred medium is homogeneous, the regions of higher “stiffness” are solely determined by the refinement of these meshes.

Figure 7 displays the computed results for

Ω_{e}^{m a x}

and

γ_{e}^{n}

by considering the four referred meshes and the discussed imp–exp approach and by demonstrating its simplicity to determine the implicit (colored,

γ_{e}^{n} \neq 0

) and explicit (white,

γ_{e}^{n} = 0

) subdomains of the model. For these different spatial discretizations, different optimal time-step values are computed, and the model is then subdivided as follows, considering the percentages of elements in each subdomain: (i) discretization 1—90.36% explicit and 9.64% implicit, for

Δ t = 0.14926

s; (ii) discretization 2—90.78% explicit and 9.22% implicit, for

Δ t = 0.10052

s; (iii) discretization 3—88.72% explicit and 11.28% implicit, for

Δ t = 0.07593

s; and (iv) discretization 4—88.45% explicit and 11.55% implicit, for

Δ t = 0.06203

Figure 8, on the other hand, illustrates the computed results for

Δ t_{e}

and

Δ t_{i}

, by considering the four referred meshes and the discussed exp–exp approach and by demonstrating its versatility to determine the multiple time-step subdomains of the model. In this case, the obtained divisions may be characterized as follows, considering the percentages of elements for each subdomain: (i) discretization 1 (4 subdomains)—1.52% for

Δ t_{1} = 0.04603 s

, 18.12% for

Δ t_{2} = 0.09207 s

, 80.32% for

Δ t_{3} = 0.18415 s

, and 0.04% for

Δ t_{4} = 0.36830 s

; (ii) discretization 2 (4 subdomains)—0.73% for

Δ t_{1} = 0.02337 s

, 7.12% for

Δ t_{2} = 0.04674 s

, 35.18% for

Δ t_{3} = 0.09348

s, and 56.97% for

Δ t_{4} = 0.18696 s

; (iii) discretization 3 (5 subdomains)—0.45% for

Δ t_{1} = 0.01249 s

, 4.28% for

Δ t_{2} = 0.02499 s

, 14.69% for

Δ t_{3} = 0.04998 s

, 80.36% for

Δ t_{4} = 0.09996 s

, and 0.22% for

Δ t_{5} = 0.19993 s

; and (iv) discretization 4 (5 subdomains)—

0.33 %

for

Δ t_{1} = 0.00811 s

, 3.05% for

Δ t_{2} = 0.01623 s

, 9.35% for

Δ t_{3} = 0.03246 s

, 30.71% for

Δ t_{4} = 0.06493 s

, and 56.56% for

Δ t_{5} = 0.12987 s

The analytical response for this model may be found in [26]. The computed time–history results at a point located 10 m horizontally away from the applied source are presented in Figure 9, for discretization 4. As one may observe, the discussed multi-domain methodologies provide considerably more accurate responses than standard techniques, as well as more adequately dissipating spurious numerical oscillations, due to their adaptive design. In fact, as this figure illustrates, they also provide better responses than their equivalent purely implicit or purely explicit approaches, indicating that not only may these hybrid procedures enhance the computational performances of these techniques but that they may also improve their accuracies.

Figure 10 displays the convergence curves for the selected time integration procedures and discretizations. As depicted in this figure, the discussed multi-domain techniques yield lower error values for discretization 1 than standard techniques yield to the much finer discretization 4, highlighting the considerably better performance of the discussed hybrid formulations. In Table 1, a detailed description of the performance of each adopted solution technique is provided, for each adopted spatial discretization. As this table indicates, the exp–exp methodology attains the highest level of accuracy across all the adopted discretizations, with error values that may become even more than three times lower than those computed by standard techniques. The multi-domain approaches also exhibit the best computational times, providing, in this case, analyses that may become more than four times faster than those performed by standard formulations. Additionally, it is interesting to observe that, for discretizations 1 and 2, lower CPU times are observed for the exp–exp methodology, whereas, for discretizations 3 and 4, better efficiency is provided by the imp–exp approach. By analyzing Table 1, one may also perceive, by comparing the data from discretizations 1 and 4, that the discussed hybrid formulations are able to provide errors that are approximately twice lower than those of standard procedures considering computational efforts more than 10 times lower (sometimes even 20 times lower) than those of these approaches, describing a huge gain in performance.

Thus, as this first example illustrates, the discussed adaptive hybrid approaches may become much more effective than standard formulations for solving wave propagation problems, especially if a great amount of spurious oscillations tend to occur on their computed responses. In this case, since the reported dissipative procedure “tracks” the higher frequency range of the model and numerical damping is only spatially and temporally applied if necessary, the discussed adaptive approaches become very effective in providing accurate responses. Moreover, as previously reported (see Section 2.3), explicit elements enable greater numerical dissipation than implicit elements, following the discussed adaptive approach; thus, the referred exp–exp formulation tends to provide more accurate responses than the imp–exp procedure once a greater amount of numerical dissipation is required by the analyzed model, as is the case here.

For the second application that is studied in this subsection, a more extensive comparison between the imp–exp and exp–exp approaches is provided. In this second example, a rectangular heterogeneous rod comprising various combinations of two distinct materials is examined. In this context, multiple model configurations are scrutinized herein, with the proportion of each material ranging (directly or inversely) from 5% to 95% of the composition of the rod. This variation results in a corresponding adjustment in the rigid region of the rod. For all the analyses that are carried out here, the so-called “material 1” maintains the same physical properties, which may be specified by an elastic modulus of

100 k P a

and density of

1 k g / m^{3}

. The density of the so-called “material 2” remains the same as that of material 1; however, its elastic modulus is systematically modified to generate progressively increasing ratios of wave propagation velocities (

c_{2} / c_{1}

) between the two materials. This leads to differences between the two considered velocities ranging from two to eight times. To spatially discretize the rod, a structured finite element mesh comprising 32,000 (400 × 80) linear square elements is considered.

The rod is fixed at its left border (

x = 0

) and subjected to a prescribed axial traction acting on its right border (

x = L

), which is defined by

P (t) = A \sin (\frac{π t}{T}) [H (t) - H (t - T)]

, where A and T describe the amplitude and the duration of the applied traction, respectively. The length of the rod is defined by

L = 1 m

and

a

(

0 \leq a \leq L

) specifies the position of the (vertical) interface between the two materials that compose the model. The analytical answer for the axial displacements of the rod may be found in [27]. As previously observed, in this study,

a

varies from 0.05 L to 0.95 L, with increments of 0.05 L, and

T = (2 / 5) [a / c_{1} + (L - a) / c_{2}]

is adopted.

Initially, the two studied multi-domain techniques are once again analyzed in conjunction with the referred classic implicit and explicit methodologies, considering

c_{2} / c_{1} = 4

and

c_{2} / c_{1} = 6

, for

a = L / 2

. In this context, time–history results of the computed axial displacements at the middle of the rod (

x = L / 2

) are depicted in Figure 11 and Figure 12 and a detailed description of the performance of each considered technique for these configurations is provided in Table 2. As one may observe in these figures and table and as is expected, particularly good results are provided by the exp–exp approach for

c_{2} / c_{1} = 4

, since, for this configuration, the discussed exp–exp formulation allows all elements of the model to be analyzed considering their critical time-step values, providing a spatial/temporal discretization for which the referred non-dissipative explicit approach (which is equivalent to the CD) is able to reproduce the analytical response of the considered model [32]. In fact, for this studied rod, which behaves like a one-dimensional application, very good responses are expected for the discussed exp–exp approach once

c_{2} / c_{1} = 2^{μ}

, for an integer

μ

value.

As is further indicated in Table 2, both multi-domain methodologies show lower computational costs than the referred standard techniques, as well as very low errors. In fact, since the imp–exp approach computes an optimal time-step value for the analysis, its CPU time is supposed to always become lower than or, at least, equal to those of equivalent purely implicit or purely explicit formulations; this aspect is described in Table 2. As indicated in this table, for

c_{2} / c_{1} = 4

, the computed optimal time-step value for the imp–exp approach equals that of the critical time step of the exp technique. Thus, all elements of the model become defined as explicit and the imp–exp approach replicates the exp procedure. On the other hand, for

c_{2} / c_{1} = 6

, the computed optimal time-step value for the imp–exp formulation becomes 6 times greater than that of the critical time step of the exp technique. In this case, the first half of the model (

0 \leq x \leq a

) automatically becomes defined by implicit elements and the second half by explicit elements; the imp–exp approach is then able to become more efficient than the exp technique, as described in Table 2. It is also interesting to observe that, for

c_{2} / c_{1} = 4

, by considering equal time-step values, the CPU time of the imp technique becomes approximately 7 times greater than that of the exp approach, clearly illustrating the amount of extra effort that the considered solver routine requires for the solution of this model.

In order to better study the performances of the discussed imp–exp and exp–exp approaches, an extensive series of comparative analyses is further conducted herein, considering various configurations for the referred heterogeneous rod. In this context, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 are presented, showcasing the computed errors and CPU time results for different material distributions along the model and

c_{2} / c_{1}

relations. As one may observe in these figures, as the size of the faster region (material 2) increases, both exp and exp–exp approaches tend to provide better accuracy since more explicit elements are then analyzed considering their critical time-step values. On the other hand, an opposite configuration takes place for the imp–exp approach. In this case, as the size of the slower region (material 1) increases, more explicit elements are introduced into the analysis considering their critical time-step values, improving the accuracy of the computed responses.

Of course, there may be oscillations in these tendencies of the computed errors due to the complexity of the solution responses that are being evaluated, which depends on the amounts of materials 1 and 2 that are considered in the analyses. To illustrate these multiple possible levels of complexity regarding the computed solution responses, time–history results for the axial displacements at the middle of the rod are depicted in Figure 20, considering three different percentages of material 2 in the model, namely 10%, 50%, and 90%. As one may observe in this figure, the complexity of the response to be evaluated greatly depends on the adopted material distribution along the rod, an aspect that further interferes with the computed error results that are provided in Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19.

As Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 further illustrate, both imp–exp and exp–exp procedures always provide more efficient analyses than the exp technique (as expected); they also tend to provide more accurate responses. It is important to observe that the imp–exp approach may become very efficient if just a small part of the domain stands as a “stiff” region, allowing most of the model to be explicitly treated while a high time-step value is considered (and this relative better performance becomes more evident as this “stiff” region becomes stiffer). On the other hand, the exp–exp approach may become extremely effective once it enables more appropriate time-step values to be automatically distributed along the model, allowing the errors of the applied spatial and temporal discretization procedures to be better counterbalanced, as well as the computational effort of the time-marching formulation to be reduced. For the referred heterogeneous rod, as previously highlighted, this is particularly the case for

c_{2} / c_{1}

equals 2, 4, and 8. For these configurations, extremely accurate responses are provided by the exp–exp approach, as indicated in Figure 13, Figure 15 and Figure 19.

One last study is carried out in this subsection, considering the referred clamped rod. In this case, a homogeneous rod is regarded (

c_{2} / c_{1} = 1

), which is spatially discretized considering an unstructured FEM mesh. The adopted mesh for the analysis of this homogeneous model, which is also composed of 32,000 linear quadrangular elements, is depicted in Figure 21. This figure also depicts the computed

Ω_{e}^{m a x}

γ_{e}^{n}

Δ t_{e}

, and

Δ t_{i}

values that arise for this discrete model. In this context, Figure 21c,e describes the obtained subdomains for the imp–exp and exp–exp analyses, respectively. For the referred imp–exp approach, 97.91% of the elements of the mesh are treated explicitly, whereas 2.08% are treated implicitly. On the other hand, for the exp–exp approach, three time-step subdomains arise, for which 6.60% of the elements of the model are analyzed with

Δ t_{1} = 8.60 {\cdot 10}^{- 5} s

, 93.20% with

Δ t_{2} = 1.72 {\cdot 10}^{- 4} s

, and 0.20% with

Δ t_{3} = 3.44 {\cdot 10}^{- 4} s

Computed time–history results for the axial displacement at the middle of the homogeneous rod are presented in Figure 22, for all discussed time-marching procedures. As one may observe in this figure, once again, the referred multi-domain methodologies provide considerably more accurate responses than standard techniques and, in this specific case, this improved accuracy becomes especially evident for the imp–exp formulation. In Table 3, a detailed description of the performance of each adopted solution technique is provided, further highlighting the enhanced effectiveness of the discussed hybrid approaches. As this table indicates, in this case, the exp–exp methodology provides the lowest CPU time for the analysis, closely followed by the imp–exp formulation, demonstrating once more that both these hybrid procedures are able to provide more accurate responses than standard techniques, considering several times faster computations than these approaches (in this specific case, computations more than 9 times faster are considered).

3.2. Applied Models

In this subsection, three geological subsurface applications are studied, further illustrating the robustness and efficiency of the discussed hybrid methodologies. These geological models are synthetically designed to emulate the level of complexity found in real subsurface structures, making them well-suited for assessing the applicability of the studied techniques in industrial contexts. The first model that is discussed herein describes an acoustic representation of the Buzios region (see Figure 23a) [28], situated in the southeastern part of the Santos Basin, Brazil, which is renowned for its pre-salt deposits. The second model, on the other hand, describes an elastodynamic application and it is commonly referred to as the 2DEW model (see Figure 23b) [29]. It replicates salt formations that are found in the Gulf of Mexico. Finally, the third model that is considered herein, known as the 2004 BP model (see Figure 23c) [30], offers a composite representation that amalgamates various geological sections from diverse regions, including the western and central/east Gulf of Mexico, as well as elements observed in the Caspian, North, and Trinidad seas. This elastodynamic model encompasses intricate stratigraphy, offering a comprehensive portrayal of subsurface characteristics. Indeed, all these three applications stand as realistic test cases, allowing us to properly evaluate the performances of the discussed techniques in practical geological scenarios (of course, 3D models considering distinct applications [33,34], could also be considered). Due to the associated high computational costs of these models, purely implicit methodologies are not employed for their analyses; thus, results are presented herein just considering the referred hybrid and standard explicit procedures.

The first model, which covers an area of 16.8 km × 8 km, is discretized considering a mesh composed of 4,864,312 linear triangular elements and it is illuminated by a pulse at its upper surface, at x = 8400 m. The second model, which represents an area of 35 km × 15 km, is discretized using a mesh composed of 717,139 linear triangular elements and, similarly to the first model, it is excited by a force applied at its upper surface, at x = 17,440 m. Finally, the third model, which extends across an area of 67.5 km × 46 km, is discretized using a mesh composed of 2,574,204 linear triangular elements and it is excited by a force applied at x = 33,750 m. In all these three applications, perfectly matched layers (PMLs) [31] are employed to avoid wave reflections at the vertical and lower horizontal boundaries of the discretized domain. The thicknesses of these PMLs are (i) 0.8 km, for model 1; (ii) 1 km, for model 2; and (iii) 1 km, for model 3.

Figure 24 illustrates the computed subdomains for these three models, considering the discussed hybrid methodologies. In this case, for the imp–exp approach, the following subdivisions are obtained, considering the percentages of elements in each subdomain: (i) model 1—71.46% explicit and 28.53% implicit, for

Δ t = 8.6751 \cdot 10^{- 4} s

; (ii) model 2—83.22% explicit and 16.78% implicit, for

Δ t = 3.4235 \cdot 10^{- 4} s

; and (iii) model 3—86.47% explicit and 13.53% implicit, for

Δ t = 4.2489 \cdot 10^{- 4} s

. For the exp–exp approach, on the other hand, the obtained subdivisions may be characterized as follows, considering the percentages of elements for each computed subdomain as well: (i) model 1 (3 subdomains)—51.88% for

Δ t_{1} = 6.405 {\cdot 10}^{- 4} s

, 22.66% for

Δ t_{2} = 1.281 {\cdot 10}^{- 3} s

, and 25.46% for

Δ t_{3} = 2.562 {\cdot 10}^{- 3} s

; (ii) model 2 (4 subdomains)—76.90% for

Δ t_{1} = 2.425 {\cdot 10}^{- 3} s

, 10.54% for

Δ t_{2} = 4.851 {\cdot 10}^{- 3} s

, 8.07% for

Δ t_{3} = 9.702 {\cdot 10}^{- 3} s

, and 4.49% for

Δ t_{4} = 1.940 {\cdot 10}^{- 2} s

; and (iii) model 3 (6 subdomains)—71.32% for

Δ t_{1} = 3.241 {\cdot 10}^{- 3} s

, 17.87% for

Δ t_{2} = 6.483 {\cdot 10}^{- 3} s

, 4.13% for

Δ t_{3} = 1.296 {\cdot 10}^{- 2} s

, 3.43% for

Δ t_{4} = 2.593 {\cdot 10}^{- 2} s

, 3.23% for

Δ t_{5} = 5.186 {\cdot 10}^{- 2} s

, and 0.02% for

Δ t_{6} = 1.037 {\cdot 10}^{- 1} s

In Table 4, a detailed description of the performance of each adopted solution technique is provided, for each model. As this table indicates, the exp–exp methodology stands as the most efficient approach for models 1 and 2, whereas, for model 3, the imp–exp technique becomes the most efficient procedure (and, for this third model, both hybrid techniques provide results twice faster than the considered standard formulations). In Figure 25, Figure 26 and Figure 27, screenshots are presented, depicting the computed scalar results for model 1 (Figure 25) and the modulus of the computed displacement results (visualized on a logarithmic scale, for enhanced clarity) for models 2 (Figure 26) and 3 (Figure 27). As these figures illustrate, the discussed hybrid methodologies provide similar results to those of the EG-α; however, as indicated in Table 4, they do so considering much lower CPU times. In fact, as illustrated in this third section, the discussed multi-domain hybrid approaches stand as very effective solution procedures, allowing complex wave propagation models to be very efficiently and accurately analyzed.

4. Conclusions

In this manuscript, adaptive implicit–explicit and explicit–explicit time-marching procedures are studied, discussing their performances and taking into account several wave propagation problems and configurations. In this context, theoretical and complex applied models have been regarded, allowing us to illustrate the good accuracy, efficiency, flexibility, and robustness of the discussed techniques. As it is reported in this paper, (i) both these multi-domain adaptive procedures are based only on single-step relations, standing as simple truly self-starting techniques; (ii) they are a function of the adopted temporal and spatial discretizations, as well as of the properties of the model and computed responses, establishing a link between these features, which allows us to enhance their related accuracies; (iii) both techniques consider self-adjusted time-step values, so that their performances may be improved; (iv) they also consider locally-defined time-integration parameters, further enhancing their effectiveness; (v) they enable advanced controllable algorithmic dissipation considering adaptive optimized calculations, efficiently eliminating the contributions of spatially unresolved high-frequency modes; and (vi) they stand as entirely automated formulations, requiring no decision, expertise, nor effort from the user.

The discussed implicit–explicit approach is extremely easy to implement and to apply, allowing us to straightforwardly unify the main positive aspects of both implicit and explicit analyses (such as guaranteed stability, low computational efforts, etc.). In fact, as illustrated in this paper, this approach may become extremely effective when just some small parts of the domain are “stiffer” than others (and the comparative better performance of the technique improves as greater differences in these different regions are observed), allowing relatively large time-step values to be considered without significantly increasing the solver effort of the hybrid approach. Since the discussed implicit–explicit formulation also computes and employs an optimized time-step value, this hybrid technique always becomes more efficient (or, at least, equally efficient) than its explicit or implicit counterparts, consistently providing a more attractive solution procedure. It may also be demonstrated that, whereas the non-dissipative explicit formulation considered in the discussed implicit–explicit technique is spectrally equivalent to the central difference method, its non-dissipative implicit formulation is more accurate than the trapezoidal rule, enabling us to compute very accurate responses. In addition, since this hybrid approach allows us to consider time-step values that are greater than those related to purely explicit analyses, its correspondent explicit elements become better discretized, allowing the errors of the considered time-marching procedure to better counterbalance the errors of the adopted spatial discretization formulation, further improving the accuracy of the analysis.

The discussed explicit–explicit approach, on the other hand, stands as a more intricate formulation to implement. However, it is also more flexible, allowing several multiple subdomains to be simultaneously considered (whereas the implicit–explicit approach just allows two basic subdomains to be simultaneously regarded, namely the implicit and the explicit subdomains). It also enables a more effective dissipative approach than the implicit–explicit formulation, as well as requiring less computational memory since its effective matrix is always diagonal. As the implicit–explicit technique, the reported explicit–explicit formulation automatically renders better space/time discretizations than purely explicit approaches, allowing space/time-related errors to be better counterbalanced and enhanced accuracy obtained. Since the studied explicit–explicit approach is highly flexible, it may become very effective once multiple domains are naturally physically represented in the analyzed model, allowing them to be very straightforwardly solved considering different time-step values. The main drawback of the discussed explicit–explicit formulation is, in fact, its complexity to be implemented, requiring subdomain interfaces to be mapped and multiple time interpolations to be carried out. However, the referred technique stands as an entirely automated formulation, requiring no extra burden from the user, as well as demanding no additional input data to be provided.

As one may observe, the discussed implicit–explicit and explicit–explicit formulations provide several positive characteristics, standing as very attractive time-marching procedures for wave propagation analyses. As discussed in this paper, although one of these two formulations may be more indicated than the other for some particular applications, according to the properties of the analyzed model, both these hybrid approaches stand as highly effective formulations, consistently providing better results than standard time-integration procedures.

Author Contributions

Conceptualization, D.S.J.; methodology, D.S.J.; software, D.S.J., I.d.S.S. and Pinto, L.R.P.; validation, D.S.J., I.d.S.S. and L.R.P.; formal analysis, D.S.J.; investigation, D.S.J., I.d.S.S. and L.R.P.; writing—original draft preparation, D.S.J., I.d.S.S. and L.R.P.; writing—review and editing, D.S.J., I.d.S.S., L.R.P. and W.J.M.; supervision, D.S.J. and W.J.M.; project administration, D.S.J. and W.J.M.; funding acquisition, D.S.J. and W.J.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), PRH-ANP (Programa de Recursos Humanos da Agencia Nacional do Petróleo, Gás Natural e Biocombustíveis), and PETROBRAS (CENPES 21066).

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Conflicts of Interest

The authors declare no conflict of interest.

References

Belytschko, T.; Mullen, R. Stability of explicit-implicit mesh partitions in time integration. Int. J. Numer. Methods Eng. 1978, 12, 1575–1586. [Google Scholar] [CrossRef]
Hughes, T.J.; Pister, K.S.; Taylor, R.L. Implicit-explicit finite elements in nonlinear transient analysis. Comput. Methods Appl. Mech. Eng. 1979, 17, 159–182. [Google Scholar] [CrossRef]
Soleimani, B.; Knoth, O.; Weiner, R. IMEX peer methods for fast-wave–slow-wave problems. Appl. Numer. Math. 2017, 118, 221–237. [Google Scholar] [CrossRef]
Miranda, I.; Ferencz, R.M.; Hughes, T.J. An improved implicit-explicit time integration method for structural dynamics. Earthq. Eng. Struct. Dyn. 1989, 18, 643–653. [Google Scholar] [CrossRef]
Soares, D., Jr. A model/solution-adaptive explicit-implicit time-marching technique for wave propagation analysis. Int. J. Numer. Methods Eng. 2019, 119, 590–617. [Google Scholar] [CrossRef]
Soares, D., Jr. An enhanced explicit–implicit time-marching formulation based on fully-adaptive time-integration parameters. Comput. Methods Appl. Mech. Eng. 2023, 403, 115711. [Google Scholar] [CrossRef]
Dujardin, G.; Lafitte, P. Asymptotic behaviour of splitting schemes involving time-subcycling techniques. IMA J. Numer. Anal. 2016, 36, 1804–1841. [Google Scholar] [CrossRef]
Belytschko, T.; Smolinski, P.; Liu, W.K. Stability of multi-time step partitioned integrators for first-order finite element systems. Comput. Methods Appl. Mech. Eng. 1985, 49, 281–297. [Google Scholar] [CrossRef]
Pinto, L.R.; Soares, D., Jr.; Mansur, W.J. Elastodynamic wave propagation modelling in geological structures considering fully-adaptive explicit time-marching procedures. Soil Dyn. Earthq. Eng. 2021, 150, 106962. [Google Scholar] [CrossRef]
Soares, D., Jr.; Pinto, L.R.; Mansur, W.J. A truly-explicit time-marching formulation for elastodynamic analyses considering locally-adaptive time-integration parameters and time-step values. Int. J. Solids Struct. 2023, 271, 112260. [Google Scholar] [CrossRef]
Tamma, K.K.; Zhou, X.; Sha, D. The time dimension: A theory towards the evolution, classification, characterization and design of computational algorithms for transient/dynamic applications. Arch. Comput. Methods Eng. 2000, 7, 67–290. [Google Scholar] [CrossRef]
Butcher, J.C. A history of Runge-Kutta methods. Appl. Numer. Math. 1996, 20, 247–260. [Google Scholar] [CrossRef]
Hulbert, G.M.; Chung, J. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. Comput. Methods Appl. Mech. Eng. 1996, 137, 175–188. [Google Scholar]
Noh, G.; Bathe, K.J. An explicit time integration scheme for the analysis of wave propagations. Comput. Struct. 2013, 129, 178–193. [Google Scholar] [CrossRef]
Kim, W. An accurate two-stage explicit time integration scheme for structural dynamics and various dynamic problems. Int. J. Numer. Methods Eng. 2019, 120, 1–28. [Google Scholar] [CrossRef]
Soares, D., Jr. An improved adaptive formulation for explicit analyses of wave propagation models considering locally-defined self-adjustable time-integration parameters. Comput. Methods Appl. Mech. Eng. 2022, 399, 115324. [Google Scholar] [CrossRef]
Newmark, N.M. A method of computation for structural dynamics. J. Eng. Mech. Div. 1959, 85, 67–94. [Google Scholar]
Chung, J.; Hulbert, G. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method. ASME. J. Appl. Mech. 1993, 60, 371–375. [Google Scholar] [CrossRef]
Bathe, K.J.; Baig, M.M.I. On a composite implicit time integration procedure for nonlinear dynamics. Comput. Struct. 2005, 83, 2513–2524. [Google Scholar]
Soares, D., Jr. A material/element-defined time integration procedure for dynamic analysis. Eng. Comput. 2023, 40, 1575–1601. [Google Scholar]
Gravouil, A.; Combescure, A. Multi-time-step explicit–implicit method for non-linear structural dynamics. Int. J. Numer. Methods Eng. 2001, 50, 199–225. [Google Scholar] [CrossRef]
Noels, L.; Stainier, L.; Ponthot, J.P. Combined implicit/explicit time-integration algorithms for the numerical simulation of sheet metal forming. J. Comput. Appl. Math. 2004, 168, 331–339. [Google Scholar] [CrossRef]
Curiel Sosa, J.L.; de Souza Neto, E.; Owen, D.R.J. A combined implicit–explicit algorithm in time for non-linear finite element analysis. Commun. Numer. Methods Eng. 2006, 22, 63–75. [Google Scholar] [CrossRef]
Ongaro, G.; Shojaei, A.; Mossaiby, F.; Hermann, A.; Cyron, C.J.; Trovalusci, P. Multi-adaptive spatial discretization of bond-based peridynamics. Int. J. Fract. 2023, 244, 1–24. [Google Scholar] [CrossRef]
Kennedy, J.; Eberhart, R. Particle swarm optimization. In Proceedings of the ICNN’95-International Conference on Neural Networks, Perth, WA, Australia, 27 November–1 December 1995; Volume 4, pp. 1942–1948. [Google Scholar]
Mansur, W.J. A Time-Stepping Technique to Solve Wave propagation Problems Using the Boundary Element Method. Ph.D. Thesis, University of Southampton, Southampton, UK, 1983. [Google Scholar]
Batra, R.C.; Porfiri, M.; Spinello, D. Free and forced vibrations of a segmented bar by a meshless local Petrov–Galerkin (MLPG) formulation. Comput. Mech. 2008, 41, 473–491. [Google Scholar] [CrossRef]
Freitas, D.F.; Siqueira, R.; Cetale, M.; González, M.; Soares, D.M. 3D Velocity Model Building of Búzios oil field—Improvements and updates. In Proceedings of the IX Simpósio Brasileiro de Geofísica, Curitiba, PR, Brazil, 4–6 October 2022. [Google Scholar]
Fehler, M. SEAM update: SEAM phase I-RPSEA update: Status of simulations. Lead. Edge 2012, 31, 1424–1426. [Google Scholar] [CrossRef]
Billette, F.J.; Brandsberg-Dahl, S. The 2004 BP velocity benchmark. In Proceedings of the 67th EAGE Conference & Exhibition, Madrid, Spain, 13–16 June 2005; p. cp-1. [Google Scholar]
Basu, U. Perfectly Matched Layers for Acoustic and Transient Waves. Dam Safety Research Program U.S. Department of the Interior Bureau of Reclamation: Washington, DC, USA, 2008. [Google Scholar]
Hughes, T.J. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis; Dover Publications Inc.: New York, NY, USA, 2000. [Google Scholar]
Van Pamel, A.; Sha, G.; Rokhlin, S.I.; Lowe, M.J.S. Finite-element modelling of elastic wave propagation and scattering within heterogeneous media. Proc. R. Soc. A Math. Phys. Eng. Sci. 2017, 473, 20160738. [Google Scholar] [CrossRef]
Imperiale, A.; Chatillon, S.; Darmon, M.; Leymarie, N.; Demaldent, E. UT simulation using a fully automated 3D hybrid model: Application to planar backwall breaking defects inspection. AIP Conf. Proc. 2018, 1949, 050004. [Google Scholar]

Figure 1. (a) Time interpolation and (b) computational flowchart for the sub-cycling process.

Figure 2. Spectral radii for the discussed solution procedure (Equation (4a,b)), considering the γ parameter defined by Equation (5b) (implicit approach) and the α parameter defined by (a) Equation (6b) and (b) Equation (6c), for

Ω_{e}^{m a x}

= 2.0, 2.1, …, 3.5 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Ω_{e}^{m a x}

= 2.0, 2.1, …, 3.5 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Figure 3. Spectral radii for the discussed solution procedure (Equation (4a,b)), considering the γ parameter defined by Equation (5a) (explicit approach) and the α parameter defined by (a) Equation (6b) and (b) Equation (6c), for

Ω_{e}^{m a x}

= 0.5, 0.6, …, 2.0 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Ω_{e}^{m a x}

= 0.5, 0.6, …, 2.0 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Figure 4. Period elongation and amplitude decay errors for the discussed solution procedure (Equation (4a,b)), considering the γ parameter defined by Equation (5b) (implicit approach) and the α parameter defined by (a) Equation (6b) and (b) Equation (6c), for

Ω_{e}^{m a x}

= 2.0, 2.1, …, 3.5 (lighter to darker gray color). Results for the CD and the TR are as well depicted as black dotted and dashed lines, respectively, for reference.

Ω_{e}^{m a x}

= 2.0, 2.1, …, 3.5 (lighter to darker gray color). Results for the CD and the TR are as well depicted as black dotted and dashed lines, respectively, for reference.

Figure 5. Period elongation and amplitude decay errors for the discussed solution procedure (Equation (4a,b)), considering the γ parameter defined by Equation (5a) (explicit approach) and the α parameter defined by (a) Equation (6b) and (b) Equation (6c), for

Ω_{e}^{m a x}

= 0.5, 0.6, …, 2.0 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Ω_{e}^{m a x}

= 0.5, 0.6, …, 2.0 (lighter to darker gray color). Results for the CD and the TR are also depicted as black dotted and dashed lines, respectively, for reference.

Figure 6. Adopted spatial discretizations for the first example: (a) discretization 1 (50 k elements); (b) discretization 2 (100 k elements); (c) discretization 3 (150 k elements); and (d) discretization 4 (200k elements).

Figure 7. Computed values for (1)

Ω_{e}^{m a x}

and (2)

γ_{e}^{n}

, for the imp–exp analyses, considering (a) discretization 1; (b) discretization 2; (c) discretization 3; and (d) discretization 4.

Figure 7. Computed values for (1)

Ω_{e}^{m a x}

and (2)

γ_{e}^{n}

, for the imp–exp analyses, considering (a) discretization 1; (b) discretization 2; (c) discretization 3; and (d) discretization 4.

Figure 8. Computed values for (1)

Δ t_{e}

and (2)

Δ t_{i}

for the exp–exp analyses, considering (a) discretization 1; (b) discretization 2; (c) discretization 3; and (d) discretization 4.

Figure 8. Computed values for (1)

Δ t_{e}

and (2)

Δ t_{i}

for the exp–exp analyses, considering (a) discretization 1; (b) discretization 2; (c) discretization 3; and (d) discretization 4.

Figure 9. Time history results for

u

, at a point located 10 m horizontally away from the applied source (discretization 4), considering solutions by (a) implicit and (b) explicit methods, as well as their hybrid extensions.

Figure 9. Time history results for

u

, at a point located 10 m horizontally away from the applied source (discretization 4), considering solutions by (a) implicit and (b) explicit methods, as well as their hybrid extensions.

Figure 10. Convergence curves for the discussed time-marching procedures and discretizations.

Figure 11. Time–history results for the axial displacement at the middle of the rod, considering

c_{2} / c_{1} = 4

: (a) implicit and (b) explicit approaches, as well as their hybrid extensions.

Figure 11. Time–history results for the axial displacement at the middle of the rod, considering

c_{2} / c_{1} = 4

: (a) implicit and (b) explicit approaches, as well as their hybrid extensions.

Figure 12. Time–history results for the axial displacement at the middle of the rod, considering

c_{2} / c_{1} = 6

: (a) implicit and (b) explicit approaches, as well as their hybrid extensions.

Figure 12. Time–history results for the axial displacement at the middle of the rod, considering

c_{2} / c_{1} = 6

: (a) implicit and (b) explicit approaches, as well as their hybrid extensions.

Figure 13. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 2

Δ t

value) approaches;

c_{2} / c_{1} = 2

Figure 14. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal Δt value) approaches;

c_{2} / c_{1} = 3

c_{2} / c_{1} = 3

Figure 15. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 4

Δ t

value) approaches;

c_{2} / c_{1} = 4

Figure 16. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 5

Δ t

value) approaches;

c_{2} / c_{1} = 5

Figure 17. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 6

Δ t

value) approaches;

c_{2} / c_{1} = 6

Figure 18. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 7

Δ t

value) approaches;

c_{2} / c_{1} = 7

Figure 19. Computed errors and CPU times for different material distributions, considering the exp–exp (open circle) and the imp–exp (solid circle, lighter gray color is depicted when a purely exp solution takes place, following the computed optimal

Δ t

value) approaches;

c_{2} / c_{1} = 8

Δ t

value) approaches;

c_{2} / c_{1} = 8

Figure 20. Analytical (black), imp–exp (light purple), and exp–exp (dark purple) time–history responses (all curves are visually the same) for the axial displacements at the middle of the rod, for

c_{2} / c_{1}

equals to (a) 2, (b) 3, (c) 4, (d) 5, (e) 6, (f) 7, and (g) 8; considering a percentage of material 2 equals to (1) 10%, (2) 50%, and (3) 90%.

c_{2} / c_{1}

equals to (a) 2, (b) 3, (c) 4, (d) 5, (e) 6, (f) 7, and (g) 8; considering a percentage of material 2 equals to (1) 10%, (2) 50%, and (3) 90%.

Figure 21. (a) Adopted spatial discretization for the homogeneous rod and its computed values for: (b)

Ω_{e}^{m a x};

(c)

γ_{e}^{n}

; (d)

Δ t_{e}

; and (e)

Δ t_{i}

Figure 21. (a) Adopted spatial discretization for the homogeneous rod and its computed values for: (b)

Ω_{e}^{m a x};

(c)

γ_{e}^{n}

; (d)

Δ t_{e}

; and (e)

Δ t_{i}

Figure 22. Time–history results for the axial displacement at the middle of the homogeneous rod, considering (a) implicit and (b) explicit approaches, as well as their hybrid extensions.

Figure 23. Geological models: (a) model 1—Buzios; (b) model 2—2DEW; and (c) model 3—2004BP.

Figure 24. Subdomain divisions for the (1) imp–exp (

γ_{e}^{n}

values are depicted) and (2) exp–exp (

Δ t_{i}

values are depicted) methods, for (a) model 1; (b) model 2; and (c) model 3.

Figure 24. Subdomain divisions for the (1) imp–exp (

γ_{e}^{n}

values are depicted) and (2) exp–exp (

Δ t_{i}

values are depicted) methods, for (a) model 1; (b) model 2; and (c) model 3.

Figure 25. Computed results along the discretized domain of model 1, for the (a) EG-α, (b) imp–exp, and (c) exp–exp methods at different time instants: (1) 3 s and (2) 6 s.

Figure 26. Computed results along the discretized domain of model 2, for the (a) EG-α, (b) imp–exp, and (c) exp–exp methods at different time instants: (1) 3 s and (2) 6 s.

Figure 27. Computed results along the discretized domain of model 3, for the (a) EG-α, (b) imp–exp, and (c) exp–exp methods at different time instants: (1) 10 s; (2) 15 s; and (3) 25 s.

Table 1. Performance of the methods for the first example *.

Discretization	Method	Δt (10⁻² s)	Error (10⁻¹)	CPU Time (s)
1	TR	14.9256 (3.60)	7.78 (2.08)	19.6 (2.27)
	IG-α	14.9256 (3.60)	7.15 (1.91)	19.8 (2.29)
	IC	14.9256 (3.60)	7.31 (1.96)	20.6 (2.38)
	imp	14.9256 (3.60)	7.51 (2.01)	19.5 (2.25)
	CD	4.60382 (1.11)	7.98 (2.13)	14.5 (1.67)
	EG-α	4.14927 (1.00)	7.85 (2.10)	15.9 (1.85)
	EC	8.62073 (2.08)	7.91 (2.12)	16.1 (1.86)
	exp	4.60382 (1.11)	5.24 (1.40)	14.0 (1.61)
	imp–exp	14.9256 (3.60)	4.20 (1.12)	9.45 (1.09)
	exp–exp	36.8306 (8.88)	3.73 (1.00)	8.67 (1.00)
2	TR	10.0516 (4.77)	7.22 (2.76)	32.7 (3.02)
	IG-α	10.0516 (4.77)	6.68 (2.25)	32.6 (3.01)
	IC	10.0516 (4.77)	6.81 (2.29)	34.2 (3.17)
	imp	10.0516 (4.77)	7.01 (2.36)	32.8 (3.03)
	CD	2.33704 (1.11)	7.62 (2.57)	22.9 (2.12)
	EG-α	2.10629 (1.00)	7.39 (2.49)	24.2 (2.42)
	EC	4.37614 (2.08)	7.51 (2.53)	28.9 (2.67)
	exp	2.33704 (1.11)	5.09 (1.71)	22.7 (2.09)
	imp–exp	10.0516 (4.77)	3.46 (1.17)	11.2 (1.04)
	exp–exp	18.6963 (8.88)	2.97 (1.00)	10.8 (1.00)
3	TR	7.59321 (6.74)	6.53 (3.08)	101.4 (4.71)
	IG-α	7.59321 (6.74)	6.02 (2.84)	101.9 (4.74)
	IC	7.59321 (6.74)	6.12 (2.88)	118.9 (5.53)
	imp	7.59321 (6.74)	6.32 (2.98)	103.6 (4.81)
	CD	1.24956 (1.11)	7.03 (3.32)	50.0 (2.33)
	EG-α	1.12619 (1.00)	6.73 (3.17)	57.4 (2.67)
	EC	2.33983 (2.08)	6.92 (3.27)	64.8 (3.01)
	exp	1.24956 (1.11)	4.57 (2.15)	50.2 (2.34)
	imp–exp	7.59321 (6.74)	2.60 (1.23)	21.5 (1.00)
	exp–exp	19.9930 (17.8)	2.12 (1.00)	26.4 (1.23)
4	TR	6.20276 (8.48)	6.17 (3.55)	159.9 (3.85)
	IG-α	6.20276 (8.48)	5.60 (3.22)	160.2 (3.86)
	IC	6.20276 (8.48)	5.71 (3.29)	174.1 (4.19)
	imp	6.20276 (8.48)	5.93 (3.41)	160.2 (3.14)
	CD	0.81175 (1.11)	6.66 (3.84)	88.8 (2.14)
	EG-α	0.73160 (1.00)	6.38 (3.67)	114.9 (2.77)
	EC	1.52001 (2.08)	6.59 (3.79)	126.7 (3.05)
	exp	0.81174 (1.11)	4.18 (3.40)	89.0 (2.15)
	imp–exp	6.20276 (8.48)	2.19 (1.26)	41.5 (1.00)
	exp–exp	12.9880 (17.8)	1.74 (1.00)	49.5 (1.19)

* Relative values are provided in parenthesis.

Table 2. Performance of the methods for the second example * (

a = L / 2

Table 2. Performance of the methods for the second example * (

a = L / 2

$c_{2} / c_{1}$	Method	$Δ t (10^{- 4} s)$	$Error (10^{- 3})$	CPU Time (s)
4	TR	0.62499 (1.11)	7.96 (9.68)	121.4 (8.70)
	IG-α	0.62499 (1.11)	7.56 (9.19)	120.1 (8.61)
	IC	0.62499 (1.11)	8.36 (10.1)	230.4 (16.5)
	imp	0.62499 (1.11)	5.77 (7.01)	120.6 (8.64)
	CD	0.62499 (1.11)	5.84 (7.10)	16.6 (1.19)
	EG-α	0.56317 (1.00)	6.20 (7.53)	18.3 (1.31)
	EC	1.17046 (2.08)	6.24 (7.58)	22.4 (1.60)
	exp	0.62499 (1.11)	5.75 (6.98)	16.6 (1.19)
	imp–exp	0.62499 (1.11)	5.75 (6.98)	16.6 (1.19)
	exp–exp	2.49998 (4.44)	0.82 (1.00)	13.9 (1.00)
6	TR	2.49998 (6.66)	18.87 (3.94)	30.43 (1.84)
	IG-α	2.49998 (6.66)	22.23 (4.64)	31.52 (1.91)
	IC	2.49998 (6.66)	18.99 (3.96)	56.68 (3.43)
	imp	2.49998 (6.66)	18.50 (3.86)	30.35 (1.83)
	CD	0.41666 (1.11)	8.49 (1.77)	20.65 (1.25)
	EG-α	0.37545 (1.00)	8.70 (1.82)	25.29 (1.53)
	EC	0.78030 (2.08)	8.78 (1.83)	28.25 (1.71)
	exp	0.41666 (1.11)	8.41 (1.76)	20.78 (1.26)
	imp–exp	2.49998 (6.66)	8.93 (1.86)	16.56 (1.00)
	exp–exp	1.66664 (4.44)	4.79 (1.00)	16.54 (1.00)

* Relative values are provided in parenthesis.

Table 3. Performance of the methods for the homogeneous rod *.

Method	$Δ t (10^{- 4} s)$	$Error (10^{- 3})$	CPU Time (s)
TR	1.42924 (1.84)	8.65 (1.67)	28.8 (4.97)
IG-α	1.42924 (1.84)	8.34 (1.61)	29.9 (5.15)
IC	1.42924 (1.84)	8.49 (1.63)	53.1 (9.15)
imp	1.42924 (1.84)	8.29 (1.60)	28.6 (4.92)
CD	0.86030 (1.11)	7.42 (1.43)	8.2 (1.41)
EG-α	0.77536 (1.00)	8.01 (1.54)	9.4 (1.61)
EC	1.61092 (2.08)	7.43 (1.43)	11.0 (1.89)
exp	0.86090 (1.11)	6.41 (1.23)	8.3 (1.42)
imp–exp	1.42924 (1.84)	5.18 (1.00)	6.4 (1.11)
exp–exp	3.44120 (4.44)	6.19 (1.19)	5.8 (1.00)

* Relative values are provided in parenthesis.

Table 4. Performance of the methods for the applied models *.

Model	Method	$Δ t (10^{- 3} s)$	CPU Time (s)
1	CD	0.64051 (1.11)	5233.7 (1.68)
	EG-α	0.57727 (1.00)	5462.8 (1.75)
	EC	1.19936 (2.08)	5683.6 (1.83)
	exp	0.64051 (1.11)	5261.2 (1.69)
	imp–exp	0.86751 (1.50)	4282.5 (1.38)
	exp–exp	2.56204 (4.44)	3112.9 (1.00)
2	CD	0.24255 (1.11)	8661.5 (1.99)
	EG-α	0.21861 (1.00)	8991.6 (2.06)
	EC	0.45419 (2.08)	9511.1 (2.18)
	exp	0.24255 (1.11)	8654.5 (1.98)
	imp–exp	0.34235 (1.57)	5436.3 (1.25)
	exp–exp	1.94044 (8.88)	4362.2 (1.00)
3	CD	0.32415 (1.11)	11322.8 (2.15)
	EG-α	0.29215 (1.00)	12283.3 (2.34)
	EC	0.60698 (2.08)	14773.1 (2.81)
	exp	0.32415 (1.11)	11737.7 (2.23)
	imp–exp	0.42489 (1.45)	5259.40 (1.00)
	exp–exp	10.3729 (35.5)	5883.30 (1.11)

* Relative values are provided in parenthesis.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Soares, D., Jr.; Sales, I.d.S.; Pinto, L.R.; Mansur, W.J. A Study on Adaptive Implicit–Explicit and Explicit–Explicit Time Integration Procedures for Wave Propagation Analyses. Acoustics 2024, 6, 651-680. https://doi.org/10.3390/acoustics6030036

AMA Style

Soares D Jr., Sales IdS, Pinto LR, Mansur WJ. A Study on Adaptive Implicit–Explicit and Explicit–Explicit Time Integration Procedures for Wave Propagation Analyses. Acoustics. 2024; 6(3):651-680. https://doi.org/10.3390/acoustics6030036

Chicago/Turabian Style

Soares, Delfim, Jr., Isabelle de Souza Sales, Lucas Ruffo Pinto, and Webe João Mansur. 2024. "A Study on Adaptive Implicit–Explicit and Explicit–Explicit Time Integration Procedures for Wave Propagation Analyses" Acoustics 6, no. 3: 651-680. https://doi.org/10.3390/acoustics6030036

Article Menu

A Study on Adaptive Implicit–Explicit and Explicit–Explicit Time Integration Procedures for Wave Propagation Analyses

Abstract

1. Introduction

2. Methods

2.1. Implicit–Explicit Approach

2.2. Explicit–Explicit Approach

2.3. Properties of the Methods

3. Results and Discussions

3.1. Theoretical Models

3.2. Applied Models

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI