Fast 3D FEM-BEM coupling for dynamic soil-structure interaction

Available online at www.sciencedirect.com Available online at www.sciencedirect.com

ScienceDirect ScienceDirect Procedia Engineering 00 (2017) 000–000 Available online at www.sciencedirect.com Procedia Engineering 00 (2017) 000–000

ScienceDirect

www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia

Procedia Engineering 199 (2017) 391–396

X International Conference on Structural Dynamics, EURODYN 2017 X International Conference on Structural Dynamics, EURODYN 2017

Fast 3D FEM-BEM coupling for dynamic soil-structure interaction Fast 3D FEM-BEM coupling for dynamic soil-structure interaction Winfried Schepers Winfried Schepers

GuD Geotechnik und Dynamik Consult GmbH, Darwinstr. 13, 10589 Berlin, Germany GuD Geotechnik und Dynamik Consult GmbH, Darwinstr. 13, 10589 Berlin, Germany

Abstract Abstract We propose a fast method for solving soil-structure interaction problems in frequency domain. Finite Elements are applied to We propose fast method forboundary solving soil-structure interaction problemsthe in interface frequencybetween domain.the Finite Elements areunbounded applied to discretize theastructure, while elements are applied to discretize structure and the discretize thelayered structure, while boundary elements are applied discretize theelement interface between and fully the unbounded horizontally soil. Some mild restrictions imposed ontothe boundary mesh allowthe forstructure storing the populated horizontally layered soil. Some mild O(N) restrictions the boundary element mesh for storing the fully problem populateda N × N soil flexibility matrix in only storageimposed and in on negligible time. For solving theallow soil-structure interaction N × N soil flexibility O(N) in negligible time. For solving the soil-structure interaction problem hybrid linear equationmatrix systeminisonly solved withstorage direct and and indirect solvers. Though our current implementation of the method in aa hybrid linearFE equation system is solved with direct and the indirect solvers. Though implementation of the that method in a commercial code does not allow for fully exploiting superior properties of our the current proposed method, it is shown iterative commercial FE codetodoes notsolvers allow for for the fully exploiting thesoil-structure superior properties of the proposed method, ithere. is shown that iterative solvers are superior direct solution of the interaction problems investigated solvers are superior to direct solvers for the solution of the soil-structure interaction problems investigated here. © 2017 The Authors. Published by Elsevier Ltd. © 2017 Published Elsevier Ltd. © 2017The TheAuthors. Authors. Publishedby by Elsevier Ltd. committee of EURODYN 2017. Peer-review underresponsibility responsibility ofthe theorganizing organizing Peer-review under of committee of EURODYN 2017. Peer-review under responsibility of the organizing committee of EURODYN 2017. Keywords: Soil-structure interaction; direct and iterative solvers; substructure; frequency domain; sparse and dense matrices Keywords: Soil-structure interaction; direct and iterative solvers; substructure; frequency domain; sparse and dense matrices

1. Introduction 1. Introduction In dynamic soil-structure interaction problems two intrinsically different types of domains are involved. On the dynamic interactione.problems two intrinsically types of domains are involved. the oneInhand theresoil-structure is the finite structure, g. a building subjected to different impinging vibrations emanating from a On nearby one hand there is the finite structure, e. g. a building subjected to impinging vibrations emanating from a nearby train track on the soil surface. The Finite Element Method (FEM) is perfectly suited to solve for the displacement train track on the structure soil surface. The Finite Element Method (FEM)loads. is perfectly solve for isthethedisplacement and forces in the induced by external static or dynamic On thesuited other to hand there unbounded and forces in the structure induced by external static or dynamic loads. On the other hand there is or thea unbounded soil, having a stress free surface, being either homogeneous or layered, either an infinite half-space stratum on soil, havingrock. a stress free surface, homogeneous layered, either for an infinite half-space or a stratum on competent Several methodsbeing have either been proposed in theorpast to account the radiation of waves to infinity, competent rock. Several methods have been proposed in the past to account for the radiation of waves to infinity, e. g. Local Transmitting Boundaries [1], Perfectly Matched Layers (PML) [2], and the Boundary Element Method e. g. Local Boundaries (PML) and the of Boundary Element Method (BEM) [3], Transmitting amongst others. The order[1], in Perfectly which theMatched methods Layers are listed here [2], is the order increasing accuracy, and (BEM) [3], amongst others. The order in which the methods are listed here is the order of increasing accuracy, concurrently the order of decreasing compatibility with the fundamental concepts of the Finite Element Method. and concurrently the order of decreasing compatibility with the fundamental concepts of the Finite Element Method.

