#### Abstract

Viscoelastic artificial boundary elements are one of the most commonly used artificial boundaries when solving dynamic soil-structure interactions or near-field wave propagation problems. However, due to the lack of clear and practical stability criteria for the explicit algorithm that considers the influence of viscoelastic artificial boundary elements, the determination of the stable time increment in such numerical analyses is still a challenge. In this study, we proposed a numerical stability analysis method for the explicit algorithm with a 3D viscoelastic artificial boundary element based on the idea of a subsystem. Through this method, the artificial boundary subsystem that controls the stability of the overall numerical system is determined, and the analytical solution for the stability condition of the explicit integration algorithm with 3D viscoelastic artificial boundary elements is obtained. On this basis, the maximum time increment for solving dynamic problems with viscoelastic artificial boundary elements can be determined.

#### 1. Introduction

During the dynamic analysis of a soil-structure interaction (SSI) system or the numerical simulation of a near-field wave propagation problem, a finite domain is usually intercepted from a semi-infinite medium to establish the near-field model, and the cutoff boundaries should be manipulated to simulate the wave radiation effect. Currently, artificial boundary (also known as an absorbing boundary or a nonreflecting boundary) techniques are the most commonly used methods for wave radiation problems. Typical artificial boundary techniques include the boundary element method (BEM) [1, 2], the perfectly matched layers (PML) [3–5], the transmitting boundary [6, 7], the viscous boundary [8], and the viscoelastic boundary [9–11] techniques. Among these methods, the viscoelastic boundary has received much attention because of the advantages of high accuracy and the spatial decoupling characteristics. However, the preprocessing operations of the viscoelastic boundary in the general finite element programs are still slightly complicated, since the equivalent spring-damper systems should be applied on each boundary node of the near-field model. To solve this problem, the viscoelastic artificial boundary element technique is developed [12–14], which is accomplished by setting a layer of elements with specified material properties along the normal direction of the cutoff boundaries to work identically to the traditional viscoelastic boundary. Through the viscoelastic artificial boundary element technique, the precision of the original viscoelastic boundary can be maintained, and the preprocessing process can be further simplified; thus, it has gained widespread use in recent years [15–19].

However, when the explicit integration method is adopted for a calculation, the stiffness, viscous damping, and geometrical dimensions of the viscoelastic artificial boundary elements will change the numerical stability of the original computational system. At present, due to the lack of clear and practical stability criteria of the explicit integration algorithm considering the influence of viscoelastic artificial boundary elements, the judgment and selection of the stable time increment in such a numerical analysis is still a challenge. Therefore, when viscoelastic artificial boundary elements are adopted in dynamic analysis, unconditionally stable implicit integration algorithms are usually employed to avoid the possible numerical instability caused by viscoelastic artificial boundary elements [20, 21]. Due to the iterative solving process in the implicit algorithm, the computational burden is heavy, and the efficiency is poor, especially for large-scale three-dimensional (3D) dynamic SSI problems. Therefore, it is necessary to study the stability of the explicit integration algorithm in numerical systems containing viscoelastic artificial boundary elements.

A conventional method for the stability analysis of a numerical system including the artificial boundaries is to establish the transfer matrix of the overall system according to the time-domain integration algorithm, and then the stability condition can be obtained through the eigenvalue analysis of the transfer matrix [22]. However, it is usually difficult to derive the analytical eigenvalue solutions of the large-scale matrix, and obtaining the eigenvalues through numerical methods is a very time-consuming process. Therefore, it is difficult to determine the maximum stable time increment for the explicit calculation through this method. Guan and Liao [23] as well as Liao and Liu [24] proposed a stability analysis method based on a subsystem containing artificial boundaries to study the numerical stability of transmitting boundaries. The subsystems adopted in this method are composed of several rows and lines of nodes along the tangential direction and the normal direction of the transmitting boundary. However, the numbers of degrees of freedom in such subsystems are still relatively large; thus, the stability conditions of the artificial boundaries can only be analyzed through numerical methods. The practical analytical form of the stability criterion, however, has not been given.

