#### Abstract

The paper mainly focuses on the study of the effects of Rayleigh damping in the simulations of FRP-concrete bonded joints, thereby proposing an approach to determine the value of its appropriate Rayleigh damping. Specifically, the element tests under Mode I and Mode II fracture modes were first carried out to investigate the effects of the mass proportional Rayleigh damping and the stiffness proportional Rayleigh damping. An FRP-concrete bonded joint is then employed to further investigate the effects of Rayleigh damping on the simulation results under Mode II fracture mode. It is shown that low-frequency vibrations are produced in the simulations of the specimens loaded by Mode I loading and could be damped by the mass proportional Rayleigh damping, while high-frequency vibrations are produced in the simulations of the specimens loaded by Mode II loading and could only be damped by the stiffness proportional Rayleigh damping. It also shows that the stiffness proportional damping is essential to damp out the oscillations in such simulations, thereby improving the convergence. In addition, the procedure proposed in this paper can lead to a proper interval for the value of the stiffness proportional Rayleigh damping, beyond which an unreasonable simulation result may be obtained.

#### 1. Introduction

External bonding of a fibre-reinforced polymer (FRP) plate to the soffit of the deteriorated reinforced concrete (RC) beams, which is referred to as “FRP-plated RC beams,” has become a widely used strengthening technique. Failure of such an FRP-plated RC beam is often governed by a debonding failure mode in the form of either concrete cover separation or intermediate crack debonding (IC debonding) [1–4]. Concrete cover separation initiates at the end of the FRP soffit plate with the formation of a horizontal crack along the height of steel tension bars, and then the horizontal crack propagates from the soffit plate end to the high moment region of the beam, eventually leading to the failure of FRP-plated RC beam [5, 6]. The horizontal crack plays a critical role in the occurrence of concrete cover separation. A recent study by Fu et al. [7] shows that an inclined FRP U-jacket installed near the soffit plate end could effectively mitigate the development of the critical horizontal crack, thus successfully suppressing concrete cover separation. Once concrete cover separation is suppressed, IC debonding commonly governs the design of FRP-plated RC beams.

IC debonding initiates at the toe of a significant flexural crack in the high moment region and then propagates along FRP-concrete bonded interface toward decreasing moment direction. IC debonding eventually occurs when the highly tensioned FRP plate leads to Mode II fracture behavior in the concrete around FRP plate [8–10]. IC debonding mechanism is complicated, especially when the FRP-concrete bonded interface is subjected to a complex condition, such as a fatigue loading condition, as the debonding process of the FRP plate can be affected by the cracking behavior of concrete [11, 12].

Single shear tests on FRP-concrete bonded joints have been widely used to investigate debonding behavior of FRP-concrete bonded interface, due to its simplicity in realization [5, 13–16]. In addition to extensive experimental and analytical studies [17], mesoscale finite element modelling, in which fine elements are used and FRP elements share nodes with adjacent concrete elements, has been employed in a number of studies [18–21]. Lu et al. (2006) [20] first proposed such a mesoscale FE modelling on FRP-concrete bonded joint to develop a bond-slip model for FRP-concrete interfaces. They believed that such a mesoscale FE modelling as a necessary complementary approach to experiment studies could provide detailed information on bond behavior of FRP-concrete bonded joints and is essential for deep insights into bond behavior of FRP-concrete bonded joints. Recently, Tao and Chen (2015) [21] developed a simple but robust FE model under the frame of concrete damage plasticity model to investigate the bond behavior of FRP-concrete bonded joint. An (2015) [22] and An et al. (2015) [18] extended Tao and Chen's (2015) model [21] for investigating the bond mechanism of FRP-concrete bonded joint under a quasi-static loads. In An (2015) [22], a numerical damping system is required in a mesoscale FE modelling to eliminate the dynamic effects arising from concrete cracking in the damaging process; otherwise it is impossible to achieve a satisfactory result under quasi-static loads.