1877-7058 © 2017 The Authors. Published by Elsevier Ltd. 1877-7058 2017responsibility The Authors. of Published by Elsevier Ltd. of EURODYN 2017. Peer-review©under the organizing committee Peer-review under responsibility of the organizing committee of EURODYN 2017.

1877-7058 © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the organizing committee of EURODYN 2017. 10.1016/j.proeng.2017.09.129

Winfried Schepers / Procedia Engineering 199 (2017) 391–396 W. Schepers/ Procedia Engineering 00 (2017) 000–000

392 2

Obviously, some sacrifices are unavoidable if wave radiation to infinity is to be cast in the FEM concept. If accuracy is the priority, some remedy has to be found for the significant computational cost arising from the application of BEM for solving dynamic soil-structure interaction problems. In the following sections a method is presented which allows FEM-BEM coupling of structures at the surface of an arbitrarily layered half-space. at very low computational cost, with only mild restrictions imposed to the BEM mesh. In section 2 the equation of motion of the structure, the soil, and the coupled soil-structure system in frequency domain is briefly derived. Section 3 presents the strategy to accelerate the solution of the coupled system, and section 4 provides results from an application, followed by an outlook in section 5. 2. FEM-BEM coupling 2.1. Generic coupled equation of motion The generic linear elastic equation of motion of the structure in frequency domain is

(− Ω 2M + iΩ C + K ) ⋅ u = P − Q   

(1)

S

in which u is the vector of displacements, P the vector of external forces, Q the vector of soil-structure interaction forces, K the elastic stiffness matrix, C the damping matrix, M the mass matrix, Ω the circular frequency, and i = − 1 . Introducing the concept of a “complex valued dynamic stiffness matrix” S, partitioning by degree of freedom (DOF) inside the structure (“R”) and on the soil-structure interface (“I”) we can write

 S RR   S IR

S RI   u R   PR   0  ⋅  =   −   S II   u I   PI   QI 

(2)

From evaluation of the discretized Boundary Element integral equation of the soil domain the generic stressdisplacement relationship

v −s = G ⋅q

(3)

is obtained, with v the soil displacements, s the vector of free-field motions, q the vector of soil stresses, and G the soil flexibility matrix. Equilibrium and compatibility at the soil-structure interface can be enforced in a weak sense only due to non-matching meshes. The displacement compatibility can conveniently be written as v = Tvu ⋅ u I ,

(4)

and the equilibrium condition as (5)

Q I = TQq ⋅ q

with transformation matrices Tvu and TQq. Equations (2)-(5) can be gathered in a block matrix, and after elimination of v and QI we obtain  S RR   S IR   0

S RI S II Tvu

0   u R   PR       TQq  ⋅  u I  =  PI  ,  − G   q   s 

(6)



Winfried Schepers / Procedia Engineering 199 (2017) 391–396 W. Schepers/ Procedia Engineering 00 (2017) 000–000

393 3

which is the starting point for any further runtime analyses and optimizations of the solution of the SSI problem. Note that the hybrid matrix equation (6) contains FEM displacements as well as BEM stresses as unknowns. 2.2. FEM Implementation, BEM implementation, and coupling conditions For our implementation we used ANSYS v14 [4] as the base solver engine. The building was discretized using the ANSYS element library. All other parts of the matrix and load vector in eq. (6) are imported into ANSYS as a substructure matrix. For evaluation of the BEM integral equations we use as fundamental solutions half-space Green’s functions obtained by means of the Thin Layer Method [5], computed once in advance, and stored in a lookup table while assembling the soil flexibility matrix G. Compatibility at the soil-structure interface is enforced by mapping the FEM node displacements on the BEM node displacements. FEM element shape functions are evaluated to provide the coefficients of the mapping matrix Tvu in eq. (4). For enforcing equilibrium we obtain from work-equivalence considerations TQq = (Tvu)T · A