In this study, we propose a numerical stability analysis method for explicit integration algorithms based on the idea of subsystems. Then, three typical 3D subsystems containing the artificial boundary elements are established, and the spectral radiuses of the transfer matrixes of those subsystems are analyzed to determine the local stability condition. The artificial boundary subsystem that controls the stability of the overall system is determined through a comparison analysis, and the analytical solution for the stability condition of the explicit stepwise integration algorithm considering the 3D viscoelastic artificial boundary elements is proposed.

#### 2. The Equivalent Physical Parameters of the 3D Viscoelastic Artificial Boundary Elements

The viscoelastic artificial boundary elements are made up of a layer of equivalent elements extending outward along the normal direction at the cutoff boundaries of the near-field model, as shown in Figure 1. The outermost nodes of the artificial boundary elements are fixed, and the material parameters are designated to simulate the wave absorption effect of the viscoelastic artificial boundaries.

We assume that the hexahedron element has a length of 2*a*, a width of 2*b*, and a height of *h*. The expression of the element matrix can be derived from the principle of virtual work [14]:where is the stiffness matrix, is the elastic modulus matrix, and is the geometric matrix of the equivalent element. By fixing nodes 5–8 of the element, the stiffness matrix of the equivalent element can be obtained:wherewhere is the elastic modulus, is the shear modulus, and is the Poisson ratio of the artificial boundary element. When *h* is much smaller than *a* and *b*, the stiffness matrix of the equivalent artificial boundary element is approximately

The element stiffness matrix formed by the 3D viscoelastic artificial boundary can be expressed as [11, 14]wherewhere and are the normal and tangential stiffness of the spring, respectively; *R* is the distance from the wave source to the artificial boundary node; *G* is the shear modulus of the internal medium; and and are the tangential and normal viscoelastic artificial boundary parameters, whose recommended values are 0.67 and 1.33 for 3D conditions, respectively. To make the equivalent element have the same rigidity as the viscoelastic artificial boundary, the following formula should be established:

Then, the equivalent elastic modulus and Poisson ratio of the viscoelastic artificial boundary elements can be deduced:

The damping matrix is proportional to the stiffness matrix, which is given as follows [14]:

Then, the damping coefficient can be derived:where , , and are the mass density, shear wave velocity, and compression wave velocity, respectively, of the internal medium. The mass density of the viscoelastic artificial boundary element is zero.

Although the assumption that the thickness of the boundary element layer *h* is much smaller than the width *L* has been introduced during the derivation of the equivalent stiffness matrix of the viscoelastic artificial boundary element, further studies [13, 14] have shown that the thickness of the boundary element *h* is quite robust and can be flexibly selected with little impact on the accuracy. Since the stability of the explicit algorithm is usually controlled by the minimum element size of the model, in practical applications, it is recommended to set the thickness of the boundary element *h* to be consistent with the internal element size, namely, *h* = *L*, to exclude the effect of the boundary element size on the overall stability, as shown in Figure 2.

#### 3. Stability Analysis Method for the Explicit Integration Algorithm Based on the Idea of a Subsystem

The main purpose of the stability analysis in the explicit integration algorithm is to obtain the discretized time increment , which satisfies the stability requirement in the calculation. For a clear theoretical derivation, in the stability analysis, a homogeneous medium and a uniform discretization mesh are usually employed, as shown in Figure 3, to obtain the general stability criterion. For an inhomogeneous model with nonuniform mesh generation, the stable time increment of the entire computational system can be comprehensively determined according to the smallest mesh size and the material properties.

From the previous studies on the numerical stability of integration algorithms, a common conclusion has been reached that the numerical stability is closely related to the cutoff frequency of the computational system [25] (the highest-order natural frequency of the model). Considering that the mode corresponding to the cutoff frequency is usually presented as the localized staggered vibration of adjacent nodes, it can be deduced that the stability of the overall system can be investigated through the numerical stability analysis of a local subsystem model.

In the following sections, two examples are adopted to illustrate the above deduction. First, we focus on the one-dimensional (1D) wave propagation problem. The discrete model of a 1D shear beam and the highest-order mode shape are shown in Figure 4, where *L* is the length of the beam element; the red circles represent the finite element nodes; the dashed lines represent the highest-order model; and the yellow squares represent the nodes of the vibration mode, which remain stationary during vibration. In addition, the shear modulus, the mass density, and the cross-sectional area of the beam model are designated as *G*, *ρ*, and *A*, respectively.

