FE simulation of steady state wave motion in solids combined with a PML approach

(1)

ScienceDirect

Available online at www.sciencedirect.com

Procedia Engineering 199 (2017) 1556–1561

Peer-review under responsibility of the organizing committee of EURODYN 2017.

10.1016/j.proeng.2017.09.054

1877-7058 Available online at www.sciencedirect.com

ScienceDirect

www.elsevier.com/locate/procedia

X International Conference on Structural Dynamics, EURODYN 2017

FE simulation of steady state wave motion in solids combined with a PML approach

Joonsang Park

^a,

* and Amir M. Kaynia

^a

aNorwegian Geotechnical Institute, Sognsveien 72, N-0855, Oslo, Norway

Abstract

This study presents a finite element (FE) framework for simulation of steady state wave motion, combined with a so-called PML approach. The steady state wave motion of current interest are not only time-harmonic load but can also be constant-speed moving load. We adopt an existing PML approach (from low-frequency electromagnetic wave) to apply to the steady state FE simulation of constant-speed moving load vibration as well as temporal harmonic load vibration. In the paper, we will present the PML formulation in both the frequency and steady-state time domains. To validate the approach, then we compare the simulation results with available analytical reference solutions, and with a set of real measured field data (high-speed train data). In addition, we solve more example of steady state vibrations resulting from different train speed (e.g. subcritical, critical and supercritical).

Keywords: finite element; perfectly matched layer; steady-state wave motion; time-harmonic wave motion.

1. Introduction

Since mid 1990's, the so-called perfectly-matched layer method (PML) became popular in its application to numerical simulation of wave propagation by means of discretization methods such as finite difference, finite element, finite volume, etc. In a nutshell, the PML method is to define finite-size physical domains where we transform (or stretch) the real-valued spatial coordinate system into a complex-valued one. The intention is that the incoming wave into the PML domain is absorbed and no significant outgoing wave (reflected back to the main computational domain)

* Corresponding author. Tel.: +47 932 12 453; fax: +47 2223 0448.

E-mail address: [email protected]

Available online at www.sciencedirect.com

ScienceDirect

X International Conference on Structural Dynamics, EURODYN 2017

FE simulation of steady state wave motion in solids combined with a PML approach

Joonsang Park

^a,

* and Amir M. Kaynia

^a

Abstract

1. Introduction

2 Park and Kaynia/ Procedia Engineering 00 (2017) 000–000

occur. Initially, PML was developed for the electromagnetic wave application [2]. Later, it has also been widely applied to acoustic, elastic and seismic waves e.g. [1]. PML application in the frequency domain is rather straightforward and well developed, whereas its application to the time domain is still a research topic.

This study presents a finite element (FE) framework for simulation of steady state wave motion, combined with a PML approach. The steady state wave motion of current interest results from temporal harmonic load and constant- speed moving load. The PML approach in this study was developed originally for the frequency domain application to the ultra-low-frequency (~1Hz) electromagnetic wave motion in the marine environment [5]. Hereby, we adopt the approach in order to apply to the steady state FE simulation of constant-speed traffic load vibration as well as temporal harmonic load vibration. This can be useful to analyze 3D complex and realistic models, including details of e.g.

sleepers, rails, ballasts, any inhomogeneity, etc. The FE framework is implemented into COMSOL MultiPhysics (PDE module). In the paper, we will present the PML formulation in both the frequency and steady-state time domains. To validate the approach, we compare the simulation results not only with available analytical reference solutions [3], but also with a set of real measured field data (high-speed train data) [4]. In addition, we will present more example of verifications from steady-state vibrations (e.g. subcritical, critical and supercritical speed moving loads). At the end of the paper, we will also discuss and evaluate the efficiency of the framework and propose a set of future works in order to improve further the framework.

2. A PML formulation for steady state

In this study, we propose a PML approach that is applicable to elastic wave propagation simulation, particularly for the steady state wave motion resulting from time-harmonic loads and for constant-speed moving loads. The current PML approach was originally developed for so-called ultra-low frequency electromagnetic wave application [5], and hereby is adapted for the current application to elastic waves numerical simulation. The main idea of most PML approaches is to stretch a given real coordinate system into a larger-scale coordinate system, which eventually may represent a virtually-infinite space, resulting in no artificial reflections from computational boundaries. Normally, the stretched coordinate is given in complex value for the time-harmonic loading case (i.e. frequency domain analysis).