(7)

where A is a narrow banded real valued matrix which maps the soil stresses to their equivalent nodal forces at the BEM nodes. For more details on the implementation see [6, 7, 8]. 2.3. Strategies for soil stiffness matrix assembly and for solving the matrix equation of motion The distinctive detrimental property of eq. (6) is the fact that the matrix on the left hand side is densely populated by virtue of the soil flexibility matrix G. For example, the shell mesh of the building in fig. 1a comprises approximately 220000 DOF including the slab footing, and the corresponding dynamic stiffness matrix contains 2.3 × 106 non-zero elements (fig. 1b). The soil flexibility matrix increases the number of non-zero elements to approximately 20 × 106 (fig. 1c). a)

b)

c)

Fig. 1. Vanishing sparsity of dynamic stiffness by soil flexibility. (a) finite element mesh. (b) non-zero elements of dynamic stiffness matrix of structural elements. (c) non-zero elements of dynamic soil-structure stiffness matrix

Common Gaussian elimination type FE solvers feature sophisticated algorithms to preserve the sparsity pattern during matrix decomposition, such that problem sizes remain solvable though storing the decomposed matrix would by far exceed the memory available. If, however, these algorithms are applied to densely populated matrices, they might break down. At this point the assembly of the soil flexibility matrix deserves some attention. Because numerical integration of complex valued transcendental functions has to be performed, the required computational time must not be neglected. From the literature, the Fast Multipole Method [9] and Hierarchical Matrices [10] are well known procedures which aim a reducing the number of function evaluations. Strategies to solve eq. (6) seek either to reduce the storage requirements, or the reduce the required number of floating point operation, or a combination of both. The strategy applied in [7] was to straightforwardly compute the soil-stiffness matrix from inverting G, and to arrive at an foundation stiffness matrix comprising the sparsely

Winfried Schepers / Procedia Engineering 199 (2017) 391–396 W. Schepers/ Procedia Engineering 00 (2017) 000–000

394 4

S populated footing stiffness matrix S II ≡ S R II and the fully populated soil stiffness matrix S II . Standard FE solvers can be used to solve the corresponding matrix equation. Iterative coupling detaches the soil from the structure. After applying an initial guess for either displacements or stresses at one subdomain, the solution of the subdomain matrix equation provides an estimate for the corresponding quantity at the interface of the other subdomain, and so forth. For time harmonic problems the convergence of iterative coupling methods is hardly guaranteed [11], in particular if the stiffness contrast between substructures is very small and floating substructures are present [8, 10].

3. Optimization

We propose a method to assemble and solve a soil-structure interaction problem in frequency domain as given in eq. (6). The particular features are: • Consumes only O(M) storage for the soil flexibility matrix, with M the number of elements, and O(M) additional storage for an index table • CPU time consumption for flexibility matrix assembly is negligible • Allows – and strongly calls for – using iterative solvers for the solution of a soil-structure interaction problems • Only mild restrictions are imposed on the BEM mesh geometry

The following requirements must be fulfilled to use the method: • BEM mesh comprises rectangular elements with constant shape functions, i. e. one node per BEM element, and BEM nodes are at element centers • BEM elements must be of equal size and orientation throughout • BEM elements are arranged in a rectangular grid, though not each grid position must be occupied • Distance between two BEM nodes in each coordinate direction is an integer multiple of the element width in that direction.

For the sake of brevity we derive the algorithm for a single pair of source DOF and receiver DOF only, instead of all 3 × 3 = 9 pairs. Furthermore, we assume that each position of the element grid is occupied, and we neglect for now the fact that the sign of a flexibility matrix element for some DOF pairs depends on the sign of the signed horizontal distances between source and receiver point. The method exploits the fact that for a regular BEM mesh fulfilling the conditions listed above the number of unique matrix elements is