Rayleigh damping [23] is the most widely used method to model dissipative forces in complex engineering structures under various dynamic loads. Rayleigh damping [23] is expressed as a linear combination of the mass and stiffness matrices, which are used to damp out low- and high-frequency oscillations, respectively. Such a qualitative conclusion cannot be easily used to define the Rayleigh damping coefficients in the use of modelling of FRP-concrete bonded joints under Mode II loading. Chen et al. (2015) [24] investigated the effect of Rayleigh damping on the numerical results of FRP-plated RC beams. They found that only the stiffness proportional damping coefficient that is related to the stiffness matrix is essential to damping out high-frequency oscillations, and the low-frequency oscillation is damped out automatically in an implicit dynamic approach, HHT- approach, which they used in their study [24, 25]. However, this may not be the case for simulating debonding behavior of FRP-concrete bonded joint under Mode II loading. Debonding behavior of FRP-concrete bonded joint is dominated by high-frequency responses, and damping out high-frequency responses may lead to inaccurate FE results. As a result, a detailed investigation on setting Rayleigh damping coefficients is required for good accuracy in the mesoscale modelling of FRP-concrete bonded joints.

This paper presents a detailed investigation on defining damping coefficients for mesoscale modelling of FRP-concrete bonded joints. Different values of mass and stiffness proportional damping ratios were first employed, respectively, in the element tests under Mode I and Mode II loading. Specimen III-6 in Yao et al. (2005) [5] was then simulated with different values for mass and stiffness proportional damping ratios. Based on the comparison between test results and numerical predictions on Specimen III-6 [5], an approach is proposed to define the value of stiffness proportional Rayleigh damping for mesoscale modelling of the debonding behavior of FRP-concrete bonded joints under Mode II loading.

#### 2. Finite Element Model

##### 2.1. Modelling of Concrete

The mechanical behavior of concrete in simulations was modelled with concrete damaged plasticity model [26] which requires the definitions of constitutive behavior under compression and tension and their associated damage models.

The constitutive behavior of concrete under uniaxial compression is described by the compression stress-strain relationship [27]where (MPa) is the compressive strength of concrete obtained from cylinder specimens; (MPa) is the initial tangent modulus of concrete [28]; (MPa) is the secant modulus of concrete from the origin to the peak compression stress ; is the compression strain corresponding to the peak compression stress and given as ; is the compression strain of concrete corresponding to the compression stress , and it is given asFor strain , the descending branch of the compression stress-strain relationship is described as follows:withThe constitutive behavior of concrete under uniaxial tension is therefore described using Hordijk's (1991) model, a model that is represented in terms of tension stress versus crack opening displacement as follows:in which (MPa) is the tensile stress normal to the crack direction under opening fracture mode; (MPa) is the concrete uniaxial tensile strength; and are material constants determined from the tensile tests of concrete; (mm) is the crack opening displacement; (mm) is the crack opening displacement at the complete loss of tensile stress and given aswhere (N/m) is the tensile fracture energy of concrete and is estimated through [29]in which (mm) is the maximum aggregate size. In this study, is assumed to be 20 mm if no test data is available.

The constitutive behavior under tension in (5) can also be expressed in terms of tension stress versus cracking strain without introducing mesh size sensitivity through the crack band theory proposed by Bažant and Cedolin [30, 31]. This theory relates the crack opening displacement, , to the cracking strain, , through the characteristic length, :where is the characteristic length of square element with four integration points (i.e., CPS4 in ABAQUS [26]) and is set to a value (i.e., is the side length of the square element) [32].

The damaging behavior of concrete is described through the degradation of its elastic modulus asin which the damage factor is assumed as a function of the inelastic deformation, that is, inelastic strain in compression or cracking strain in tension. The scalar stiffness degradation variable ( for compression and for tension) is related to the plastic deformation through [33]where ranges from zero (undamaged material) to one (fully damaged material); is the plastic strain; the ratio of the plastic strain with stiffness degradation to that without stiffness degradation is suggested as 0.7 [33], a value that is found to fit well with experimental data of cyclic tests in both compression and tension cases [22].

##### 2.2. Modelling of FRP