However, for the constant-speed loading case, which is simulated in the time domain, the stretched coordinate is in real value. The stretched coordinate of the current study can be expressed by the following simple formula.

xax h

  (1)

where  and x are the stretched and unstretched coordinates, respectively, and h is the (representative) element size in a PML domain. Finally, a is a PML parameter that should be determined according to the wavelength, element size (h), and the PML domain length or size (H) in the unstretched coordinate. Then, the total stretched length (L) of the PML domain can also be given as below.

L  H aH h (2)

Based on numerical investigation, it is found out that L4_P_{, where}_P is the P-wave wavelength. Once knowing L, H, and h, we calculate the parameter a. For the time-harmonic load case, we multiply it by a complex factor of (1- i). For the case of constant-speed moving load, we use the stretching as given above.

2.1. Validation

To demonstrate the performance of the proposed PML approach in the study, we simulate time-harmonic wave propagation in a homogeneous half space (_=1000kg/m³_{, C}_S=1000m/s, CP=2000m/s, 1% damping) that is subjected to a vertical point load, which is converted afterwards into the time-domain seismogram via the fast Fourier transformation. Fig. 1 (a) and (b) show the geometry and mesh of the FE model. Note that the main computational domain is of very small size in the current model. Namely, the size is 1m×1m×60m, which is not of any realistic interest but sufficient enough for the current purpose. The source and receive of interest are placed between x=0m and x=60m, which are also the points at the interfaces between the main domain and PML, along the top surface of the homogeneous half space and on the symmetric plane (y=0 plane). The geometry conditions given in the current model is rather extreme (i.e. very slim main domain), so it can demonstrate the high performance of the proposed PML approach. Plots (a) and (b) in Fig. 2 show the vertical displacement along the x-axis from 0 to 60m for two frequencies

(2)

Joonsang Park et al. / Procedia Engineering 199 (2017) 1556–1561 1557

ScienceDirect

X International Conference on Structural Dynamics, EURODYN 2017

FE simulation of steady state wave motion in solids combined with a PML approach

Joonsang Park

^a,

* and Amir M. Kaynia

^a

Abstract

1. Introduction

ScienceDirect

X International Conference on Structural Dynamics, EURODYN 2017

FE simulation of steady state wave motion in solids combined with a PML approach

Joonsang Park

^a,

* and Amir M. Kaynia

^a

Abstract

1. Introduction

occur. Initially, PML was developed for the electromagnetic wave application [2]. Later, it has also been widely applied to acoustic, elastic and seismic waves e.g. [1]. PML application in the frequency domain is rather straightforward and well developed, whereas its application to the time domain is still a research topic.

This study presents a finite element (FE) framework for simulation of steady state wave motion, combined with a PML approach. The steady state wave motion of current interest results from temporal harmonic load and constant- speed moving load. The PML approach in this study was developed originally for the frequency domain application to the ultra-low-frequency (~1Hz) electromagnetic wave motion in the marine environment [5]. Hereby, we adopt the approach in order to apply to the steady state FE simulation of constant-speed traffic load vibration as well as temporal harmonic load vibration. This can be useful to analyze 3D complex and realistic models, including details of e.g.

sleepers, rails, ballasts, any inhomogeneity, etc. The FE framework is implemented into COMSOL MultiPhysics (PDE module). In the paper, we will present the PML formulation in both the frequency and steady-state time domains. To validate the approach, we compare the simulation results not only with available analytical reference solutions [3], but also with a set of real measured field data (high-speed train data) [4]. In addition, we will present more example of verifications from steady-state vibrations (e.g. subcritical, critical and supercritical speed moving loads). At the end of the paper, we will also discuss and evaluate the efficiency of the framework and propose a set of future works in order to improve further the framework.

2. A PML formulation for steady state

In this study, we propose a PML approach that is applicable to elastic wave propagation simulation, particularly for the steady state wave motion resulting from time-harmonic loads and for constant-speed moving loads. The current PML approach was originally developed for so-called ultra-low frequency electromagnetic wave application [5], and hereby is adapted for the current application to elastic waves numerical simulation. The main idea of most PML approaches is to stretch a given real coordinate system into a larger-scale coordinate system, which eventually may represent a virtually-infinite space, resulting in no artificial reflections from computational boundaries. Normally, the stretched coordinate is given in complex value for the time-harmonic loading case (i.e. frequency domain analysis).