Then, two subsystems are intercepted from the shear beam model for analysis. Subsystem A is composed of two half elements between two adjacent nodes of the vibration mode with fixed constraints imposed on both ends, as shown in Figure 5(a). Since the nodes of the vibration mode are fixed, the vibration mode of subsystem A is identical to the localized shape of the highest-order mode of the overall beam model. Therefore, the natural frequency of subsystem A is consistent with the highest-order natural frequency of the entire shear beam. Subsystem B is composed of two beam elements, with fixed constraints on the two endpoints. Therefore, there is only one degree of freedom in subsystem B, as shown in Figure 5(b). Because of the different vibration modes between subsystem B and the overall beam model, the natural frequency of subsystem B is different from the cutoff frequency of the original system.

**(a)**

**(b)**

Since there is a classical solution for the stability condition in a 1D undamped linear elastic system, to compare the stability conditions of the subsystems with those of the overall model, the damping of the subsystem is not considered here. According to the concept of the shear beam, the shear stiffness , mass , and natural frequency of the 1D subsystem A can be obtained as follows:where *C*_{S} is the velocity of the shear wave and .

In addition, according to the shear stiffness and mass of subsystem A, the single degree of freedom motion equation can be established. When the central difference method is adopted for the time-domain stepwise integration calculation, the following condition should be met to ensure the stability of the explicit calculation:

By substituting (11) into (12), the numerical stability condition of subsystem A can be obtained as follows:

Equation (13) is the exact stability condition of the central difference algorithm in the 1D wave propagation problem, indicating that the stability condition of the numerical integration algorithm in a 1D continuous medium can be obtained through the stability analysis of an A-type subsystem.

Through similar analysis procedures, the shear stiffness , mass , natural frequency , and stability condition in the center differential method of subsystem B can be expressed as follows:

A comparison between (13) and (14) shows that the stability condition of subsystem B is loose compared to that of subsystem A, which means that the stability condition of subsystem B provides an upper limit estimation of the maximum stable time increment of the entire calculation system.

Then, the previous deduction in the two-dimensional (2D) condition is verified. The 2D plane model and the corresponding highest-order mode are shown in Figure 6, in which the 4-node isoparametric element with a side length of *L* is adopted for model discretization. The mass density and the elastic modulus are designated as and *E*. There is a classical solution of the stability condition in the 2D undamped linear elastic system, namely, [25], where is the velocity of the compression wave. The premise of the above classical solution is that the Poisson ratio is 0, leading to . To compare the numerical stability conditions of the subsystems with the classical solution, the Poisson ratio of the system is set to zero. Similarly, the red circles represent the finite element nodes, and the yellow squares represent the nodes of the vibration mode.

Two 2D subsystems are intercepted from the original computational model, as shown in Figure 7. Subsystem A is composed of four quarters of an element. The four corner nodes of subsystem A (nodes 6–9) are the nodes of the vibration mode and remain stationary during the analysis. The normal constraints are applied on the four side nodes (nodes 2–5), which can accurately simulate the highest-order vibration mode shape. The above boundary conditions of subsystem A can be expressed as follows:where *u* is the nodal displacement, the subscripts *x* and *y* denote the *x*- and *y*-directions, respectively, and the subscripts 1–6 denote the node numbers in Figure 7.

**(a)**

**(b)**

The motion equation for the 2D subsystem of subsystem A can be expressed as follows:

Then, the natural frequency and the stability condition in the central differential method of the 2D subsystem of subsystem A can be calculated as

Equation (17) is the exact stability condition of the central difference algorithm for 2D wave problems, indicating that the stability condition of the numerical integration algorithm in 2D continuous medium can be obtained through the stability analysis of the 2D subsystem of subsystem A.

The other subsystem shown in Figure 7(b) (subsystem B) consists of four adjacent elements, with fixed constraints imposed on all the outer nodes (nodes 2–9). Then, the motion equation of subsystem B can be expressed as

Then, the natural frequency and the stability condition in the central differential method of subsystem B can be calculated as

It can be seen from the comparison between (17) and (19) that, for 2D conditions, the maximum stable incremental time step of subsystem B also provides an upper limit estimation of the numerical stability condition of the entire calculation system.