The FRP plate was simulated in this study using a 2D plane stress element available in ABAQUS [26] (i.e., CPS4), with its thickness set to 1 mm. Its Young's modulus is accordingly obtained by ensuring that the axial stiffness of FRP in the simulations is the same as that in the tests because the debonding behavior only occurs in the body of concrete around the FRP-concrete bonded joint [19, 21, 34].

##### 2.3. Dynamic Approach

In essence, concrete cracking is a dynamic progress that will induce local dynamic process in the middle of an overall static response, leading to the static approaches being ineffective for accurate solutions. It is for this reason that a dynamic approach proves to be more robust than conventional Newton-Raphson and arc-length techniques in overcoming the convergence difficulty arising from concrete cracking [24]. As a result, the implicit dynamic approach [24] is adopted herein to simulate the debonding behavior of FRP-concrete bonded joint under Mode II loads.

The motion equation for a dynamic problem is expressed aswhere , , and are the mass, damping, and stiffness matrices, respectively; , , and are the acceleration, velocity, and displacement vectors, respectively; and is the vector of external forces. Accordingly, and represent the inertia force and the damping force, respectively, while is the internal force arising from the deformation of material.

##### 2.4. Approach to Define Rayleigh Damping

Rayleigh damping [23] is defined by a damping matrix formed as a linear combination of the mass and stiffness matrices and given aswhere is the mass proportional damping coefficient, is the mass matrix, is the stiffness proportional damping coefficient, and is the stiffness matrix. For a given mode , the critical damping ratio in Rayleigh damping is expressed in terms of damping coefficients and aswhere is the circular frequency of mode ; is the coefficient to define the mass proportional Rayleigh damping ratio ; is used to define the stiffness proportional Rayleigh damping ratio . The mass proportional damping ratio damps low-frequency response whereas the stiffness proportional damping ratio damps high-frequency response [26], because a reasonable value for the mass proportional damping ratio could only be achieved with the value in the low-frequency domain while the value for the stiffness proportional damping ratio could only be achieved with the value in the high-frequency domain.

It has been known for long that the mass proportional damping ratio damps low-frequency response whereas the stiffness proportional damping ratio damps high-frequency response [26]. Despite that, there has not been a clear definition about low- and high-frequency responses, let alone obtaining the meaningful values for the mass and stiffness proportional coefficients. It is for this reason that a procedure is proposed in this study to find out a proper interval for both the mass proportional damping ratio and the stiffness proportional damping ratio . Specifically, the procedure is stated as follows:(1)Determine the fundamental circular frequency of the model with its own boundary conditions.(2)Calculate value and value for the desired damping ratio, which varies from 0.0005% to 50% with a multiple of 10.(3)Use these obtained values alone to assign Rayleigh damping ratio , respectively, in different simulations.(4)Run these simulations with different Rayleigh damping ratios , and then extract the load-displacement curves from these simulations when the simulation is accomplished.(5)Find out the curves that share similar results, whose minimum and maximum values for their corresponding damping ratio form the reasonable interval for stiffness proportional damping ratio . Any value within this interval is capable of removing the oscillations induced by dynamic effect without sacrificing accuracy in simulation results.

#### 3. Element Tests under Different Fracture Modes

Element tests with a 2D concrete block of 310 mm in length and 100 mm in height, as shown in Figure 1, were conducted to investigate the effects of Rayleigh damping on the simulation results under different fracture modes, that is, Mode I and Mode II fracture modes. The concrete block was meshed into quadrilateral elements with the side length as 10 mm. The concrete damaged plasticity model [26] was used to model the cracking behavior of concrete. In order to ensure that failure only occurs in the center layer (Figure 1(a)), the strength of concrete elements in the center layer was set at 20.7 MPa, which is 10% lower than that of other concrete elements.

**(a)**

**(b)**

##### 3.1. Element Tests under Mode I Loading

To conduct element test under Mode I loading, the aforementioned concrete block was fixed in both directions (i.e., and directions) at the node at the left upper corner, and the other nodes at the left end were fixed in the horizontal direction (i.e., direction) (Figure 1(a)), while the nodes at the right end of the block were applied with a displacement load.