However, for the constant-speed loading case, which is simulated in the time domain, the stretched coordinate is in real value. The stretched coordinate of the current study can be expressed by the following simple formula.

xax h

  (1)

where  and x are the stretched and unstretched coordinates, respectively, and h is the (representative) element size in a PML domain. Finally, a is a PML parameter that should be determined according to the wavelength, element size (h), and the PML domain length or size (H) in the unstretched coordinate. Then, the total stretched length (L) of the PML domain can also be given as below.

L  H aH h (2)

Based on numerical investigation, it is found out that L4_P_{, where}_P is the P-wave wavelength. Once knowing L, H, and h, we calculate the parameter a. For the time-harmonic load case, we multiply it by a complex factor of (1- i). For the case of constant-speed moving load, we use the stretching as given above.

2.1. Validation

To demonstrate the performance of the proposed PML approach in the study, we simulate time-harmonic wave propagation in a homogeneous half space (_=1000kg/m³_{, C}_S=1000m/s, CP=2000m/s, 1% damping) that is subjected to a vertical point load, which is converted afterwards into the time-domain seismogram via the fast Fourier transformation. Fig. 1 (a) and (b) show the geometry and mesh of the FE model. Note that the main computational domain is of very small size in the current model. Namely, the size is 1m×1m×60m, which is not of any realistic interest but sufficient enough for the current purpose. The source and receive of interest are placed between x=0m and x=60m, which are also the points at the interfaces between the main domain and PML, along the top surface of the homogeneous half space and on the symmetric plane (y=0 plane). The geometry conditions given in the current model is rather extreme (i.e. very slim main domain), so it can demonstrate the high performance of the proposed PML approach. Plots (a) and (b) in Fig. 2 show the vertical displacement along the x-axis from 0 to 60m for two frequencies

(3)

1558 Joonsang Park et al. / Procedia Engineering 199 (2017) 1556–1561Author name / Procedia Engineering 00 (2017) 000–000 3

of 50 and 200Hz, respectively. The FE results (blue lines) are compared with a reference solution (red dots) [3]. It is clearly shown that the proposed PML approach performs extremely well.

(a) (b)

Fig. 1. FE model used for the validation: (a) geometry for the whole FE model, and (b) FE mesh zoomed in near the main domain of 1m×1m×60m.

(a) (b)

Fig. 2. Complex-valued vertical displacement along the x-axis from 0 to 60m for two frequencies of (a) 50Hz and (b) 200Hz.

Finally, Fig. 3 compares the results in the time domain for the FE and reference solutions. The time domain solutions are obtained via the fast Fourier transform. The two solutions use the exactly same frequency sampling (i.e.

fmin, fmax, f). Figure 4 shows the comparison between the two solutions, which demonstrate very good agreement.

3. Comparison with measured data

In this section, we attempt to simulate a real case of high-speed train, and compare the FE results with the measured data. The real case measured data was acquired in Ledsgaard (south Sweden) with the X-2000 train model. More details of the real measurement case study can be found in [4]. Fig. 4 shows the load configuration of the X-2000 train, where we apply a point load per bogie. Table 1 shows the soil parameters used in the FE simulation. All this information is taken from [4]. Fig. 4 also shows the FE model applied in the current example, where all the geometrical details of the rail system is implemented. Note that we again apply the symmetric plane approach (on y=0 plane) in

order to save the computational time. The simulation is done directly in the time domain (i.e. no FFT is applied). The train load is applied along the top edges of the rail in the FE model and the response of interest is recorded in the center of one of the ballasts. We consider two different train speeds of V=70km/h (ca. 20m/s) and V=200km/h (ca 55m/s).

Fig. 3. Seismograms of vertical displacement at x=60m, comparing the FE and reference solutions.

Fig. 4. (Left) Train loading configuration for X-2000 run in the Ledsgaard test [4]. Note that the train moves from left to right in the simulation, and each load corresponds each bogie; (Right) FE model applied for the Ledsgaard example.