N = na ⋅ n b with n a and n b being the number of nodes in X- and Y-direction, respectively, while the matrix comprises a total of N² matrix elements, see fig. 2. Hence, the storage requirement for the N × N flexibility matrix is O(N) only. For the assembly of the full matrix an index matrix is required. Having computed the horizontal distances between a source point and a receiver point in terms of integer multiples of the element widths, these integer number can be used as the array index to obtain the storage location of the corresponding unique flexibility matrix element. The index matrix is also used to keep track of the sign of the flexibility matrix elements. The pseudo code for assembly of the full flexibility matrix can be found in Appendix A. Solving eq. (6) with a direct solver would nullify the advantages of the low storage requirements, because the full flexibility matrix must be assembled in order to use generic solver routines. To exploit the very low storage requirements for the flexibility matrix, iterative solvers are the best choice. Because iterative solvers iterate over the results of matrix-vector products, for optimal performance a subroutine must be provided which exploits the sparsity of the FEM matrix as well as the low storage requirement of the flexibility matrix to provide iterates for the solver.



W. Schepers/ Procedia Engineering 00 (2017) 000–000 Winfried Schepers / Procedia Engineering 199 (2017) 391–396

5 395

Fig. 2. BEM mesh with na = 2 and nb = 4, unique distances between nodes, unique flexibility matrix elements, and index matrix for fast and data sparse assembly of the flexibility matrix

At this point we unfortunately reached the limits of the user programmable features offered by commercial FE software to manipulate the solution process and the FEM matrices. First, for solving eq. (6) with a direct solver, the flexibility matrix G must be passed to ANSYS as a sub-block of a substructure matrix alongside TQq and Tvu and a block filled with zeroes. Yet, ANSYS expects any substructure matrix to be a densely populated, and does not exploit the strong sparsity of the transformation matrices. Furthermore, ANSYS does not provide an interface to implement a user defined matrix-vector product subroutine for its iterative solvers, nor an interface to implement more equation solvers. The implemented iterative solvers can only be applied to eq. (6). 4. Application

Due to the limitations of our current implementation environment with ANSYS as the base solver engine, we are not yet able to fully exploit the advantages of the proposed method. For a preliminary check of the convergence properties when using an iterative solver, we applied ANSYS’ ICCG solver (Incomplete Cholesky Conjugate Gradient) to solve eq. (6). We compare the performance with the solution of the same soil-structure problem, but solved by first inverting the flexibility matrix, following the procedure in [7]. The problem was an embankment on an homogeneous half-space excited by some forces applied to the half-space surface. The results are summarized in Appendix B. The time required to assemble the soil flexibility matrix is negligible. It can be observed that the iterative solver performed well despite the large number of non-zero elements. 5. Summary and outlook

We proposed a fast method for solving a soil-structure interaction method in frequency domain. Finite Elements are applied for discretizing the structure, while boundary elements are applied to discretize the interface between the structure and the unbounded horizontally layered soil. Some mild restrictions imposed on the boundary element mesh allows for storing the fully populated N × N soil flexibility matrix in only O(N) storage in negligible time. Though the implementation in a commercial FE code did not allow to fully exploit the superior properties of the proposed method, it was shown that iterative solvers are superior to direct solvers for the solution of the soilstructure interaction problem. For fully taking advantage of the low storage requirements the method will be cast in a framework for solving complex valued linear equation systems, which allows for the implementation of sparse substructure matrices to be used with direct solvers, and which also allows to implement dedicated subroutines for matrix-vector products as well as dedicated preconditioners to be used with iterative solvers.

Winfried Schepers / Procedia Engineering 199 (2017) 391–396 W. Schepers/ Procedia Engineering 00 (2017) 000–000

396 6

Appendix A. Flexibility matrix assembly pseudo code 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020

Let da = Element width in X Let db = Element width in Y Let na = Element count in X Let nb = Element count in Y Let N = na × nb