The dynamic approach was employed in this element’s tests with the damping ratio varying from 0, 0.05%, 5%, to 50%. The values of mass and stiffness proportional damping coefficients are obtained from their corresponding damping ratios through (13) with the fundamental circular frequency of the concrete block 2961 rad/s, which is obtained through frequency analysis in ABAQUS.

Different values of mass proportional damping coefficient , ranging from 0, 2.96, 296.08, to 2960.83, are first employed alone corresponding to the mass proportional damping ratio of 0, 0.05%, 5%, and 50% in the element test, thereby investigating its effect on the simulation results. The stress is defined as the applied tensile force divided by the area of the cross section of the block, while the displacement is the relative displacement in the direction of node pairs in the left and right edges of the center layer. The stress-displacement curves from the numerical results are shown in Figure 2. It is found that the mass proportional damping coefficient has no effect on the stress-displacement curve, leading to the stress-displacement curve of the case with zero damping being exactly the same as that of the other cases. This element’s tests indicate that there is no need to use mass proportional damping for improvement of convergence, due to unnoticeable oscillation in simulations under Mode I loading and the inherent presence of numerical damping in HHT-*α* method [24–26]. Also, no negative effect is found on the simulation results.

Different values of stiffness proportional damping coefficient , ranging from 3.38 × 10^{−7}, 3.38 × 10^{−5}, to 3.38 × 10^{−4}, which correspond to the stiffness proportional damping ratio of 0, 0.05%, 5%, and 50%, respectively, are then employed alone as Rayleigh damping ratio in this element test, thereby investigating its effect on the results. The predicted stress-displacement responses are shown in Figure 2. It is found that the stress-displacement curves are almost the same as that in the cases with the stiffness proportional damping ratio smaller than 0.05%, whereas the peak stress increases with the increase in the stiffness proportional damping ratio when it is greater than 5%. This numerical investigation indicates that the results from mesoscale modelling of Mode I fracture would be reliable only when the stiffness proportional damping ratio ranges between 0.00% and 0.05%; otherwise inaccurate results may be achieved.

In summary, the oscillations induced in the element test under Mode I fracture mode should belong to low-frequency domain, and there is no need to use mass proportional damping for improvement of convergence because of unnoticeable oscillation in simulations under Mode I loading and the inherent presence of numerical damping in HHT- method. On the other hand, an unreasonable value for the stiffness proportional damping coefficient may lead to an inaccurate result in simulations.

##### 3.2. Element Tests under Mode II Loading

Once again the same concrete block is taken to do element tests under Mode II loading with its left and right ends fixed in the horizontal direction (i.e., ) and with its left half at the bottom side fixed in the vertical direction (i.e., ). The load was uniformly applied on the right half of the top side with a displacement control technique (Figure 1(b)), leading to the two-dimensional concrete block under pure shear loading.

Different mass and stiffness proportional Rayleigh damping ratios were used, respectively, to investigate their effects on the simulation results under Mode II loading. The damping ratio varies from 0, 0.05%, 0.5%, 5%, to 50% for both stiffness and mass proportional Rayleigh damping ratios. The values of mass and stiffness proportional damping coefficients corresponding to different damping ratios were determined through (13) with the fundamental circular frequency of the concrete block 9073 rad/s, which was obtained through a frequency analysis on this concrete block under Mode II loading using ABAQUS.

Finite element analyses of the concrete block under Mode II loading were first conducted with the mass proportional damping alone. Different values of 0, 90.7, 907.29, and 9072.92, which represent Rayleigh damping ratios of 0, 0.05%, 0.5%, 5%, and 50%, respectively, were assigned to the mass proportional damping. The averaged shear stress and relative displacement in the vertical direction node pairs in the left and right edges of the center layer of concrete are, respectively, abbreviated as* stress* and* displacement* in this section for ease of description. The stress-displacement curves from the relevant analyses are shown in Figure 3. Occurrence of Mode II fracture in the center layer of concrete led to significant oscillations, but such oscillations cannot be effectively damped out by applying the mass proportional damping ratio even when it is increased to 50%. It is indicated that the oscillations arising in simulations of concrete under Mode II loading would belong to high-frequency vibrations.

