© S. Li et al., Published by EDP Sciences 2025

1 Introduction

Bridge cranes are handling devices located usually overhead in the buildings, which are employed to lift and transfer heavy and critical loads. In nuclear power plants (NPPs), main cranes are in the reactor building and the turbine building. A bridge crane is a permanent part of pressurized water reactors: it is not only used during the operational phase of the power station for refuelling and maintenance, but also to install heavy equipment during construction. For example, the giant polar crane in the first reactor building of Hinkley Point C, a new NPP built in UK, will be used to lift the reactor vessel and steam generators during construction and after for refuelling and maintenance when the power station is operational. The polar crane is almost 45 m long and weighs 750 tons. It can rotate 360° on circular rails inside the dome and will remain part of the nuclear plant for the duration of its operating life.

The seismic behavior of bridge cranes is an important safety issue for industrial facilities, especially nuclear plants. Within the context of seismic Level 1 Probabilistic Safety Assessment (PSA) studies performed on nuclear power plants, bridge crane's failure has been identified as a significant contributor in the probability of core meltdown [1]. Therefore, the seismic response of bridge cranes is an essential issue for nuclear safety, requiring experimental data and validated numerical simulations. However, modelling the non-linear dynamic behaviour of bridge cranes under earthquake is a challenging exercise in engineering, due to non-linear effects, such as frictional contacts and impacts between the trolley wheels and the girder beams during an earthquake, as well as the wheel-rail interaction between the end truck beams and the runway beams. This paper focuses on the non-linear dynamic behavior of a 1/5 scale overhead crane under seismic loading. Three-dimensional numerical simulation will be set up in order to predict the seismic behavior of the bridge crane, in comparison to the experimental results. The mock-up of the crane, including runway beams, end truck beams, girder beams and a trolley, was tested in 2015 on the shake table AZALEE of French Sustainable Energies and Atomic Energy Commission (CEA). The dominant failure modes of a crane concern the uplift and the derailment, requiring an accurate modelling of impacts and frictional contacts experienced by the wheel-rail interaction during the seismic excitation. Thus, in order to well reproduce the seismic response of the mock-up, an efficient non-linear transient dynamic analysis is necessary with an adequate modelling of the impacts and frictional contacts in the three-dimensional detailed mesh, because the classical linear analysis always considers that the components of the bridge crane are rigidly attached to each other [2]. The absence of the frictional contacts in linear analysis also leads to the overestimation of the transmitted efforts to the supporting structure. With respect to more and more rigorous seismic design requirements, advanced non-linear dynamic numerical techniques are required to take into account such non-linear effects. Moreover, a full 3D approach is more suitable to describe the contacting body surfaces for predicting impact and sliding contacts as well as stresses in the bodies. It allows also to avoid rotational degrees of freedom as present in beam finite elements. Nevertheless, only some simplified numerical models in this field have been employed in the literature due to the important computation time required by refined mesh with 3D finite elements. Sarh et al. [2] used a simplified numerical model composed of Timoshenko beam elements to analyze the earthquake behavior of a bridge crane, with the possibility of sliding. Otani et al. [3] studied the vertical dynamic response of an overhead crane under earthquake by 1/8 scale model excitation test and time history analysis. Schukin and Vayandrakh [4] developed a coupled FE Model between polar crane and supporting ring wall to simulate the seismic behavior of a polar crane, including shell elements, beam elements and linear springs. Suzuki et al. [5] focused on the bi-dimensional dynamic response of an overhead crane under earthquake, using beam elements and gap elements to consider slip and leap in ABAQUS. Feau et al. [6] established a numerical model with beam elements to investigate the seismic response of a 1/5 scale bridge crane and the sliding contact is taken into account via a penalty approach. Nguyen et al. [7,8] investigated experimentally and numerically the seismic response of a 1/20 scale container crane through shake table test and carried out finite element analyses using beam elements. The derailment and the uplift of the crane, corresponding to the dominant failure modes during the shake table test, are reproduced by introducing a support element at the bottom leg of the container crane, characterized by a friction coefficient combined with springs. Kettler et al. [9,10] carried out static computations in order to explore wheel eccentricity on the runway beam, by employing 3D elements and hard-contact option with a friction coefficient between the wheels and rails, using the software Abaqus. Zhang et al. [11] considered the contact-impact problem as a linear complementarity problem (LCP) in order to treat non-smooth dynamics of a polar crane using an implicit time integration in MATLAB environment. The global behaviour of the containment vessel and its overhead polar crane was obtained by following a flexible multibody approach with modal synthesis for the different components of the problem. In the literature, though some comparisons were carried out between experimental and numerical results, only approximate agreement can be achieved, due to the well-known sensitivity of the non-smooth dynamic response to input uncertainties. In addition, some key features of the dynamic response of bridge crane were not analyzed in the literature, such as lateral impacts/contacts, sliding behavior and efforts transmitted to the anchorages.

In this paper, with the aim of predicting numerically the seismic non-linear behavior of an experimentally tested bridge crane in comparison with the experimental results provided by the SOCRAT Benchmark [12], an efficient explicit time integrator, called CD-Lagrange scheme [13–16], based on the Central Difference (CD) scheme and a velocity-impulse formulation [17–20], is adopted to reproduce the impacts and frictional contacts experienced by the wheel-rail interaction during earthquake. Concerning the resolution of the contact constraints, the advantage of CD-Lagrange algorithm is that there is no iteration in the time step and no need to define additional numerical parameters, on the contrary of penalty methods and implicit time integrations which are widely used in finite element analysis to solve contact problems [21–23]. As a detailed 3D mesh of the mock-up is considered, a subdomain strategy is set up to reduce the computation time by coupling CD-Lagrange scheme with Newmark implicit time integration scheme with a large time scale. The coupled non-smooth transient explicit/implicit problem will be dealt in MATLAB environment with a dual approach using non-linear substructuring in the framework of Heterogeneous Asynchronous Time Integrator (HATI), which allows us to choose appropriate time integrators and time step sizes according to the subdomains [24–27]. In the case of the crane under seismic excitation, only small impact/contact patches will be dealt with the CD-Lagrange scheme, associated with a fine time step satisfying the CFL condition, while the unconditionally stable Newmark implicit scheme can be adopted with a large time step in the other parts to reduce the computation time [28]. In addition, thanks to the hybrid approach, the introduction of the Rayleigh viscous matrix will have no negative effect on the critical time step of Central Difference scheme in the explicit subdomain [29], because Rayleigh damping will be adopted only in the implicit subdomain, without reducing the critical time step of the explicit subdomain. Here, the HATI strategy is applied to the mock-up of the bridge crane for reproducing the shake table tests. Numerical results are compared with experimental results provided by the SOCRAT benchmark in terms of vertical displacement of beams, pseudo-acceleration spectra, relative sliding displacements and forces transmitted to the supports. The different comparisons between numerical and experimental results give an overall satisfactory agreement, demonstrating the ability of the proposed numerical model to reproduce the non-linear seismic behavior of bridge cranes.