The results of the 1D and 2D subsystem analyses show that (1) if the cutoff frequency of the entire system and the corresponding vibration mode can be accurately determined, the partial subsystem can be intercepted from the entire system according to the characteristics of the highest-order mode with specified constraints applied on the boundary nodes of the subsystem to simulate the shape of the highest-order vibration mode. Then, the maximum stable time increment of the entire system can be obtained through the stability analysis of the subsystem. (2) If it is difficult to determine the highest-order natural frequency and the corresponding mode of the entire system, a subsystem composed of all the elements connected to a single node can be intercepted to perform the stability analysis, whose maximum time increment provides an upper approximation of the stability condition of the overall system.

When the viscoelastic artificial boundary elements are applied in the computational system, the numerically stable condition of the internal domain can be obtained. The cutoff frequency in the artificial boundary region and the corresponding vibration mode, however, is difficult to determine. In this situation, a B-type subsystem can be employed for the numerical stability analysis.

Furthermore, for the stability analysis of an inhomogeneous subsystem, the motion equation can be written in the following form according to the explicit stepwise integration algorithm:rwhere the vector **U** consists of the nodal acceleration, velocity, and displacement of the subsystem; **Q** is the external force vector; **A** is the transfer matrix, which is related to the discrete steps in the space domain and time domain, as well as the material parameters of the subsystem; and the superscript *p* denotes the natural number, .

The following conditions should be met to ensure the stability of the numerical integration algorithm [26]:(1)If transfer matrix **A** has no repeated eigenvalues, then (), where is the spectral radius and denotes the *i*th eigenvalue of transfer matrix **A**(2)If transfer matrix **A** has repeated eigenvalues, the moduli of the repeated eigenvalues should be less than 1

Therefore, the numerical stability analysis of the subsystem is based on the generation of the transfer matrix and the calculation of the spectral radius.

#### 4. Stability Analysis of 3D Artificial Boundary Subsystems

When the 3D viscoelastic artificial boundary elements are employed in the numerical models, the B-type subsystem can be classified into three types according to the interception region, as shown in Figure 8: (a) the surface boundary subsystem, which is intercepted from the lateral or bottom surface of the model and consists of four artificial boundary elements and four internal elements; (b) the edge boundary subsystem, which is intercepted from the lateral or bottom edge regions and consists of six artificial boundary elements and two internal elements; and (c) the corner boundary subsystem is intercepted from the corner region and consists of seven artificial boundary elements and one internal element. The meshes in the model are uniformly generated, and the side length is *L*.

The boundary condition of these three kinds of subsystems is obtained by fixing the 26 surface nodes of the system and letting the remaining central node be free. The subsystems with all fixed surface nodes may have a higher natural vibration frequency, which can ensure that the numerical stability conditions of the subsystem are as close as possible to the overall system. In addition, there is only one free node in such a subsystem, which is helpful for deriving the analytical solution of the stability condition.

For the internal medium, the elastic modulus, shear modulus, mass density, and Poisson ratio are designated *E*, *G*, , and . Damping is not considered in the internal region. Then, the stability condition of the internal domain can be calculated by [25]

For the viscoelastic artificial boundary elements, the elastic modulus , the damping coefficient , the Poisson ratio , and the mass density can be determined according to (8) and (10).

Then, each type of subsystem is focused on, and the corresponding stability analysis is conducted in the following sections.

##### 4.1. Stability Analysis of the Surface Boundary Subsystem

Consider a regular hexahedral element from the surface boundary subsystem (Figure 8(a)), as shown in Figure 9.

The concentrated mass matrix of the regular hexahedral element of the internal medium can be expressed as

The stiffness matrix of the regular hexahedral element of the internal medium can be expressed aswhere

Node 1: ; Node 2: ; Node 3: ; Node 4: ; Node 5: ; Node 6: ; Node 7: ; and Node 8: .

The internal medium element contains the mass matrix and the stiffness matrix, while the artificial boundary element contains the damping matrix in addition to the mass matrix and stiffness matrix. The damping matrix can be obtained by (9) and (10). According to this method, the mass, stiffness, and damping matrix of a single internal medium element and artificial boundary element can be obtained. Then, the eight elements in the subsystem can be assembled according to the node number (the node number is shown in Figure 9). The fixed boundary conditions are applied to the surface nodes of the subsystem. Then, the motion equation of the central node can be expressed as follows:where the subscripts *x*, *y*, and *z* denote the *x*-, *y*-, and *z*-directions, respectively.