Different values of stiffness proportional damping coefficient , ranging from 1.1 × 10^{−7}, 1.1 × 10^{−6}, 1.1 × 10^{−5}, to 1.1 × 10^{−4}, are also employed alone to investigate its effect on Mode II fracture of concrete. The shear stress-displacement curves from the relevant FE analyses are shown in Figure 3. Some oscillations due to Mode II fracture of concrete are also observed in the stress-displacement curve of the cases with the stiffness proportional damping ratio lower than 0.05%, but the oscillations were effectively damped out in the cases with stiffness proportional damping ratio equal to or larger than 0.5%. The stress-displacement curves tend to be higher than their counterpart with a smaller stiffness proportional damping ratio, when it is set to a value equal to or larger than 50%. Specifically, the ultimate stress in the case with the stiffness proportional damping ratio of 50% is shown larger by 12% than its counterparts with the stiffness proportional damping ratio of 0.05%.

In the aforementioned FE investigations it is found that the mass proportional Rayleigh damping, , is not needed to damp out the oscillations due to Mode I fracture of concrete due to unnoticeable oscillation in simulations under Mode I loading and the inherent presence of numerical damping in HHT- method, and it is ineffective to damp out the oscillations as a result of Mode II fracture of concrete. On the other hand, the stiffness proportional Rayleigh damping, , is essential to regularize oscillations due to Mode II fracture of concrete, which is characterized by high-frequency oscillations. In addition, the stiffness proportional Rayleigh damping, , tends to have a significant effect on the overall response of FE analysis when it is large enough (e.g., 5%), whereas the mass proportional Rayleigh damping, , has a minor effect on the overall response of FE analysis even when it is larger than 50%.

#### 4. Mesoscale Modelling on FRP-Concrete Bonded Joint Tests

##### 4.1. General

Specimen III-6 in Yao et al. (2005) is employed as a reference case of the specimens of FRP-concrete bonded joint to further investigate the effect of stiffness proportional Rayleigh damping on the simulation results under Mode II loading. Specimen III-6 consists of a 100 (width) × 150 (height) × 350 (length) mm concrete prism bonded with a 100 (width) × 0.165 (thickness) × 100 (length) mm CFRP plate on its top surface. The compressive strength of concrete is 27.1 MPa, and the elastic modulus of FRP plate is 256 GPa. This specimen was fitted into a setup as shown in Figure 4 and loaded through applying a displacement load at the FRP plate.

A 2D model is constructed to simulate the mechanical behavior of Specimen III-6, whose concrete and FRP parts are both simulated using a 4-node plane stress element (CPS4) with full-integration scheme. They are connected with each other by sharing nodes, with the presence of adhesive layer ignored in the geometrical modelling. As per that in the corresponding physical tests, the out-of-plane thicknesses of both FRP and concrete parts for Specimen III-6 are both set to 100 mm.

As shown in Figure 5, the size of concrete part in these simulations is set as 344 × 150 mm, based on the consideration of insensibility of the length of concrete part on simulation results and the consideration of the preference for square element in simulations while considering the characteristic length of the square elements in the use of damaged plasticity model. On the other hand, the size of FRP part is set as 100 × 1 mm by assuming its thickness as 1 mm rather than its nominal thickness for the sake of convenience in the geometrical modelling.

In order to balance the computational accuracy and efforts, the concrete part is divided into five parts (i.e., A1, A2, A4, A8, and B8 zones) that are meshed with different sizes of elements (Figure 5). The A1 zone with a dimension of 136 mm × 24 mm is meshed with square elements with side length of 1 mm, because this is the region where concrete cracking possibly occurs under loading. The A8 and B8 zones, where less cracks of concrete occur, are set with the dimensions of 344 × 120 mm and 208 × 24 mm, respectively, and are both meshed with square elements with side length of 8 mm. The A2 and A4 zones are the transition zones, which connect A1 zone and A8 zone, are set with the dimensions of 344 × 2 mm and 344 mm × 4 mm, respectively, and are meshed using square elements of 2 mm and 4 mm in side length, respectively. Adjacent elements in the same zone are connected with each other through sharing nodes, while the neighbouring elements in different zones are connected to each other through the tie constraints.