The paper is organised as follows. First, the implementation of the subdomain approach in MATLAB environment is briefly reminded on the basis of the 3D mesh of the crane mock-up. It concerns the CD-Lagrange scheme to model for the wheel-rail interaction in explicit patches and the HATI framework to couple explicit and implicit time integrations schemes with their own time step. Before the comparison between numerical and experimental results from the SOCRAT Benchmark, a numerical validation under gravity loading is given by comparing HATI results with reference results obtained by full explicit computation using only the CD-Lagrange scheme. Then, the numerical/experimental comparison is carried for different dynamic tests: an impulsive test according to the transverse direction of the crane for calibrating the friction coefficient of the CD-Lagrange algorithm, a three-dimensional seismic test with a moderate acceleration level and finally a three-dimensional seismic test with a high acceleration level.

2 Multi-time step explicit/implicit subdomain coupling with frictional contacts and high-frequency Rayleigh damping

In this section, the implementation of the hybrid approach including frictional contacts and Rayleigh damping to study the non-linear seismic behavior of bridge cranes is presented. The weak form of the coupled hybrid problem after time and space discretization is written, by focusing in particular on the CD-Lagrange scheme for frictional contact/impact problems. Small explicit partitions are considered, including the wheels and small parts of the rails, dealt with the CD-Lagrange scheme, whereas the main partition of the problem is dealt with an implicit time integration with a large time step and a viscous Rayleigh damping with stiffness matrix contribution. The GC method [24,25] is employed in order to couple the CD-Lagrange scheme with the classical Newmark implicit time integration scheme, the Constant Average Acceleration (CAA) scheme. The GC method and its extended versions have been successfully applied in non-linear dynamics and hybrid testing [30–33]. Indeed, prescribing velocity continuity at the interface, it has been proved that the GC coupling method is stable for all time integrators belonging to Newmark family (implicit and explicit), with the time steps depending on the partitions.

2.1 The coupled hybrid asynchronous problem