Table 1. Soil parameters applied to the FE simulation.

Soil layer Thickness [m]

Density [kg/m³]

CS [m/s] CP [m/s]

V= 70 [km/h] V= 200 [km/h] V= 70 [km/h] V= 200 [km/h]

Crust 1.1 1500 72 65 500 500

Organic clay 3.0 1260 41 33 500 500

Clay 4.5 1475 65 60 1500 1500

Clay 6.0 1475 87 85 1500 1500

Half space - 1475 100 100 1500 1500

Fig. 5 (a) and (b) compare the calculated FE results (blue lines) together with the measured data (red lines) for the V=70km/h (~20m/s) and V=200km/h (~55m/s) cases, respectively. The agreement between the two results is not so perfect in magnitude, but the FE result is capturing the general behavior fairly well considering all the uncertainties in the FE model inputs.

(4)

of 50 and 200Hz, respectively. The FE results (blue lines) are compared with a reference solution (red dots) [3]. It is clearly shown that the proposed PML approach performs extremely well.

(a) (b)

Fig. 1. FE model used for the validation: (a) geometry for the whole FE model, and (b) FE mesh zoomed in near the main domain of 1m×1m×60m.

(a) (b)

Fig. 2. Complex-valued vertical displacement along the x-axis from 0 to 60m for two frequencies of (a) 50Hz and (b) 200Hz.

Finally, Fig. 3 compares the results in the time domain for the FE and reference solutions. The time domain solutions are obtained via the fast Fourier transform. The two solutions use the exactly same frequency sampling (i.e.

fmin, fmax, f). Figure 4 shows the comparison between the two solutions, which demonstrate very good agreement.

3. Comparison with measured data

In this section, we attempt to simulate a real case of high-speed train, and compare the FE results with the measured data. The real case measured data was acquired in Ledsgaard (south Sweden) with the X-2000 train model. More details of the real measurement case study can be found in [4]. Fig. 4 shows the load configuration of the X-2000 train, where we apply a point load per bogie. Table 1 shows the soil parameters used in the FE simulation. All this information is taken from [4]. Fig. 4 also shows the FE model applied in the current example, where all the geometrical details of the rail system is implemented. Note that we again apply the symmetric plane approach (on y=0 plane) in

order to save the computational time. The simulation is done directly in the time domain (i.e. no FFT is applied). The train load is applied along the top edges of the rail in the FE model and the response of interest is recorded in the center of one of the ballasts. We consider two different train speeds of V=70km/h (ca. 20m/s) and V=200km/h (ca 55m/s).

Fig. 3. Seismograms of vertical displacement at x=60m, comparing the FE and reference solutions.

Fig. 4. (Left) Train loading configuration for X-2000 run in the Ledsgaard test [4]. Note that the train moves from left to right in the simulation, and each load corresponds each bogie; (Right) FE model applied for the Ledsgaard example.

Table 1. Soil parameters applied to the FE simulation.

Soil layer Thickness [m]

Density [kg/m³]

CS [m/s] CP [m/s]

V= 70 [km/h] V= 200 [km/h] V= 70 [km/h] V= 200 [km/h]

Crust 1.1 1500 72 65 500 500

Organic clay 3.0 1260 41 33 500 500

Clay 4.5 1475 65 60 1500 1500

Clay 6.0 1475 87 85 1500 1500

Half space - 1475 100 100 1500 1500

Fig. 5 (a) and (b) compare the calculated FE results (blue lines) together with the measured data (red lines) for the V=70km/h (~20m/s) and V=200km/h (~55m/s) cases, respectively. The agreement between the two results is not so perfect in magnitude, but the FE result is capturing the general behavior fairly well considering all the uncertainties in the FE model inputs.

(5)

1560 Joonsang Park et al. / Procedia Engineering 199 (2017) 1556–1561

Author name / Procedia Engineering 00 (2017) 000–000 5

(a) (b)

Fig. 5. . Comparison of the calculated FE results (blue lines) together with the measured data (red lines) for V=70km/h (~20m/s) (a) and 200km/h (~55m/s) (b) cases, respectively.

4. Train speed effect