The constraints in physical tests are simulated by constraining the related nodes in the vertical, horizontal or both directions, as shown in Figure 5. Specifically, except the nodes at the top-left corner that are restrained in the vertical direction to prevent uplifting behavior of the specimen while loading, all the nodes at the bottom are restrained in the vertical direction (i.e., direction) while nodes at the lower-right corner are restrained in the horizontal direction (i.e., direction) to simulate the horizontal supports in physical tests. The specimen is loaded through a displacement load at the right end of the FRP plate. In addition, the implicit dynamic approach described in the previous sections is employed in this study for the solutions.

##### 4.2. Effect of Stiffness Proportional Damping

In the simulations of Specimen III-6 in Yao et al. (2005), different values of the stiffness proportional damping ratio are taken to investigate its effect on the predicted response in the relevant simulations. The value of stiffness proportional damping ratio varies from 0.0005% to 50% at a multiple of 10. With the fundamental frequency of 9973 rad/s, which was determined through a frequency analysis, the stiffness proportional damping coefficient varies from 1 × 10^{−9} to 1 × 10^{−4} at a multiple of 10 corresponding to the desired damping ratio .

It is found that a satisfactory simulation result is achieved with the predicted failure mode consistent with that in physical tests, as shown in Figure 7, with an error of 5% (see Figure 6) between the ultimate debonding loads in the physical tests and that in the cases with stiffness proportional Rayleigh damping ratio set to 0.05%, 0.5%, and 5%. On the other hand, in Figure 6 diverged simulation results with serious oscillations are found in the cases with stiffness proportional Rayleigh damping ratio set to 0.0005% and 0.005%.

**(a)**

**(b)**

In addition, the predicted load from the case with the stiffness damping ratio of 50% is significantly higher than those in other cases. As a result, it is believed that the range between 0.05% and 5% is a reasonable interval for the stiffness proportional damping ratio to produce reliable FE results for Specimen III-6 in Yao et al. (2005).

#### 5. Discussions

The aforementioned studies indicate that the mass proportional damping coefficient does not have any negative effect on the accuracy of the simulation results but is not effective to damp out the numerical oscillations induced by Mode II loading, whose mechanical response belongs to high-frequency domain.

Significant effects of the stiffness proportional damping coefficient are found on the accuracy of simulation results, and there always exists an interval for the stiffness proportional damping coefficient to obtain an accurate simulation result. When a value used exceeds its upper limit, the simulation results will deviate from that in its corresponding physcial tests, because the undamaged elasticity matrix is employed in ABAQUS rather than the matrix that it should use in severe damage state while using damaged plasticity model to avoid negative damping in numerical analysis. Specifically, the element test under Mode II fracture mode with stiffness proportional damping ratios of 5% and 50% (see Figure 2) and the simulations of Specimen III-6 with stiffness proportional damping ratio 50% (see Figure 6) are the examples of such a case. On the other hand, when a value employed exceeds its lower limit the simulation results are observed with different degrees of oscillations in their load-displacement curves.

Another less attractive feature of stiffness proportional Rayleigh damping system [35] is that the value of the stiffness proportional damping ratio decreases with specimen deteriorating. Specifically, it decreases with the decrease of response frequency (see (13)), and frequency also decreases with the cracks propagating due to the stiffness degradation of the specimen, as indicated in where is the mass of specimen and is the stiffness of the specimen. This feature could be illustrated by the fact that, in the simulations of FRP-concrete bonded joint with no or insufficient stiffness proportional damping, the difficulty of convergence is only encountered when a certain length of crack has been formed at the joint, but not in the initial stage due to the presence of inherent numerical damping in HHT- method [24–26]. At the same time when FRP-concrete bonded joints deteriorated, the mass of the already debonded portion also becomes larger, as indicated in (14), a change that also makes the value of frequency smaller, thereby making the stiffness proportional damping ratio smaller in the process of damaging at the joint, as indicated in (13). Due to the inherent presence of the aforementioned feature, the value of the stiffness proportional damping coefficient is needed to be set to a relatively larger value, thereby satisfying the increasing demand for the value of the stiffness proportional Rayleigh damping in the process of crack propagation in the simulations, as stated in (14).