It can be seen from (25) that the motion equations in the *x*-, *y*-, and *z*-directions are decoupled and identical to each other. Therefore, the stability analysis can be performed in only one direction. Taking the *x*-direction as an example, the explicit time-domain stepwise integration format based on the central difference method can be expressed as follows [27]:

By combining (25) and (26), the following expression can be obtained:where

The eigenvalues of transfer matrix **A** are

By substituting (29) into the stability condition , the following expression can be obtained:

Equation (30) provides the analytical solution for the stability condition of the surface boundary subsystem. The stability condition is related to the wave velocities *C*_{P} and *C*_{S}, the mass density *ρ*, the Poisson ratio *μ* of the internal media, and the element size *L*.

##### 4.2. Stability Analysis of the Edge Boundary Subsystem

For the edge boundary subsystem, following the same procedure of matrix assembly, the motion equation of the central node can be expressed as follows:

It can be seen from (31) that the motion equations in different directions are coupled with each other, and the stability analysis of the overall subsystem is required. By combining (26) and (31), the following expression can be obtained:where

The maximum eigenvalue of the transfer matrix **A** can be solved as follows:

By substituting (34) into the stability condition , the analytical solution of the stability condition of the edge boundary subsystem can be obtained:

##### 4.3. Stability Analysis of the Corner Boundary Subsystem

For the corner boundary subsystem, following the same procedure matrix assembly, the motion equation of the central node can be expressed as follows:

The motion equations of the corner boundary subsystem in different directions are also coupled with each other. By combining (26) and (36), the following expression can be obtained:where

The maximum eigenvalue of the transfer matrix **A** can be solved as follows:

By substituting (39) into the stability condition , the analytical solution of the stability condition of the corner boundary subsystem can be obtained:

#### 5. Validation

A 3D homogeneous elastic half space model is established, as shown in Figure 10. The dimensions of the model are 100 m × 100 m × 50 m, and the mesh size is 2 m × 2 m × 2 m. The mass density of the internal medium *ρ* is 2500 kg/m^{3}, the shear wave velocity *C*_{S} is 500 m/s, the compression wave velocity *C*_{P} is 935 m/s, and the Poisson ratio is 0.3. The viscoelastic artificial boundary elements are applied on the cutoff boundaries to absorb the outgoing waves, and their material parameters are calculated through (8) and (10). A concentrated impulsive force is imposed on the geometric center of the model along the *z*-direction, and its time history is shown in Figure 11.

According to (21), (30), (35), and (40), the maximum stable time increments for the internal region , the surface boundary subsystem , the edge boundary subsystem , and the corner boundary subsystem in this model can be calculated, whose values are 0.0021 s, 0.00176 s, 0.00071 s, and 0.00032 s, respectively. The central difference method is adopted, and the accuracy of the above stability analysis can be verified by specifying different fixed time increments and conducting explicit dynamic analyses.

It is important to note that the minimum distances from the geometric center of the model to the surfaces *L*_{S}, to the edges *L*_{E}, and to the corners *L*_{C} gradually increase (*L*_{S} < *L*_{E} < *L*_{C}); thus, the time for the wave crest arriving at the surfaces, the edges, and the corners also increases one by one.

First, the maximum stable time increment of the internal domain is adopted for the explicit analysis ( = 0.0021 s). When the pulse wave propagates within the internal region, the calculation can be executed accurately. However, once the waves reach the bottom surface of the model, the stability condition of the surface boundary subsystem is no longer satisfied (); thus, instability appears near the center of the bottom surface, as shown in Figure 12(a).

**(a)**

**(b)**

**(c)**

As we reduce the time increment to the maximum stable time increment of the surface boundary subsystem ( = 0.00176 s), the process of wave propagation in the internal region and around the bottom surface can be precisely simulated. However, computational instability occurs when the waves reach the bottom edges, as shown in Figure 12(b), since the stability condition of the edge boundary subsystem is not met ().