Here we focus our interest on the effect of train speed on the dynamic response of the railway system and in the soil layers (subcritical, critical and supercritical speeds). The physical model of interest is almost the same as one used in Section 3. However, all the soil layers are now assumed to have the same material properties of density=1500kg/m³, S-wave velocity (Cs)=55m/s, P-wave velocity (Cp)=110m/s. This type of soil has a Rayleigh wave velocity of around 51.3m/s. We vary the train speed (V=20, 50, 80m/s) to see its influence on the response distribution as the train moves.

Note that the three train speeds correspond to the subcritical, critical and supercritical cases.

Fig. 6 shows four plots. The first three snapshots (a-c) illustrate the vertical-displacement (color) and 3D deformed- shape at a given particular time for the three different train speeds of V=20, 50 and 80m/s, respectively. The particular time for each snapshot is chosen so that the moving loads are more or less at the same position for easy comparison.

The last plot (d) compares, more directly, the three different train speed cases in terms of the vertical-velocity response snapshot along the railway. Note that the x-axis in plot (d) is the relative distance between the first train load (the right-end load of 320kN in Fig. 4) and the others. It is easily observable that as the train speed increases, the dynamic response in the railway system and the soil increase significantly, particularly for the case of V=80m/s. The increase happens not only near the top surface but also in the soil layers, which can be seen more clearly in plots (a-c). At the same time, we can also notice that the response is significantly vanishing as soon as the steady state wave motion train enters into the PML domains. This confirms that the PML approach proposed in the study works very satisfactorily for the purpose of numerical simulation of the steady state wave motion induced by moving train.

5. Conclusion

We have proposed a PML approach, which can be effectively implemented in any FE framework in order to simulate both the time-harmonic and steady state transient wave propagation. The PML approach consists of a simple stretching function in a form of power law and with only few input parameters. To validate the approach, we have compared several simulation results with 1) an analytical solution based on the continuum wave equation and 2) a set of measured field data. Finally, we have also made a series of numerical simulations with varying the train speed so that we can take into account of the effect of the train speed on the dynamic response of the railway system and the soil layers. It has been clearly observed that as the train speed increases, the dynamic response increases and when the train speed exceeds beyond the Rayleigh critical speed, the whole system vibrates more significantly than the subcritical speed case.

(a) (b)

(c) (d)

Fig. 6. (a) Vertical displacement snapshot for V=20m/s; (b) Vertical displacement snapshot for V=50m/s; (c) Vertical displacement snapshot for V=80m/s; (d) comparison of vertical velocity snapshots for V=20m, 50 and 80/s.

Through various numerical examples and tests, we have confirmed that the PML approach proposed in the study produces excellent results in comparison to reference solutions and measured field data. However, so far we have applied the PML approach to 1) frequency-domain and 2) steady state transient wave motion cases. Additional research is needed for the pure (non-steady state) transient wave motion which requires that the PML approach should be modified and improved. Work is underway by the authors on this development.

Acknowledgements

The authors would like to acknowledge partial support for this study by the research project DESTination Rail - Decision Support Tool for Rail Infrastructure Managers, funded by the European Commission, Grant Agreement 636285 (H2020-MG-2014-2015).

References

[1] J.M. Barbosa, J. Park, and E. Kausel, Perfectly Matched Layers in the Thin Layer Method, Journal of Computer Methods In Applied Mechanics and Engineering, 217–220 (2012) 262–274.

[2] J.-P. Bérenger, (1994) A perfectly matched layer for the absorption of electromagnetic waves, J. Computat. Phys. 114 (2) (1994) 185–200 [3] E. Kausel and J.M. Roësset, Stiffness matrices for layered soils, Bulletin of the Seismological Society of America, 71(6) (1984) 1743-1761.

[4] A.M. Kaynia, C. Madshus, and P. Zackrisson, Ground Vibration from High Speed Trains: Prediction and Countermeasure. J. Geotech. &

Geoenv. Engrg, ASCE, 126(6) (2000) 531-537.

[5] J. Park, T.I. Bjørnarå, and B.A: Farrelly, Absorbing boundary domain for CSEM 3D modelling, Proceedings of COMSOL conference (2010), November 17-19, Paris, France

(6)

(a) (b)

Fig. 5. . Comparison of the calculated FE results (blue lines) together with the measured data (red lines) for V=70km/h (~20m/s) (a) and 200km/h (~55m/s) (b) cases, respectively.