The coupled hybrid asynchronous problem can be written in the following discrete form in space and time with the imposed velocity continuity at the interface between the explicit and implicit subdomains: $$\{\begin{array}{l}\begin{array}{c}\hfill M_1{{\displaystyle \stackrel{\mathrm{..}}{U}}}_1^j+C_1{{\displaystyle \dot{U}}}_1^j+K_1{U}_1^j=F_1^{ext,j}+F_1^{contact,j}+L_1^T{\lambda }^j\text{for}j=1,\dots ,m\hfill \end{array}\\ M_2{{\displaystyle \stackrel{\mathrm{..}}{U}}}_2^m+C_2{{\displaystyle \dot{U}}}_2^m+K_2{U}_2^m=F_2^{ext,m}+L_2^T{\lambda }^m\\ L_1{{\displaystyle \dot{U}}}_1^j+L_2{{\displaystyle \dot{U}}}_2^j=0\text{for}j=1,\dots ,m\\ \begin{array}{l}\begin{array}{l}\text{Normal contact conditions}\hfill \end{array}\hfill \\ \text{Tangential frictional conditions}\hfill \end{array}\end{array}$$(1)

where M₁ and M₂ are the mass matrices, K₁ and K₂ are the stiffness matrices, C₁ and C₂ are the Rayleigh damping matrices. The mass matrix M₁ in the explicit subdomain is lumped to make the matrix diagonal. Thanks to the hybrid approach, Rayleigh viscous matrix is only applied in the implicit subdomain (C₁=0). The viscous Rayleigh damping can be defined as C₂ = α_MM₂ + α_KK₂, with both mass matrix and stiffness matrix. The relationship between the damping ratio ξ and two Rayleigh damping coefficients is given by $\xi =\frac{{\alpha }_M}{2\omega }+\frac{{\alpha }_K\omega }{2}$, where ω is the angular frequency. $F_1^{ext}$ and $F_2^{ext}$ are the external forces and $F_1^{contact}$ is the contact force deduced from the normal-contact and tangential-frictional conditions expressed as Lagrange multipliers in CD-Lagrange scheme reminded in the following subsection.

As shown in Figure 1, the implicit subdomain equilibrium is prescribed at the time t_m at the end of the large time Δt₂, whereas the explicit subdomain equilibrium is prescribed at every fine time scale t_j = jΔt₁ (j = 1, … , m), m being the time step ratio between two partitions. The velocity continuity at the interface between explicit and implicit subdomains is ensured at the fine time scale by introducing Lagrange multiplier λ. L₁ and L₂ are the Boolean constraint matrices at the interface with matching meshes as considered in this paper, which pick out the degrees of freedom at the interface between explicit and implicit subdomains to impose the velocity continuity. Concerning non-matching meshes, the mortar approach needs to be employed; L₁ and L₂ will be non-Boolean constraint matrices [34,35].

Fig. 1 Relationship between Δt1 and Δt2, m being the time step ratio between the two partitions.

2.2 Time integration in the explicit subdomain with frictional contacts and impacts

In the following, the CD-Lagrange scheme is applied in the explicit subdomain with impacts and frictional contacts. In fact, the Central Difference scheme will be combined with the non-smooth dynamic approach in which Lagrange multipliers are introduced to impose the kinematic constraints. The normal contact conditions and the tangential frictional conditions expressed in velocity are imposed in a weak way, by introducing Lagrange multipliers. The approach is classified as an event-capturing method which allows to take into account an accumulation of impacts in a given time step, contrarily to the family of event-tracking time integration methods which require the accurate detection of the impact times. The advantage of CD-Lagrange time integrator is that there is no iteration in the time step and no need to define additional numerical parameter. The velocity-impulse formulation [16–19] is employed to treat the frictional contact. Multiplying explicit subdomain equilibrium by its time step Δt = t_j − t_j−1 and using the Central Difference approximation for the acceleration: ${{\displaystyle \ddot{U}}}^j=\left({{\displaystyle \dot{U}}}^{j+\frac{1}{2}}-{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}\right)/\text{Δ}t$, the following dis-crete non-smooth dynamic equation can be obtained: $$M{{\displaystyle \dot{U}}}^{j+\frac{1}{2}}=M{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}+\text{Δ}t\left(F^{ext,j}-C{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}-K{U}^j\right)+I^j$$(2)

where M is a lumped mass matrix obtained by the classical row sum technique [36], ${{\displaystyle \dot{U}}}^{j-\frac{1}{2}}$ is the known mid-step velocity at the current time step, ${{\displaystyle \dot{U}}}^{j+\frac{1}{2}}$ is the mid-step velocity at the next time step, F^ext,j is the external force at the end of the time step and I^j is an impulse term related to the impact/contact. Here, to reduce the amount of notation, we have omitted the subscript 1, which identifies the explicit subdomain in the global discrete coupling problem given in equation (1). In equation (2), we have also omitted the interface force $L_1^T{\lambda }^j$, corresponding to the Lagrange multiplier introduced at the interface between the explicit and implicit subdomains, which will be dealt in the following subsection according to the GC method. The displacement at the end of the time step U^j can be calculated from the mid-step velocity ${{\displaystyle \dot{U}}}^{j-\frac{1}{2}}$, known on the current time step [t_j−1, t_j] in the Central Difference scheme: $$U^j=U^{j-1}+\text{Δ}t{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}$$(3)

In the case of bridge cranes under seismic excitation, the impulse term I^j in equation (2) is represented by three Lagrange multiplier vectors ${\lambda }_N^{j+\frac{1}{2}}$, ${\lambda }_H^{j+\frac{1}{2}}$ and ${\lambda }_T^{j+\frac{1}{2}}$ to ensure, respectively, vertical impact/contact conditions, lateral impact/contact conditions and tangential friction between wheels and rails. Three Boolean constraint matrices select the appropriate degrees of freedom involved in impact/contact conditions in terms of nodal velocities: $L_N^j$ for vertical impact/contact, $L_H^j$ is for lateral impact/contact and $L_T^j$ for tangential friction. The impulse term can be rewritten as follows: $$I^j={\left(L_N^j\right)}^T{\lambda }_N^{j+\frac{1}{2}}+{\left(L_H^j\right)}^T{\lambda }_H^{j+\frac{1}{2}}+{\left(L_T^j\right)}^T{\lambda }_T^{j+\frac{1}{2}}.$$(4)

The discrete non-smooth dynamic equation is completed by the HSM conditions, written in terms of normal gap rate and Lagrange multipliers. For all pairs of nodes potentially in contact, with the index l = 1, … , N_c and N_c being the total number of pairs of nodes, we have to satisfy the following contact conditions concerning the gap and the gap rate: $$\{\begin{array}{c}\hfill \text{if}{g}_{N,j}^l>0\text{then}{\lambda }_{N,j+\frac{1}{2}}^l=0\hfill \\ \hfill \text{if}{g}_{N,j}^l\le 0\text{then}\{\begin{array}{c}\hfill {{\displaystyle \dot{g}}}_{N,j+\frac{1}{2}}^l\ge 0\hfill \\ \hfill {\lambda }_{N,j+\frac{1}{2}}^l\ge 0\hfill \\ \hfill {{\displaystyle \dot{g}}}_{N,j+\frac{1}{2}}^l{\lambda }_{N,j+\frac{1}{2}}^l=0\hfill \end{array}\hfill \end{array}$$(5)

where $g_{N,j}^l$is the gap function, ${{\displaystyle \dot{g}}}_{N,j+1/2}^l$ is the gap rate, λ_N is the contact impulse.

Then, the vertical impact/contact case is focused to explain how to compute the impulse ${\left(L_N^j\right)}^T{\lambda }_N^{j+\frac{1}{2}}$ in equation (4). The lateral impact/contact case can be dealt in the same manner as the vertical impact/contact. When the vertical contact is detected, the discrete kinematic constraint in terms of velocity is built by writing the gap rate between the nodes at the contact interface equal to zero: ${{\displaystyle \dot{g}}}_N^{j+\frac{1}{2}}=L_N^j{{\displaystyle \dot{U}}}^{j+\frac{1}{2}}=0$. In other words, there is no normal relative vertical velocity at the contact interface. The local impact/contact equation can be written as follows, after multiplying the non-smooth equation by $L_N^j{M}^{-1}$: $$H_N{\lambda }_N^{j+\frac{1}{2}}=b_N^{j+\frac{1}{2}}$$(6) $$\{\begin{array}{c}\hfill \hfill \\ \hfill \begin{array}{l}{H}_N=L_N^j{M}^{-1}{\left(L_N^j\right)}^T\\ b_N^{j+\frac{1}{2}}=-L_N^j\left({{\displaystyle \dot{U}}}^{j-\frac{1}{2}}+\text{Δ}t{M}^{-1}\left(F^{ext,j}-C{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}-K{U}^j\right)\right)\end{array}\hfill \end{array}$$(7)

where H_N is the interface operator and the right hand-side vector $b_N^{j+\frac{1}{2}}$ can be calculated with the known quantities. The obtained Lagrange multiplier should be checked to make sure that it is positive for all pairs of contact nodes, because between the nodes in contact only compression is possible. On the contrary, the Lagrange multiplier will be set to zero. It should be also noted that the Lagrange multiplier calculated by equation (6) can correspond to an impact with a jump in velocity or a contact without jump in velocity.

Regarding the friction impulse ${\left(L_T^j\right)}^T{\lambda }_T^{j+\frac{1}{2}}$ in equation (4), it is assumed that the friction is in a stick state at the beginning, meaning that the relative tangential velocity is set to zero: ${{\displaystyle \stackrel{}{v}}}_T^{j+\frac{1}{2}}=L_T^j{{\displaystyle \dot{U}}}^{j+\frac{1}{2}}=0$. The local equation to calculate the tangential friction Lagrange multiplier can be obtained in a similar way as the vertical impact/contact case: $$H_T{\lambda }_T^{j+\frac{1}{2}}=b_T^{j+\frac{1}{2}}$$(8) $$\{\begin{array}{l}{H}_T=L_T^j{M}^{-1}{\left(L_T^j\right)}^T\\ b_T^{j+\frac{1}{2}}=-L_T^j\left({{\displaystyle \dot{U}}}^{j-\frac{1}{2}}+\text{Δ}t{M}^{-1}\left(F^{ext,j}-C{{\displaystyle \dot{U}}}^{j-\frac{1}{2}}-K{U}^j\right)\right)\end{array}.$$(9)

The tangential friction Lagrange multiplier should be also checked afterwards as shown below in Equation (10), in order to make sure that it respects the Coulomb friction law, with the index l = 1, … , N_c, N_c being the total number for the pairs of contact nodes. Let us define the constraint matrix $L_{T,l}^j$ which picks out only the degrees of freedom of the pair of nodes with the index l. $$\{\begin{array}{l}\text{if }║{\lambda }_{T,l}^{j+\frac{1}{2}}║<\mu {\lambda }_{N,l}^{j+\frac{1}{2}}\text{then}{L}_{T,l}^j{{\displaystyle \dot{U}}}^{j+\frac{1}{2}}=0\left(\text{stick state}\right)\\ ║{\lambda }_{T,l}^{j+\frac{1}{2}}║\ge \mu {\lambda }_{N,l}^{j+\frac{1}{2}}\text{then}{\lambda }_{T,l}^{j+\frac{1}{2}}=-\frac{b_{T,l}}{║b_{T,l}║}\mu {\lambda }_{N,l}^{j+\frac{1}{2}}\text{and}{L}_{T,l}^j{{\displaystyle \dot{U}}}^{j+\frac{1}{2}}\ne 0\left(\text{slip state}\right)\end{array}.$$(10)

Here, µ is the Coulomb friction coefficient. As written in Equation (10), the calculated friction impulse will be compared with the critical value defined by the coefficient of friction and the normal impulse. In the case that the calculated friction impulse is smaller than the critical value, the calculated friction impulse is valid, corresponding to the stick state in which bodies in contact remain stuck. Contrarily, if the calculated friction impulse is larger than the critical value, the critical value should be imposed, indicating that there is already relative movement between the contact nodes.

2.3 HATI procedure between explicit and implicit subdomains

In the implicit subdomain, different from the explicit subdomain with impact/contact, Newmark implicit time integration (CAA) is adopted with a large time step, given by Δt₂ = [t₀,t_m], thanks to its unconditional stability. By introducing Newmark implicit scheme predictors $U_2^{0,p}=U_2^0+\text{Δ}t{{\displaystyle \dot{U}}}_2^0+\left(\frac{1}{2}-{\beta }_2\right)\text{Δ}{t}^2{{\displaystyle \ddot{U}}}_2^0$ and ${{\displaystyle \dot{U}}}_2^{0,p}={{\displaystyle \dot{U}}}_2^0+\left(1-{\gamma }_2\right)\text{Δ}t{{\displaystyle \ddot{U}}}_2^0$ in terms of displacement and velocity into the implicit subdomain equilibrium, we obtain the effective mass matrix ${{\displaystyle \tilde{M}}}_2$, which is composed of the mass matrix, the stiffness matrix and the damping matrix as given below: $${{\displaystyle \tilde{M}}}_2{{\displaystyle \stackrel{\mathrm{..}}{U}}}_2^m=F_2^{ext,m}-C_2{{\displaystyle \dot{U}}}_2^{0,p}-K_2{U}_2^{0,p}+L_2^T{\lambda }^m$$(11) $${{\displaystyle \tilde{M}}}_2=M_2+{\beta }_2\text{Δ}{t}_2^2{K}_2+{\gamma }_2\text{Δ}{t}_2{C}_2\left({\gamma }_2=0.5,{\beta }_2=0.25\right).$$(12)

To solve the explicit/implicit coupling problem following GC coupling method, the kinematic quantities of the implicit subdomain are decomposed into two parts: free quantities and linked quantities. The free quantities can be easily obtained with the internal and external forces, while the linked quantities should be calculated from the interface forces determined by the Lagrange multiplier vector λ. Consequently, the implicit subdomain equation can be divided into two parts: $${{\displaystyle \tilde{M}}}_2{{\displaystyle \stackrel{\mathrm{..}}{U}}}_2^{m,link}=L_2^T{\lambda }^m$$(13) $${{\displaystyle \tilde{M}}}_2{{\displaystyle \stackrel{\mathrm{..}}{U}}}_2^{m,free}=F_2^{ext,m}-C_2{{\displaystyle \dot{U}}}_2^{0,p}-K_2{U}_2^{0,p}.$$(14)

Taking into consideration that the velocity continuity at the interface between explicit and implicit subdomains is ensured at the fine time scale, given by the time step Δt₁ = [t_j−1,t_j]. The free quantities of the implicit subdomain need to be interpolated between the beginning and the end of the large time step. At the coupling interface, the velocities are split into free and linked velocities, leading to the following velocity continuity equation: $$L_1{{\displaystyle \dot{U}}}_1^{link,j}+L_2{{\displaystyle \dot{U}}}_2^{link,j}=-L_1{{\displaystyle \dot{U}}}_1^{free,j}-L_2{{\displaystyle \dot{U}}}_2^{free,j}.$$(15)

On the basis of GC coupling method, it has been derived that, by replacing the linked velocities by accelerations, the above velocity continuity equation leads to the following reduced-size interface problem [25,26]: $$H{\lambda }^j=b^j$$(16)

where H is the already known interface operator and the right-hand side vector b^j depends only on the free velocities computed in explicit and implicit subdomains, as defined below: $$\{\begin{array}{c}\hfill H=\frac{1}{2}\text{Δ}{t}_1{L}_1{M}_1^{-1}{L}_1^T+\frac{1}{2}\text{Δ}{t}_2{L}_2{{\displaystyle \tilde{M}}}_2^{-1}{L}_2^T\hfill \\ \hfill b^j=-L_1{{\displaystyle \dot{U}}}_1^{free,j}-L_2{{\displaystyle \dot{U}}}_2^{free,j}\hfill \end{array}.$$(17)

Finally, the Lagrange multiplier vector can be calculated as well as the linked quantities $L_1^T{\lambda }^j$ in equation (1). After that, all the quantities (displacement, acceleration and velocity) related to the explicit subdomain 1 are known at the fine time scale t_j. We can also proceed to the next fine time scale up to the last fine time step Δt₁ = [t_j−1,t_j] from which the interface force $L_2^T{\lambda }^m$ is computed and is employed to solve the implicit subdomain 2 as given in equation (1).

3 Non-linear transient dynamic analysis of the experimentally tested bridge crane

The capacity of the proposed numerical model to predict the nonlinear seismic response of the experimentally tested bridge crane is investigated on the basis of the experimental results provided by the SOCRAT international benchmark held from 2020 to 2022 [12]. A 1/5 scale mock-up of an overhead crane was tested on the AZALEE shake table operated by Sustainable Energy and Atomic Energy Commission (CEA) in Saclay, France. The mock-up is a simplified 1/5 scale model of a 22.5 m long overhead bridge crane and composed of different components: trolley, girder beams, end truck beams and runway beams, illustrated in Figure 2.

The key aspect for numerically reproducing the response of the bridge crane subjected to the shake table test is to simulate the frictional contacts between the wheels and rails. As presented before, frictional contacts will be dealt though the use of the explicit CD-Lagrange scheme in very localized areas (8 explicit patches compose the explicit subdomain), and the main part of the mesh is dealt with a classical implicit CAA scheme. It permits us to reduce the computation time and to easily take into account the Rayleigh viscous damping without any deleterious effect on the critical time step related to the explicit subdomain. In this section, we briefly present the geometric and material characteristics of the simulated 1/5 scale mock-up of the bridge crane. Next, the HATI simulation is first validated with respect to a full explicit simulation and the HATI potential for reducing the CPU time is underlined. Then, the numerical results from the HATI simulation in MATLAB are compared with the experimental results in terms of displacements, pseudo-acceleration spectra and transmitted forces at the anchorages.

Fig. 2 1/5 scale mock-up of the bridge crane on the shake table.

3.1 Finite element model of the crane mock-up

The finite element model of the bridge crane is established according to the geometric and mechanical parameters of the experimentally tested bridge crane on the shake table. The 3D mesh, based on hexahedral elements to model the runway beams, end truck beams, girder beams and the trolley, has been set up as given in Figure 3, with axes X, Y and Z, being the longitudinal direction, transverse direction and vertical direction, respectively. By taking advantage of the hybrid approach versatility, an explicit time integrator is employed to treat the small impact/contact areas in red patches, shown in Figure 3, associated with a fine time step satisfying the stability condition to simulate high frequency phenomena of multi-impacts and frictional contacts. Newmark implicit time integrator is chosen outside the contact areas, displayed in blue in Figure 3, to adopt a large time step as well as viscous damping with stiffness matrix contribution in order to filter out the high frequency spurious vibrations. As shown in the finite element model of the bridge crane, the explicit subdomain is minimized in comparison to the implicit subdomain, in order to have a better computational efficiency and limit the effect of the removal of viscous damping in the explicit subdomain.

At the beginning, it is assumed that the trolley is located in the middle of the girder beams and the girder beams are in the middle of the runway beams. As shown in Figure 4, the trolley of the bridge crane is supported by four wheels and each end truck beam of the bridge crane is supported by two wheels. The wheels 1 to 4 concern the trolley, whereas the other four wheels concern the end truck beams.

The different components of the bridge crane are assumed to be elastic. The material characteristics of the steel are defined as follows: density ρ=7800kg/m³, Young's modulus E = 2 × 10¹¹ Pa, Poisson's ratio ν= 0.3. The geometric parameters of the different beams are given in Table 1. The initial lateral gap is defined as 2mm between the wheels and the rails. At the frictional contact interface, the pairs of points in contact are the bottom points of the wheels in contact with the points on the top surface of the rails. Concerning the boundary conditions, the ends of the runway beams are clamped on the four plates above the load cells. The support conditions, which involve load cells anchored on the shake table, are simplified in the numerical model, by considering that the stiffness of the load cell is high with respect to the stiffness of the mock-up, justifying the clamped boundary conditions on the bottom surface of the end of the runway beams.

The first eigenmodes of the numerical model are illustrated in Figure 5, with the numerical frequencies compared to the experimental results: it includes a first eigenmode related to a longitudinal displacement of the girder beams (along the axis of the girder beams), a second eigenmode related to a transverse displacement and the last one related to the vertical motion. It can be observed that numerical eigenfrequencies are in good agreement with respect to the experimental frequencies identified through the application of a white noise during the shake table campaign. Nonetheless, the second eigenmode characterized by the transverse deformation is more complicated to be captured. It should be noted that clamped conditions are assumed between the wheels and the rails in order to conduct the modal analysis of the model. The following transient nonlinear computations take into account the frictional contact which may have an effect on the frequency content of the structural vibration.

Fig. 3 Three-dimensional mesh divided into an explicit subdomain, including 8 small impact/contact areas (in red), and a large implicit subdomain (in blue) outside the impact/contact areas.

Fig. 4 Sliding configuration of the wheels.

Fig. 5 The first eigenfrequencies and eigenmodes of the numerical model: (a) numerical frequency 7.54 Hz, experimental frequency 7.6 Hz (b) numerical frequency 8.59 Hz, experimental frequency 8.3 Hz (c) numerical frequency 13.5 Hz, experimental frequency 13.3 Hz.

3.2 Multi-time step explicit/implicit simulation of the bridge crane under gravity loading

In the transient seismic analysis of the experimentally tested bridge crane, the finite element model is firstly subjected to the gravity loading before the seismic excitation, the time-dependent gravity is characterized by a progressive linear increase to the nominal value between 0s and 0.1s and remains constant from 0.1s. In these first numerical computations, HATI framework is validated against full explicit computation in order to evaluate the accuracy and efficiency of the approach. Here, all the wheels are defined as fixed in the sense that the wheels can slide and jump but not roll on the rails. This corresponds to the sliding configuration, as depicted in Figure 4. The friction coefficient of the wheel-rail interaction in the sliding configuration is assumed to be equal to 0.1 for all the wheels. In multi-time step explicit/implicit (E/I) simulations, the CD-Lagrange algorithm is employed in the explicit subdomain to reproduce the impacts and frictional contacts experienced by the wheel-rail interaction during earthquake with a fine time step Δt₁ equal to 1 × 10⁻⁶ s in this numerical model on the basis of the CFL stability condition. Newmark implicit time integration (CAA scheme) is applied in the implicit subdomain without impact/contact with a large time scale Δt₂ = mΔt₁. The time step ratio m between the explicit and implicit subdomains is set to be 100 to reduce the computation time, because the implicit scheme is unconditionally stable. stiffness proportional damping matrix is only adopted in the implicit subdomain, allowing us to not influence the stable time step in the explicit subdomain. In the full explicit reference case using CD-Lagrange scheme, the introduction of stiffness proportional damping matrix will reduce significantly the critical time step [29]. The decrease of the critical time step in the Central Difference scheme can be calculated on the basis of the following equations [29]: $$\text{Δ}t\le \frac{2}{{\omega }_{max}}\left(\sqrt{1+{\xi }_{max}^2}-{\xi }_{max}\right)$$(18) $${\xi }_{max}=\frac{{\alpha }_M}{2{\omega }_{max}}+\frac{{\alpha }_K{\omega }_{max}}{2}$$(19)

where ω_max is the maximum angular frequency of the finite element model by solving the following eigenvalue problem: det(K−ω²M) = 0. For example, considering stiffness proportional damping matrix (C = α_KK) for a damping ratio 1% in the numerical model at 50 Hz will decrease 200 times the critical time step of the full explicit reference case according to equations (18) and (19). Instead, the hybrid approach enables us to keep the same critical time step in the explicit subdomain as summarized in Table 2. The two first cases correspond to full explicit computations without or with stiffness proportional damping matrix (ξ = 0.01). It can be seen that, when adopting stiffness proportional damping matrix to filter out the spurious high frequency, the time step is drastically reduced. In the case of multi-time step explicit/implicit simulation, stiffness proportional damping matrix is only adopted in the implicit subdomain, without affecting the stable time step in the explicit subdomain. Moreover, a larger time step can be employed in the implicit subdomain so as to save the computation time.

During the gravity application, an important Rayleigh damping with mass matrix contribution is applied in addition to the stiffness matrix contribution so as to reach a stabilized state. It should be noted that the mass matrix contribution of the Rayleigh damping has almost no negative influence on the stable time step of the Central Difference scheme. Indeed, as it can be seen in equation (19), the term $\frac{{\alpha }_M}{2{\omega }_{max}}$ is near to zero if ω_max is very large, which is the case for the 3D detailed model under consideration. The vertical displacement calculated by E/I simulation in the middle of one girder beam is plotted Figure 6, compared with the full explicit computation. It is shown that a stabilized state is reproduced due to the introduced Rayleigh damping. The results obtained by E/I simulation agree well the full explicit reference, validating the hybrid coupling method. Concerning the computation time, as reported in Table 3, the E/I simulation is 5 times quicker than the full explicit reference. The case with Rayleigh damping including both mass and stiffness contributions, targeting a damping ratio 1% at 50 Hz, is added to give more insight into the computational efficiency of the hybrid method. If the stiffness matrix contribution is introduced in Rayleigh damping, the HATI simulation turns out to be much quicker than the full explicit reference, reaching a computational gain more than 1000 times. It highlights the computational efficiency of the hybrid approach, especially in the case that high-frequency filtering provided by the stiffness matrix is introduced.

Fig. 6 Vertical displacement in the middle of one girder beam.

3.3 Multi-time step explicit/implicit simulation of the bridge crane under an impulsive signal

The objective of this calculation is to calibrate the friction coefficient between the wheels and the rails by adopting the previous HATI strategy, because the friction coefficient is an important parameter for the non-linear behavior of the bridge crane under earthquake. The ratio between the implicit and explicit time steps is still equal to 100. The considered signal, illustrated in Figure 7, corresponds to an impulsive signal. It is applied through the shake table only in the transverse direction Y, that is according to the axes of the runway beams, generating a slip of the end truck along the rails of the runway beams. Displacements sensors are fixed at the ends of runway beams during the experimental test to measure the relative displacements in Y direction of the four end truck beam wheels in relation to the runway beams.

In this test, all the wheels are considered as fixed, as done in the sliding configuration shown in Figure 4, in order to calibrate the friction coefficient for wheels in their fixed configuration, that is when rolling is prevented. After the gravity loading detailed in the previous test, the unidirectional impulsive signal is applied to the numerical model. Numerically, the frictional contact is dealt with the CD-Lagrange scheme in 8 small explicit patches. A friction coefficient equal to 0.1 is considered to have a first comparison between multi-time step explicit/implicit simulation and full explicit simulation. Rayleigh damping is not introduced in this comparison, because we need to have reference results provided by the full explicit calculation. As discussed in Section 3, in the case of high-frequency Rayleigh damping with stiffness matrix contribution, such a three-dimensional full explicit simulation during 5 seconds is too time-consuming.

He relative displacements of the four end truck beam wheels in relation to the runway beams in Y direction during the impulsive signal calculated by E/I simulation are shown in Figure 8, compared with the full explicit results. The displacement starts from zero, increases to a maximum (positive Y direction) and then decreases to a minimum (negative Y direction) and finally stops at this final sliding value of about 2.4 cm. It is observed that the results of E/I simulation match the reference full explicit results, which validates the HATI approach under impulsive loading in the presence of frictional contacts. Because of the symmetry of the model, it is noted that all the four end truck beam wheels have the same sliding displacement.

Secondly, in order to calibrate the friction coefficient between the wheels and the rails in the numerical model with respect to the experimental results, the influence of friction coefficient on the sliding behavior is analyzed, by using three different friction coefficients μ equal to 0.05, 0.1 and 0.2. In the numerical model, Rayleigh damping including stiffness contributions is adopted in the implicit subdomain to filter out the high-frequency oscillations, without reducing the critical time step in the explicit domain. To calibrate the Rayleigh parameter associated with the stiffness matrix, we consider a target damping ratio 5% at 50 Hz. Figure 9 shows the sliding of the four end truck beam wheels calculated by E/I simulation in Y direction using different friction coefficients, compared with the sliding measured during the experimental test. Again, due to the symmetry of our model, all the wheels have the same behavior, thus, only the numerical relative displacements of the wheel 5 obtained by E/I simulation are shown in Figure 9. On the contrary, it is noted that experimental results have an asymmetric behavior. The sliding of the wheels 5 and 6 is about 1 cm. Concerning the wheels 7 and 8, the sliding is about 3 cm. This asymmetric behavior observed during the experimental test may be due to the deterioration of the wheel-rail interaction surface in different tests, leading to the change of the friction coefficient during the wheel-rail interaction. It could be reproduced by introducing different friction coefficients depending on the left and right runway beams. Nonetheless, the friction coefficients could also change according to the different tests. As a consequence, instead of trying to define uncertain friction coefficients, a constant friction coefficient is considered for the wheel/rail frictional contact in the following simulations. From Figure 9, the value of 0.1 is adopted, providing a relative displacement close to the mean of experimental (left and right) values.

Fig. 7 Impulsive signal in Y direction in terms of acceleration and its Fourier transform.

Fig. 8 The comparison between the E/I simulation and the full explicit reference for the sliding of the four end truck beam wheels in Y direction.

Fig. 9 Comparison between numerical and experimental relative displacements of the four wheels in Y direction.

3.4 Multi-time step explicit/implicit simulation of experimental tests under seismic excitations at 0.5g

Three-directional seismic excitation is considered here. The shake table tests employed ground motion recordings made during the Chüetsu-oki earthquake in July 2007 at the Kashiwazaki-Kariwa nuclear power plant [11]: the signals were measured at a height of 22 meters in the turbine building of unit 5 of the power plant and were adapted to coincide with the main frequency of the crane mock-up in the Y direction. After a first gravity stage, previously described, the three-directional seismic excitation is applied to the finite element model of the bridge crane. The acceleration signals are plotted in Figure 10, with PGA (peak ground acceleration) equal to 0.5g in X and Y directions. It has to be noted that the three-directional input is provided by the Benchmark SOCRAT, corresponding to a mean value of the recorded accelerations at the crane supports during the shake table tests. The same sliding configuration, given in Figure 4, in which all the 8 wheels are prevented from rolling (fixed wheels), is used for this test. The friction coefficient for wheel/rail contacts is equal to 0.1. The accelerations on the observation points provided by the SOCRAT Benchmark were recorded by the accelerometers placed on the girders and runway beams, as depicted in Figures 11 and 12. A tri-axial accelerometer is installed on the middle of the upper surface of the trolley. In earthquake engineering, response spectra are widely used to provide the graphical representation of the response of a structure to a given earthquake. They describe the maximum peak response of a structure in terms of displacement, velocity or acceleration as a function of period or frequency. Here, pseudo-acceleration response spectra will be computed from accelerations recorded on the observation points with a critical damping ratio equal to 5%, classically considered for crane-bridge under earthquake loading. The pseudo-acceleration spectra are calculated by multiplying the maximum displacement spectra by the square of the angular frequency.

The pseudo-acceleration response spectra on the observation point of the girder beam computed by E/I simulation are compared with response spectra of the experimental results and the input signals in Figure 13. The first peak in Figure 13, at about 8 Hz, is close to the eigenmode in the X direction. The numerical response spectrum in the Z direction is a very good agreement with experimental results, with a peak at about 11 Hz. The reproduction of the seismic crane behavior in the transverse direction is less well reproduced. It can be noted that eigenmodes are visible on the response spectra peaks. The high-frequency branch observed in the experimental and numerical spectra, which corresponds to the frequencies higher than 50 Hz, comes from the numerous vertical impacts and frictional contacts of the wheels on rails, as discussed in Ambiel et al. [15,16], for the case of pounding between steel-concrete composite buildings.

The pseudo-acceleration response spectra computed from the accelerations recorded on the observation point of the runway beam are also compared to the results of E/I simulation, as shown in Figure 14. It is observed that experimental and numerical response spectra are in very good agreement in Y direction, because this direction corresponds to the strong axial stiffness of the runway beams, thus the recorded signal is close to the applied signal at the end of the runway beams. In X and Z directions, the numerical and experimental spectra have similar trends.

For the trolley, the numerical and experimental response spectra, plotted in Figure 15, are in satisfactory agreement in X direction. In Z direction, the seismic response predicted by the numerical model is quite good, whereas the response in Y direction of the trolley is poorly reproduced. The difference in Y direction could be explained by numerous lateral impacts/contacts between the wheels of the trolley and the rails of the girder beams, as illustrated in Figure 16. In our numerical model, the initial lateral gap is equal to 2mm between the wheels of the trolley and the rails of the girder beams in Y direction. The lateral relative displacements in Y direction of the wheels of the trolley with respect to the girder beams using explicit/implicit simulation are shown in Figure 17. As observed, the wheels 3 and 4 have the same sliding in Y direction as the wheels 1 and 2, due to the symmetry of our numerical model. In addition, as illustrated in Figure 17, when the lateral relative displacements reach the initial gap which is equal to 2mm, lateral impacts/contacts are triggered, preventing the trolley from the movement in Y direction. An asymmetric behavior was observed in the experimental campaign, with a skewing response of the trolley. This implies that the lateral gaps between different wheels and rails may also be different during the tests, which could modify the pseudo-acceleration response spectrum in Y direction. The reason why the lateral gaps are not the same in the experimental test is uncertain: the possible fabrication imperfections or the possible plastic deformations experienced by the wheels during the lateral impacts/contacts. Eventually, we could use asymmetric lateral gaps to improve the results in this test in the Y direction, but the gaps could also change according to the different tests. Rather than adopting random parameters in each test, which are not accurately identified in the experimental tests, we prefer to keep the initial symmetric gap of 2 mm in the different simulations as an average value.

Fig. 10 Three-directional seismic excitation applied to the finite element model.

Fig. 11 Observation point on the girder beam

Fig. 12 Observation point on the runway beam.

Fig. 13 Pseudo-acceleration response spectra computed from accelerations on the observation point of the girder beam.

Fig. 14 Pseudo-acceleration response spectra computed from accelerations on the observation point of the runway beam.

Fig. 15 Pseudo-acceleration response spectra computed from accelerations recorded at the upper surface center of the trolley.

Fig. 16 Vertical and lateral contacts during wheel-rail interaction.

Fig. 17 The relative displacement of the trolley wheels with respect to the girder beams in Y direction.

3.5 Multi-time step explicit/implicit simulation of experimental tests under seismic excitations at 1g

A three-dimensional seismic excitation with a PGA equal to 1g is now considered, plotted in Figure 18. For this high-level acceleration test, relative displacements of the end truck beams on the runway beams along the Y-axis are measured by wire sensors as well as displacements of the two girder beams at mid-span along the Z-axis. The configuration for the wheels is mixed, characterized by 4 sliding wheels and 4 rolling wheels, as shown in Figure 19. As in the previous computation, the friction coefficient for sliding wheels is 0.1, while a very low friction coefficient is chosen for the rolling wheels, taken as equal to 0.01. Again, after a first gravity stage, the three-directional seismic excitation is applied to the finite element model of the crane.

In Figure 20, the vertical displacements of two girder beams at mid-span obtained by E/I simulation are compared to the displacements measured during the experimental campaign. Here, the Pcharge1H1 corresponds to the observation point indicated in Figure 11 and Pcharge2H1 corresponds to the observation point located in the middle on the other girder beam. As shown, the results obtained by E/I simulation have a good agreement with experimental results. The sliding of end truck beam wheels on the rails of the runway beams, calculated by E/I simulation, is compared to the experimental results in Figure 21. Because of the symmetry of our numerical model, all the wheels have the same behavior and only the numerical relative displacement of the wheel 6 is plotted. For the experimental results, the wheel 5 has the same behavior as the wheel 6 and the wheel 7 has the same behavior as the wheel 8. Consequently, to simplify the drawing, only the relative displacements of the wheels 6 and 8 are plotted. Globally speaking, the sliding stages of the end truck beams versus time, in positive and negative Y directions, are well reproduced, although the final position of the crane is not well captured by the explicit/implicit simulation. Again, from experimental data, it is noted the asymmetric dynamic response of the crane, with higher displacement on the left runway than on the right runway. This observed trend for a skew response of the crane on the rails cannot be captured by the model which assumes a symmetric model and a constant friction coefficient for the wheel-rail interaction. As a consequence, asymmetric behavior is prevented in our numerical model, in which all sliding displacements of end truck wheels are identical.

When the loading is beyond design conditions, the behavior of the anchoring system is a primary cause of failure for bridge cranes. Consequently, it is crucial to estimate the efforts transmitted to the anchorages, which are the reaction forces at the interface between the bridge and its supporting structure. It was one of the main motivations for the Benchmark SOCRAT and this is why four load cells are installed between the end of runway beams and the shake table. The forces recorded at the load cell at the left support A shown in Figure 19 are compared with the numerical results in Figure 22. The Fourier spectra of the reaction forces in three directions are also provided in Figure 23. Despite some discrepancies, the numerical model provides a reasonable estimation for the reaction forces in the three directions. In addition, the vertical relative displacement of the end truck wheels with respect to the rails of the runway beams calculated by E/I simulation is plotted in Figure 24 in comparison to the the reaction force in Z direction. As observed, the detachment and recontact between the wheels and the rails are consistent with the peaks of the reaction force in Z direction.

Finally, the discrete energy balance is plotted in Figure 25, corresponding to the summation of the energies in the two partitions. The implicit partition includes the kinetic, internal, external energies as well as the damping energy dissipated through the viscous damping. The explicit partition includes the kinetic, internal, external energies as well as dissipated frictional energy due to the frictional sliding of the wheels on the rails. As expected, the spurious normal contact energy can be neglected in comparison to the other energies. From Figure 25, it can be highlighted that the external energy coming from the earthquake excitation is mainly dissipated through viscous damping in the material bulk through the introduction of the damping matrix, and to a lesser extent, by the friction between the wheels and rails.

Fig. 18 Three-directional seismic excitation applied to the finite element model (1 g).

Fig. 19 Mixed configuration of the wheels.

Fig. 20 Comparisons between numerical and experimental vertical displacements of the girder beams at mid-span.

Fig. 21 Sliding of end truck beam on the runway beam along the Y-axis.

Fig. 22 Reaction forces at the support A in different directions.

Fig. 23 Fourier spectra of reaction forces in the three directions.

Fig. 24 Detachment and recontact compared with the peaks of the reaction force in Z direction.

Fig. 25 Discrete energy balance for the whole model including the explicit and implicit partitions: Wext (external energy due to the earthquake excitation), Wvisc (damping energy dissipated by the damping matrix), Wcontact (spurious energy at the normal contact), Wfriction (frictional energy), Wkin (kinetic energy), Wint (internal energy) and Wbalance (check of the discrete energy balance equal to zero).

4 Conclusion

In the global framework of safety demonstration for nuclear reactors, the issue of non-linear dynamic behavior of crane bridges under earthquake needs to be considered, because bridge crane's failure could be a significant contributor in the probability of core meltdown. In order to satisfy more and more rigorous seismic design requirements, advanced non-linear dynamic numerical techniques are necessary to take into account the non-linear effects, such as the frictional contacts and impacts. In this paper, three-dimensional non-linear seismic analysis of an experimentally tested bridge crane is carried out in MATLAB environment using non-linear substructuring for non-smooth transient dynamics to investigate the capacity of the numerical model to predict the non-linear earthquake behavior of the bridge crane. The explicit CD-Lagrange scheme has been adopted in small contact patches whereas the remaining of the mesh is treated with a classical Newmark implicit time integration associated with a large time step. The different comparisons between numerical and experimental results give a satisfactory agreement in terms of response spectra, displacements, sliding of the mock-up on the runway beams and forces transmitted to the anchorage systems, demonstrating the ability of the proposed numerical model to reproduce the non-linear dynamic behavior of bridge cranes under earthquake excitations. The dynamic response of an overhead crane subjected to earthquake is mainly governed by the impacts and frictional contacts experienced between the wheels and the rails on which the wheels can roll or slide. It has been highlighted the high sensitivity of friction coefficient, making difficult the accurate prediction of the sliding stages of the crane mock-up versus time. It is expected that the proposed hybrid approach can provide an efficient method to improve the knowledge about the seismic non-linear behaviour of bridge cranes in nuclear plants. In particular, the numerical approach could be employed in future studies to explore, via parametric studies, the influence on the dynamic response of the friction coefficient, the gap between wheels and rails and the numerical damping. Finally, the influence of differential motions of the crane supports should also be further analysed on the basis of the proposed numerical model.

Funding

This research received no external funding.

Conflicts of interest

The authors declare no competing interests.

Data availability statement

All data analyzed during this study are included in the present article.

Author contribution statement

Conceptualization, S.L. and M.B.; Methodology, S.L. M.B. and A.G.; Software, S.L. and M.B.; Validation, S.L. and M.B.; Formal Analysis, S.L.; Investigation, S.L. and M.B.; Data Curation, S.L. and M.B.; Writing − Original Draft Preparation, S.L. and M.B.; Writing − Review & Editing, S.L. M.B. and A.G.; Visualization, S.L.; All authors have read and agreed to the published version of the manuscript.

References

All Tables

All Figures

	Fig. 1 Relationship between Δt1 and Δt2, m being the time step ratio between the two partitions.
In the text

	Fig. 2 1/5 scale mock-up of the bridge crane on the shake table.
In the text

	Fig. 3 Three-dimensional mesh divided into an explicit subdomain, including 8 small impact/contact areas (in red), and a large implicit subdomain (in blue) outside the impact/contact areas.
In the text

	Fig. 4 Sliding configuration of the wheels.
In the text

	Fig. 5 The first eigenfrequencies and eigenmodes of the numerical model: (a) numerical frequency 7.54 Hz, experimental frequency 7.6 Hz (b) numerical frequency 8.59 Hz, experimental frequency 8.3 Hz (c) numerical frequency 13.5 Hz, experimental frequency 13.3 Hz.
In the text

	Fig. 6 Vertical displacement in the middle of one girder beam.
In the text

	Fig. 7 Impulsive signal in Y direction in terms of acceleration and its Fourier transform.
In the text

	Fig. 8 The comparison between the E/I simulation and the full explicit reference for the sliding of the four end truck beam wheels in Y direction.
In the text

	Fig. 9 Comparison between numerical and experimental relative displacements of the four wheels in Y direction.
In the text

	Fig. 10 Three-directional seismic excitation applied to the finite element model.
In the text

	Fig. 11 Observation point on the girder beam
In the text

	Fig. 12 Observation point on the runway beam.
In the text

	Fig. 13 Pseudo-acceleration response spectra computed from accelerations on the observation point of the girder beam.
In the text

	Fig. 14 Pseudo-acceleration response spectra computed from accelerations on the observation point of the runway beam.
In the text

	Fig. 15 Pseudo-acceleration response spectra computed from accelerations recorded at the upper surface center of the trolley.
In the text

	Fig. 16 Vertical and lateral contacts during wheel-rail interaction.
In the text

	Fig. 17 The relative displacement of the trolley wheels with respect to the girder beams in Y direction.
In the text

	Fig. 18 Three-directional seismic excitation applied to the finite element model (1 g).
In the text

	Fig. 19 Mixed configuration of the wheels.
In the text

	Fig. 20 Comparisons between numerical and experimental vertical displacements of the girder beams at mid-span.
In the text

	Fig. 21 Sliding of end truck beam on the runway beam along the Y-axis.
In the text

	Fig. 22 Reaction forces at the support A in different directions.
In the text

	Fig. 23 Fourier spectra of reaction forces in the three directions.
In the text

	Fig. 24 Detachment and recontact compared with the peaks of the reaction force in Z direction.
In the text

Fig. 25 Discrete energy balance for the whole model including the explicit and implicit partitions: Wext (external energy due to the earthquake excitation), Wvisc (damping energy dissipated by the damping matrix), Wcontact (spurious energy at the normal contact), Wfriction (frictional energy), Wkin (kinetic energy), Wint (internal energy) and Wbalance (check of the discrete energy balance equal to zero).

In the text