In addition, it is also found that the issues arising from dynamic loads in the element test under Mode II fracture mode are generally not as serious as that in simulations of FRP-concrete bonded joint. Specifically, in the element tests oscillations in the load-displacement curves are not observed until the loads in all the load-displacement curves dropped to a value close to zero, as shown in Figure 3, so the occurrence of such oscillations has no effect on the value of ultimate load in simulations of traditional structural components even with insufficient stiffness proportional Rayleigh damping ratio. On the other hand, in simulations of FRP-concrete bonded joint oscillations are observed when the loading is increasing (see Figure 6) so that the ultimate load may be influenced by the oscillations, thereby leading to a possible failure of relevant simulations. That is because the loads applied onto the FRP end are not resisted by the whole part of FRP-concrete bonded joint, but only by part of the joint, which is called effective length [36, 37]. It is for this reason that the debonding behavior of FRP-concrete bonded joint, in essence, is a series of localized debonding phenomena that lead to a plateau in the load-displacement curves, thereby making the oscillations observed in the area close to ultimate load and affecting the accuracy of the ultimate load. That is also the reason that the issues arising from underdamping of high-frequency response in simulations of FRP-concrete bonded joint are more serious than that in element tests under Mode II fracture mode.

In summary, stiffness proportional damping coefficient plays a delicate but pivotal role in damping out high-frequency oscillations arising in simulations loaded under Mode II fracture mode, especially in simulations like FRP-concrete bonded joint. It is for this reason that a reasonable interval for stiffness proportional damping ratio is needed in the relevant simulations to avoid the aforementioned issues induced by both over- and underdamping of high-frequency response.

#### 6. Conclusions

This paper has been concerned with the effects of Rayleigh damping on the simulation results under different loading modes, that is, Mode I and Mode II loading modes, especially on the simulation results of FRP-concrete bonded joint. Simple element tests, respectively, under Mode I or Mode II loading are first conducted with different values of mass and stiffness proportional ratios to investigate their effects on the simulation results under different loading modes. It is found that the frequency of vibration in the simulations of the specimen under dynamic loads is highly connected to the loading modes which are more easily to be distinguished than the frequency of vibration. It is for this reason that it could be concluded that the mass proportional damping ratio should be used to damp out the oscillations in the simulations of specimens loaded under Mode I loading mode while the stiffness proportional damping ratio should be used to damp out the oscillations in the simulations of specimens loaded under Mode II loading mode. Furthermore, it is also found that an inappropriate value for the stiffness proportional Rayleigh damping may lead to an unreasonable result in simulations, while the mass proportional Rayleigh damping may not.

Specimen III-6 conducted by Yao et al. [5] is then simulated using a mesoscale FE model to further investigate the effect of the mass and stiffness proportional Rayleigh damping ratios in simulations under Mode II fracture mode. It is found that the stiffness proportional Rayleigh damping ratio within a specified range, which is obtained through the approach suggested in Section 2.4, is essential to remove the oscillations without satisficing the accuracy of simulation results. Specifically, when the damping ratio is set to a value larger than its upper limit, the load in the simulation results is shown larger than that it should be in physical tests. On the other hand, when the damping ratio is set to a value lower than its lower limit, it is ineffective to remove the oscillations arising in the simulations, the same as that in the case with the mass proportional damping ratio alone.

#### Conflicts of Interest

The authors declare no competing financial interests.

#### Acknowledgments

The authors would like to acknowledge financial support received from Science Foundation of China University of Petroleum, Beijing (Project nos. 2462015YJRC026 and C201601) and also from China Scholarships Council/University of Edinburgh Scholarships for the first author’s PhD study. In addition, this work is accomplished under direction of Professor Jian-Fei Chen from Queen’s University at Belfast and Professor Pankaj Pankaj from the University of Edinburgh, so the authors are also grateful to their kind help and directions.