4. Train speed effect

Here we focus our interest on the effect of train speed on the dynamic response of the railway system and in the soil layers (subcritical, critical and supercritical speeds). The physical model of interest is almost the same as one used in Section 3. However, all the soil layers are now assumed to have the same material properties of density=1500kg/m³, S-wave velocity (Cs)=55m/s, P-wave velocity (Cp)=110m/s. This type of soil has a Rayleigh wave velocity of around 51.3m/s. We vary the train speed (V=20, 50, 80m/s) to see its influence on the response distribution as the train moves.

Note that the three train speeds correspond to the subcritical, critical and supercritical cases.

Fig. 6 shows four plots. The first three snapshots (a-c) illustrate the vertical-displacement (color) and 3D deformed- shape at a given particular time for the three different train speeds of V=20, 50 and 80m/s, respectively. The particular time for each snapshot is chosen so that the moving loads are more or less at the same position for easy comparison.

The last plot (d) compares, more directly, the three different train speed cases in terms of the vertical-velocity response snapshot along the railway. Note that the x-axis in plot (d) is the relative distance between the first train load (the right-end load of 320kN in Fig. 4) and the others. It is easily observable that as the train speed increases, the dynamic response in the railway system and the soil increase significantly, particularly for the case of V=80m/s. The increase happens not only near the top surface but also in the soil layers, which can be seen more clearly in plots (a-c). At the same time, we can also notice that the response is significantly vanishing as soon as the steady state wave motion train enters into the PML domains. This confirms that the PML approach proposed in the study works very satisfactorily for the purpose of numerical simulation of the steady state wave motion induced by moving train.

5. Conclusion

We have proposed a PML approach, which can be effectively implemented in any FE framework in order to simulate both the time-harmonic and steady state transient wave propagation. The PML approach consists of a simple stretching function in a form of power law and with only few input parameters. To validate the approach, we have compared several simulation results with 1) an analytical solution based on the continuum wave equation and 2) a set of measured field data. Finally, we have also made a series of numerical simulations with varying the train speed so that we can take into account of the effect of the train speed on the dynamic response of the railway system and the soil layers. It has been clearly observed that as the train speed increases, the dynamic response increases and when the train speed exceeds beyond the Rayleigh critical speed, the whole system vibrates more significantly than the subcritical speed case.

(a) (b)

(c) (d)

Fig. 6. (a) Vertical displacement snapshot for V=20m/s; (b) Vertical displacement snapshot for V=50m/s; (c) Vertical displacement snapshot for V=80m/s; (d) comparison of vertical velocity snapshots for V=20m, 50 and 80/s.

Through various numerical examples and tests, we have confirmed that the PML approach proposed in the study produces excellent results in comparison to reference solutions and measured field data. However, so far we have applied the PML approach to 1) frequency-domain and 2) steady state transient wave motion cases. Additional research is needed for the pure (non-steady state) transient wave motion which requires that the PML approach should be modified and improved. Work is underway by the authors on this development.

Acknowledgements

The authors would like to acknowledge partial support for this study by the research project DESTination Rail - Decision Support Tool for Rail Infrastructure Managers, funded by the European Commission, Grant Agreement 636285 (H2020-MG-2014-2015).

References

[1] J.M. Barbosa, J. Park, and E. Kausel, Perfectly Matched Layers in the Thin Layer Method, Journal of Computer Methods In Applied Mechanics and Engineering, 217–220 (2012) 262–274.

[2] J.-P. Bérenger, (1994) A perfectly matched layer for the absorption of electromagnetic waves, J. Computat. Phys. 114 (2) (1994) 185–200 [3] E. Kausel and J.M. Roësset, Stiffness matrices for layered soils, Bulletin of the Seismological Society of America, 71(6) (1984) 1743-1761.

[4] A.M. Kaynia, C. Madshus, and P. Zackrisson, Ground Vibration from High Speed Trains: Prediction and Countermeasure. J. Geotech. &

Geoenv. Engrg, ASCE, 126(6) (2000) 531-537.

[5] J. Park, T.I. Bjørnarå, and B.A: Farrelly, Absorbing boundary domain for CSEM 3D modelling, Proceedings of COMSOL conference (2010), November 17-19, Paris, France