Then, the time increment is further reduced to the maximum stable time increment of the edge boundary subsystem ( = 0.00071 s). This time the simulation of the wave propagation can be accurately conducted until the waves arrive at the bottom corners of the model. The instability appears in the corner point, as shown in Figure 12(c), on account of the unsatisfied stability condition of the corner boundary subsystem ().

Finally, the maximum stable time increment of the corner boundary subsystem is adopted for the explicit analysis ( = 0.00032 s). The wave propagation process in the calculation model is shown in Figure 13, and the vertical displacements of observation points A, B, C, and D are shown in Figure 14. It can be seen that the dynamic calculation of the entire model can be successfully executed, and the outgoing waves are effectively absorbed by the artificial boundaries.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The above dynamic analyses verify the precision of the stability conditions proposed in this study. From the comparison between the stability conditions of the internal domain and different substructures, it can be deduced that the numerical stability of the computational system is controlled by the corner region of the model. To validate this deduction, we successively change the important model parameters that would affect the numerical stability, namely, continuously changing the mass density of the internal medium *ρ* from 1000 kg/m^{3} to 2500 kg/m^{3}, the velocity of the shear wave of the internal medium *C*_{S} from 300 m/s to 1000 m/s, the element size *L* from 1 m to 10 m, the model size *R* from 1 m to 1000 m, and the Poisson ratio of the internal medium *μ* from 0.3 to 0.5, and the maximum time increments are investigated under different model parameters. The results are shown in Figure 15.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Under different model parameters, the corner boundary subsystem provides the strictest stability condition, followed by the edge boundary subsystem. In addition, the surface boundary subsystem provides stricter stability conditions than the internal domain under different mass densities, velocities of shear waves, element sizes, and model sizes. However, when the Poisson ratio is relatively large, the stability condition of the internal medium is stricter than that of the surface subsystem. When the shear wave velocity is constant, the compression wave velocity can be expressed as follows:

When the Poisson ratio increases, the compression wave velocity of the internal medium increases rapidly. According to (21), a larger compression wave velocity will lead to stricter stability conditions in the internal medium.

From the above analysis, it can be demonstrated that the numerical stability of the system containing the viscoelastic artificial boundary elements is controlled by the corner region of the model. The stability condition of the corner subsystem given by (40) is a sufficient condition for the numerical stability of the entire model.

#### 6. Conclusion

In this study, a numerical stability analysis method for the explicit algorithm based on the idea of a subsystem is proposed. Through this method, the artificial boundary subsystem that controls the stability of the overall system is determined, and the analytical solution for the stability condition of the explicit integration algorithm considering the 3D viscoelastic artificial boundary elements is derived. Based on the results, the following conclusions can be obtained:(1)For the stability problem of large-scale numerical calculations, the stability analysis of the explicit time-domain stepwise integration algorithm can be performed on the local subsystem. When the internal medium satisfies the stability condition of the algorithm, the numerical stability of the overall system is controlled by the corner region. This stability condition is a sufficient condition for the stability of the overall system, not a necessary and sufficient condition. However, the results of the example verification show that the difference between the sufficient condition and the necessary and sufficient condition is very small. Thus, the sufficient condition can basically meet the needs of the stability analysis in practical applications.(2)When the viscoelastic artificial boundary elements are adopted, the stability of the overall system is different from the internal medium stability due to the influence of the quality, stiffness, and damping of the artificial boundary elements. The former has more stringent numerical stability conditions and requires smaller integration steps to meet overall system stability conditions compared to that of the latter.(3)The stability conditions of the local subsystem are derived based on the format of the central difference, which is commonly used in dynamic numerical analysis. It is important to note that the stability condition is closely related to the integral format. If another time domain step integration format is applied, the method proposed in this paper can also be used by intercepting the representative subsystems and deriving the corresponding transfer matrix and spectral radius. Then, the stability conditions based on the other integral format can be obtained according to a similar process introduced in this paper.(4)The stability conditions given in this paper are derived from regular hexahedron elements and are equally applicable to rectangular hexahedron elements. Since the stability condition is affected by the minimum size of the elements, the minimum side length of the rectangular elements can be used as a parameter to derive the stability condition.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China (nos. 51878384 and U1839201) and the National Key Research and Development Program of China (no. 2018YFC1504305). Financial support from these organizations is gratefully acknowledged.