021 022 023 024 025 026 027 028 029 030 031 032 033 034 035 036 037 038 039 040 041

Allocate index matrix J with dimensions (0 ... na – 1)×(0 ... nb – 1) × 2 ! Initialize index matrix Let Jn,m = {0,0} Allocate flexibility matrix G with dimension N×N ! Loop over rows of flexibility matrix DO i = 1 TO N ! Loop over all columns of lower triangular part of flexibility matrix DO j = 1 TO i ! Compute indices of distance between source node i and receiver node j

Let n = |xi – xj|/da Let m = |yi – yj|/db ! Compute flexibility matrix IF Jn,m = {0,0} Compute Gi,j by integration of Green’s functions Let Jn,m = {i · SIGN(xi – xj), j · SIGN(yi–yj)} ELSE Let k = Jn,m,1 Let l = Jn,m,2 Let Gi,j = SIGN(k) × SIGN(l) × Gk,l ENDIF ENDDO (j) ENDDO (i) DO i = 1 TO N DO j = i + 1 TO N Let Gi,j = Gj,i ENDDO ENDDO

Appendix B. Model size and runtime performance of example problem Size indicator

Solution strategy Solution of after [6] eq. (6)

Runtime performance

Solution strategy after [6]

Solution of eq. (6)

FEM DOFs

25893

25893

Assembly of flexibility matrix

1s

1s

FEM DOFs at soil-structure interface

3024

3024

Inversion of flexibility matrix

2406 s

n/a

BEM DOFs

11160

11160

Number of iterations of ICCG

440

1524

Total number of DOFs

25893

37053

Wall clock time

3148 s

2862 s

1945 MB

12606 MB

RAM used

References [1] Kausel E. (1988), Local transmitting boundaries, J. Eng. Mech. 114(6) 1011–1027 [2] Kausel E., de Oliveira Barbosa J. M. (2012) PMLs: A direct approach. Int J Numer Meth Eng 90(3):343–352. doi: 10.1002/nme.3322 [3] Bode C., Hirschauer R., Savidis S. A. (2002) Soil–structure interaction in the time domain using halfspace Green's functions. Soil Dyn Earthq Eng 22(4):283–295. doi: 10.1016/S0267-7261(02)00020-9 [4] ANSYS. www.ansys.com [5] Kausel E. (1986) Wave propagation in anisotropic layered media. Int J Numer Meth Eng 23(8):1567–1578. doi: 10.1002/nme.1620230811 [6] Savidis S. A., Hirschauer R., Bode C., Schepers W. (2002) 3D-Simulation of Dynamic Interaction Between Track and Layered Subground. In: Popp Karl (Eds.) System dynamics and long term behaviour of railway vehicles, track and subgrade. Springer, pp 431–450 [7] Hirschauer, R. (2001) Kopplung von Finiten Elementen mit Rand-Elementen zur Berechnung der dynamischen Baugrund-BauwerkInteraktion. Dissertation, Technische Universität Berlin. ISBN 3-7983-1883-2 (in German) [8] Schepers, W. (2014) Fast solution methods for large-scale soil-structure interaction problems in frequency domain. Dissertation, Technische Universität Berlin. DOI 10.14279/depositonce-4240 (in German) [9] Grasso E., Chaillat S., Bonnet M., Semblat J.-F. (2012) Application of the multi-level time-harmonic fast multipole BEM to 3-D viscoelastodynamics. Eng Anal Bound Elem 36(5):744–758. doi: 10.1016/j.enganabound.2011.11.015 [10] Coulier P., François S., Lombaert G., Degrande G. (2014) Coupled finite element – hierarchical boundary element methods for dynamic soil–structure interaction in the frequency domain. Int J Numer Meth Eng 97(7):505–530. doi: 10.1002/nme.4597 [11] Collino F., Ghanemi S., Joly P. (2000) Domain decomposition method for harmonic wave propagation: a general presentation. Comput Method Appl M 184(2-4):171–211. doi: 10.1016/S0045-7825(99)00228-5