Mechanism research of slip effect between frictional laminated beams

Abstract

For laminated beams connected by Coulomb friction, interlayer slippage occurs when interfacial shear stress exceeds the resistant friction stress. Then, the physical properties of the laminated beams will change and may even cause structural damage. In this article, the law of interlayer slippage of laminated beams is obtained by mechanical derivation, and the finite element model is used for comparison verification. First, the internal shear force calculation formula of the laminated beam considering the interlaminar friction is derived from the segment micro-element method. Second, interlayer slippage laws of the frictional laminated beams in both horizontal and longitudinal direction are derived according to the hierarchical slip determination conditions. Third, according to the state quantity of different boundary conditions, the transfer matrix method is used to solve the longitudinal length of interlayer slip. Then, the design of the algorithm program is completed by MATLAB. Finally, based on the comparison between the finite element model calculation results and the calculation results of the algorithm program, the accuracy of the proposed method is verified. The analysis results indicate that slippage destroys the integrity of the laminated beam section and reduces the bending moment of inertia of the slip section of the beam. The influence of slip effect on the frictional laminated beams should be considered in deflection calculation.

Introduction

For bridge structures, the main part is the beam. The beam can be considered as a laminated structure with myriad of fine fiber layers. The bending mechanism of beam is equal to the joint shear resistance of each fiber layer. When the interlayer shear force of Coulomb’s frictional laminated beam exceeds the maximum friction resistance between layers, the interlayer slippage will occur. Taking this into consideration, the slippage phenomenon of laminated structures in practical engineering can be explained. As slippage phenomenon is one key issue for bridge structures, many researchers in the field of civil engineering focus on this problem.

In order to study the interlayer slippage phenomenon of the laminated beam, the researchers did lots of work from two perspectives, one is the theory about beams and the other is the simulation model about beams.

For the theory about beams, some researchers introduced theory to the issue, and some are developed currently theory. Murakami ¹ presented a Timoshenko beam theory with built-in interlayer slip, in which he recovered Newmark’s Euler–Bernoulli beam theory ² by introducing the kinematic constraint of zero shear strain and provided a failure modes model of laminated beams with delamination. Hjiaj et al. ³ were concerned about the large deformation analysis of steel–concrete composite beams and divided the composite beams into the integral part and the sliding part to study the influence of shear flexibility on the overall nonlinear behavior of the composite beams by the co-rotational framework. Aydogdu ⁴ proposed a new theory of shear deformation for the analysis of static and dynamic behavior of composite beams and compared numerical models under different boundary conditions. Monetto ⁵ analyzed the mechanical behavior of a three-layer beam in a way of establishing beam model by simulating each layer of the composite beam into a linear elastic Euler–Bernoulli beam and setting the connection stiffness in the transverse direction much more stronger than in longitudinal direction. Based on layerwise displacement theory, Liu and Jeffers ⁶ proposed a multilayered shell formulation within the framework of isogeometric analysis. Through numerical case, the accuracy of the proposed formulation was proved. When the beam acts as shell or plate, the Kirchhoff theory was introduced. To analyze the Kirchhoff plate, the rational Bezier triangles in isogeometric analysis coupled with a feature-preserving automatic meshing algorithm was proposed by Liu and Jeffers. ⁷

For the simulation models about beams, there are also many scientists who have put great effort to solve different kinds of issues about the interlayer slip. Karama et al. ⁸ deduced an exponential numerical model of the composite beam considering the transverse shear stress by the principle of virtual work and analyzed the mechanical behavior of the composite beam according to the model. Krawczyk et al. ⁹ deduced a geometric nonlinear finite element analysis model of a laminated beam undergoing large deformation, subjected to small strain and moderate interlayer slip by assuming that there was only local interlayer interaction in the laminated beam and proposed that the interface shear stress could be used to define the layer interactions. Queiroz et al. ¹⁰ have established an ANSYS finite element model (FEM) to analyze and study the longitudinal slip of steel–concrete composite beams under local interaction and full interaction. Jiang et al. ¹¹ proposed a method to calculate the distribution of the section moment and the shear force of a steel–concrete composite beam based on global–local bending model. Yang et al. ¹² simplified the calculation formula for the bending stress and interfacial shear force of the steel–concrete composite beam with interlayer slippage by introducing the stress increasing factor η. Later, in order to capture the plasticity growth in large-deformation frames, a formulation which is used to utilize a layer-based discretization of a distributed plasticity isogeometric analysis frame model was proposed. Through several numerical cases, the performance of the proposed model was assessed. ¹³

Even though the theory and models are proposed to deal with the issue about the interlayer slip, many of the researchers focus on the effects of interlayer slip on the overall dynamic properties of the structure or on the mechanical properties of the slip component. Less literature discusses the law of interlayer slip of multilayered beams.

This article focuses on the law of interlayer slip between the laminated beam layers and the derivation of related quantitative analysis. First, in order to study the mechanism of slip effect, the shear force in contact layer is derived based on the section analysis method. Then, the horizontal interlayer slippage rule according to Coulomb’s law of friction was presented. Third, based on the transfer matrix method (TMM) and the Laplace transfer method, the program with MATLAB used to calculate the length of longitudinal slippage between frictional laminated beams was proposed. Finally, the ANSYS numerical simulation results are compared to verify the accuracy and efficiency of the proposed program.

The shear force in contact layer

Take the bending segment of the laminated beam with N layers as the research object, as shown in Figure 1. The cross sections of the frictional laminated beams with variable width and equal height are shown in Figure 2.

Figure 1.

Bending segment of the laminated beam with N layers.

Figure 2.

Cross sections of the laminated beam with variable width and equal height.

When the laminated beam is subjected to forces as a whole structure without interlayer slippage, the shear stress distribution of the segment is uniform. According to theorem of conjugate shearing stress, ¹⁴ the shear stress in the beam is

$$\tau =\frac{F_S{S}_z^{*}}{I_{z}b}$$

(1)

where $F_s$ is the shear force acting on the cross section. $S_z^{*}$ is the static moment of the area $A^{*}$ to the neutral axis. $A^{*}$ is the area portion with a distance y from the neutral axis on the cross section. $I_Z$ is the moment of inertia. $b$ is the width of the cross section

$$S_z^{*}=\underset{A^{*}}{\int }\mathrm{ydA}$$

(2)

For the frictional laminated beams with interlayer slippage, the effect of the interlayer friction between layers should be considered. Assuming that interlayer slippage has occurred in the laminated beam, the cross section of the beam with a height of $h^P$ remains intact, and the frictional force acting on the adjacent surface of the upper and lower surface is $f_t$ and $f_b$ , respectively, as shown in Figure 1. Taking infinitesimal element of length $\mathrm{dx}$ of the laminated beam segment as the research object, ¹⁵ the free-body diagram is shown in Figure 3.

Figure 3.

Free-body diagram of one infinitesimal element of the laminated beam.

According to the Euler–Bernoulli beam theory, the equilibrium of the infinitesimal element of the laminated beam in the longitudinal direction conforms to the following equations

$${\mathrm{dF}}_N^P=(f_t-f_b)\mathrm{dx}$$

(3)

$$d{M}^P=F_S^P\mathrm{dx}-(f_t+f_b)\frac{h^P}{2}\mathrm{dx}=E^P{I}^{P}d{\kappa }^P$$

(4)

where $E^P$ is the elastic modulus of the laminated beam element. $I^P$ is the moment of inertia of the beam segment relative to the neutral axis. ${\kappa }^P$ is bending curvature of the neutral axis of the laminated beam.

For the volume element of the shaded part of the upper laminated beam element, as shown in Figure 4, the equilibrium in the longitudinal direction requires that

$${\mathrm{dF}}_N^{*}=(f_t-{\tau }_t^{*}{b}^{*})\mathrm{dx}$$

(5)

where $b^{*}$ is the width of the cross section with a y* distance from the neutral axis.

Figure 4.

The development process and rule of slippage in the laminated beam with N layers: (a) state 1: without slippage, (b) state 2: first slippage occurs, (c) state 3: slippage develops, and (d) state 4: slippages occur in all contact layers.

$F_N^{*}$ can be deduced from the integral of axial stress on the section of the laminated beam element

$$F_N^{*}=\int {\sigma }^{*}\mathrm{dA}=\underset{A^{*}}{\int }(-E^P{\kappa }^{P}y+\frac{F_N^P}{A})\mathrm{dA}$$

(6)

where ${\sigma }^{*}$ is the axial stress of the cross section with y* distance from the neutral axis.

Equation (6) is equal to the differential form as follows

$${\mathrm{dF}}_N^{*}=\underset{A^{*}}{\int }(-E^P(d{\kappa }^P)y+\frac{{\mathrm{dF}}_N^P}{A})\mathrm{dA}$$

(7)

Combining equations (5) and (7)

$$\begin{array}{c}\multicolumn{1}{c}{{\tau }_t^{*}{b}^{*}-f_t=\underset{A^{*}}{\int }(E^{P}y\frac{d{\kappa }^P}{\mathrm{dx}}-\frac{1}{A}\frac{{\mathrm{dF}}_N^P}{\mathrm{dx}})\mathrm{dA}}\\ \multicolumn{1}{c}{=E^P\frac{d{\kappa }^P}{\mathrm{dx}}\underset{A^{*}}{\int }\mathrm{ydA}-\frac{A^{*}}{A}\frac{{\mathrm{dF}}_N^P}{\mathrm{dx}}}\end{array}$$

(8)

By substituting equations (2) and (3) into equation (8), we deduce the shear stress ${\tau }_t^{*}$ of the upper laminated beam segment at a y* distance from the neutral axis, as shown in equation (9)

$${\tau }_t^{*}=\frac{1}{b^{*}}[E^P{S}^{*}\frac{d{\kappa }^P}{\mathrm{dx}}-(f_t-f_b)\frac{A^{*}}{A}+f_t]$$

(9)

where $A^{*}$ is the cross-sectional area of the volume element, that is, the shaded part of the beam element shown in Figure 3. $A$ is the area of the total cross section.

Similarly, the shear stress ${\tau }_b^{*}$ of the bottom laminated beam segment at a y* distance from the neutral axis can be deduced, as shown in equation (10)

$${\tau }_b^{*}=\frac{1}{b^{*}}[E^P{S}^{*}\frac{d{\kappa }^P}{\mathrm{dx}}+(f_t-f_b)\frac{A^{*}}{A}+f_b]$$

(10)

The relative slippage of frictional laminated beam with variable width and equal height is shown in Figure 2. The relative slippage occurs mainly between layers. To accurately illustrate the position of slippage, the definition of contact layer is introduced.

The definition of the ith contact layer represents the contact surface between the ith layer beam and the (i + 1)th layer beam. The static moment and the cross-sectional area in the contact layer can be calculated through following equations

$$S_i^{*}=\frac{h^2}{2}\sum _{i=1}^i(N-2i+1)b_i$$

(11)

$${A_i}^{*}=\sum _{i=1}^{i}h\times b_i;A=\sum _{i=1}^{N}h\times b_i$$

(12)

When $i\le N/2$ , the shear stress of the ith contact layer is

$$\begin{array}{cc}\multicolumn{1}{c}{{\tau }_t^i}& =\frac{1}{b_i}[E^P\frac{h^2}{2}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(N-2i+1)b_i-(f_t-f_b)\frac{\sum _{i=1}^{i}h\times b_i}{\sum _{i=1}^{N}h\times b_i}+f_t]\hfill \\ \multicolumn{1}{c}{}& =\frac{E^P{h}^2}{2b_i}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(N-2i+1)b_i-\frac{(f_t-f_b)}{b_i}\frac{\sum _{i=1}^i{b}_i}{\sum _{i=1}^N{b}_i}+\frac{f_t}{b_i}\hfill \end{array}$$

(13)

When $i>N/2$ , the shear stress of the ith contact layer is

$$\begin{array}{cc}\multicolumn{1}{c}{{\tau }_b^i}& =\frac{1}{b_i}[E^P\frac{h^2}{2}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(N-2i+1)b_i+(f_t-f_b)\frac{\sum _{i=1}^{i}h\times b_i}{\sum _{i=1}^{N}h\times b_i}+f_b]\hfill \\ \multicolumn{1}{c}{}& =\frac{E^P{h}^2}{2b_i}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(N-2i+1)b_i+\frac{(f_t-f_b)}{b_i}\frac{\sum _{i=1}^i{b}_i}{\sum _{i=1}^N{b}_i}+\frac{f_b}{b_i}\hfill \end{array}$$

(14)

Horizontal interlayer slippage rule

According to Coulomb’s law of friction, ¹⁶ the maximum static friction stress is equal to the kinetic friction stress, and the kinetic friction stress is equal to the product of the normal contact stress $p$ and the friction coefficient $\mu$ . The ultimate friction resistance of the ith contact layer $f_{\mathrm{cr}}(i)$ is equal to the maximum static friction stress

$$f_{\mathrm{cr}}(i)=\{\begin{array}{c}\frac{\mu V}{L{b}_i},1\le i\le \frac{N}{2}\hfill \\ \multicolumn{1}{c}{\frac{\mu V}{L{b}_{i+1}},\frac{N}{2}<i\le N-1}\end{array}$$

(15)

where V is the applied external load and L is the length of the laminated beam. The interlayer slip state discrimination condition satisfies the following relationship ¹⁷

$$\{\begin{array}{cc}\tau (i)<f_{cr}(i)& (non-slipstate)\\ \tau (i)=f_{cr}(i)& (criticalstate)\\ \tau (i)>f_{cr}(i)& (slipstate)\end{array}$$

(16)

State 1: when applying small external load, the laminated beam remains sectionally intact without interlayer slippage. So, the frictional forces acting on the adjacent surface of the upper and lower surface are equal to zero, that is, $f_t=f_b=0$ . The shear stress in the section can be calculated through equation (17)

$$\{\begin{array}{c}{\tau }_0^t=\frac{E^P{h}^2}{2b_i}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(N-2i+1)b_i,1\le i\le \frac{N}{2}\hfill \\ \multicolumn{1}{c}{{\tau }_0^b=\frac{E^P{h}^2}{2b_{i+1}}\frac{d\kappa }{\mathrm{dx}}\sum _{i=i+1}^N(2i+1-N)b_i,\frac{N}{2}<i\le N-1}\end{array}$$

(17)

The distribution law of the shear stress and the ultimate friction resistance in non-slip state is shown in Figure 4(a). It can be seen that the shear stress is distributed on the cross section quadratically. Therefore, the interlayer slippage occurs first in the contact layer at the neutral axis.

State 2: the first slippage in the contact layer at the neutral axis occurs as the applied external load increases. Then, the laminated beam section can be divided into top part and bottom part, respectively. Each part can be thought as an independent frictional laminated beam. The shear stress at state 2 can be calculated by equation (18)

$$\{\begin{array}{c}{\tau }_{1-\mathrm{top}}^t=\frac{E^P{h}^2}{2b_i}\frac{d\kappa }{\mathrm{dx}}\sum _{i=1}^i(\frac{N}{2}-2i+1)b_i+\frac{\mu V}{L{b}_i}\frac{\sum _{i=1}^i{b}_i}{\sum _{i=1}^{\frac{N}{2}}{b}_i}\hfill \\ \multicolumn{1}{c}{(\mathrm{while}1\le i\le \frac{N}{4})}\\ \multicolumn{1}{c}{{\tau }_{1-\mathrm{top}}^b=\frac{E^P{h}^2}{2b_i}\frac{d\kappa }{\mathrm{dx}}\sum _{i=i+1}^{\frac{N}{2}}(2i+1-\frac{N}{2})b_i+\frac{\mu V}{L{b}_i}\frac{\sum _{i=1}^i{b}_i}{\sum _{i=1}^{\frac{N}{2}}{b}_i}}\\ \multicolumn{1}{c}{(\mathrm{while}\frac{N}{4}\le i\le \frac{N}{2})}\end{array}$$

(18a)

$$\{\begin{array}{c}{\tau }_{1-\mathrm{bottom}}^t=\frac{E^P{h}^2}{2b_{i+1}}\frac{d\kappa }{\mathrm{dx}}\sum _{i=\frac{(N+2)}{2}}^i(\frac{N}{2}-2i+1)b_i+\frac{\mu V}{L{b}_{i+1}}\frac{\sum _{i=i}^N{b}_i}{\sum _{i=\frac{(N+2)}{2}}^N{b}_i},\hfill \\ \multicolumn{1}{c}{(\mathrm{while}\frac{N}{2}<i\le \frac{3N}{4})}\\ \multicolumn{1}{c}{{\tau }_{1-\mathrm{bottom}}^b=\frac{E^P{h}^2}{2b_{i+1}}\frac{d\kappa }{\mathrm{dx}}\sum _{i=i+1}^N(2i+1-\frac{N}{2})b_{i+1}+\frac{\mu V}{L{b}_{i+1}}\frac{\sum _{i=i}^N{b}_i}{\sum _{i=\frac{(N+2)}{2}}^N{b}_i},}\\ \multicolumn{1}{c}{(\mathrm{while}\frac{3N}{4}<i\le N-1)}\end{array}$$

(18b)

The distribution law of the shear stress and the ultimate friction resistance in state 2 after the first slippage is shown in Figure 4(b). It can be seen that after the first slippage, the shear stress presents two symmetrical quadratic distributions. The peak values are located at the upper and lower layers which are close to the overall neutral axis. Therefore, it can be inferred that the second slippage appears symmetrically at the (i–1)th and (i + 1)th contact layers.

State 3: as the applied external force continues to increase, the laminated beam can be divided into three parts when the second slippage occurs, as shown in Figure 4(c). Among the three parts, the top and bottom segments can be seen as single-sided laminated beam subjected to the frictional force, while the middle segment is the double-sided laminated beam subjected to the frictional force. For the single-sided laminated beam subjected to the frictional force, the shear stress can be calculated by equation (18). For the double-sided laminated beam subjected to the frictional force, the shear force can be calculated by equation (19)

$${\tau }_2^i=\frac{E^P{h}^2}{4}\frac{d\kappa }{\mathrm{dx}}(\frac{h^2}{4}-y_i^2)+\frac{\mu V}{L{b}_i}$$

(19)

where $y_i$ is the distance from the calculating shear stress place to the neutral axis of the independent section.

State 4: from the previous derivation, it can be concluded that the laminated beam will repeat State 3 until all the contact layers appear to slip and become N layers of independent beam without bond, as shown in Figure 4(d). The shear stress is distributed in a zigzag pattern across the cross section and satisfied continuity.

Longitudinal interlayer slippage rule

The relationship between the interlayer slippage and the change ratio of curvature

• 1. When the first slip occurs, the change ratio of curvature in state 1 can be deduced by combining equations (15) and (17)

$${\frac{d\kappa }{\mathrm{dx}}|}_{\mathrm{slip}-1}=\frac{2\mu V}{E^P{h}^{2}L\sum _{i=1}^{N/2}(N-2i+1)b_i}$$

(20)

• 2. When the second slip occurs, the change ratio of curvature in state 2 can be deduced by combining equations (15) and (18)

$${\frac{d\kappa }{\mathrm{dx}}|}_{\mathrm{slip}-2}=\frac{4\mu V}{E^P{h}^{2}L(N-2)\sum _{i=1}^{N/2}{b}_i}$$

(21)

• 3. The relationship between the change ratio of curvature and the interlayer slippage can be summarized according to induction deduction method

$${\frac{d\kappa }{\mathrm{dx}}|}_{\mathrm{slip}-k}=\{\begin{array}{c}\frac{2\mu V}{E^P{h}^{2}L\sum _{i=1}^{N/2}(N-2i+1)b_i},k=1\hfill \\ \multicolumn{1}{c}{\frac{4\mu V}{E^P{h}^{2}L(N+2-2k)\sum _{i=1}^{(N+4-2k)/2}{b}_i},2\le k\le \frac{N}{2}}\end{array}$$

(22)

So, the shear force at kth slippage state can be calculated based on the force balance in infinitesimal segment

$$\begin{array}{cc}\multicolumn{1}{c}{F_S}& =\sum _{n=1}^{2k}{F}_S^{n-P}\hfill \\ \multicolumn{1}{c}{}& ={E{I}_{k-\mathrm{all}}\frac{d\kappa }{\mathrm{dx}}|}_{\mathrm{slip}-k}+\frac{\mu \mathrm{Vh}}{2L}[N+2(k-1)]\hfill \end{array}$$

(23)

where $I_{k-\mathrm{all}}$ is a sum of the moment of inertia of all beam segments to their own independent neutral axis at kth slippage state. It can be calculated by equation (24)

$$I_{k-\mathrm{all}}=\sum _{n=1}^{2k}{I}_n^P=\frac{h^3}{12}(\frac{\left(\sum _{i=1}^m{b}_i\sum _{i=1}^m{b}_i(3{\beta }^2+1)\right)-3{\left(\sum _{i=1}^m(2i-1)b_i\right)}^2}{\sum _{i=1}^m{b}_i}+2\sum _{k=1}^{k-1}{b}_m)$$

(24)

where $m=(N+4-2k)/2$ and $\beta =N+5-2k-2i$ .

Calculation of interlayer slip length based on TMM

The longitudinal length of interlayer slippage can be calculated by equations (22)–(24). But, the longitudinal length calculation is related to the boundary conditions of the laminated beam. To establish the solvable function between the whole structure and the boundary conditions, the TMM is introduced.

TMM for beam element considering axial force

Assuming an arbitrarily selected i-unit of the laminated beam structure, the left-end state vector is S_i–1 and the right-end state vector is S_i. According to the mechanical equilibrium of the unit and assuming the unit and node state vector deliver continuously, the relationship between the leftmost node and the rightmost node in the whole structure can be established^18–21

$$S_i=\underset{i=1}{\overset{k_n}{\Pi }}{U}_i{U}_P{S}_0$$

(25)

where U_i is field matrix, U_P is point matrix, and k_n is the total number of units.

The overall stiffness matrix of TMM is the multiplicative form of field matrix and point matrix. This method does not change the dimension of matrix in the process of multiplication, saves computer memory, and is easy to solve. It is especially suitable for solving the mechanics problems of chain structure.

If the beam is subjected to the axial force in addition to the vertical load, the interaction between the axial force and the deflection will inevitably cause the moment balance expression function to produce an additional term.^22,23 Ignore the longitudinal displacement of the microbodies and assume that the axial force N acting on the end of the beam is a constant value. The mechanical equilibrium equation of the beam element subjected to the axial force can be deduced.

Select an infinitesimal element of a beam with the length of $l_{(i)}$ , as shown in Figure 5. Define the state variables including the vertical deflection $\omega$ at both ends of the element, the angle of rotation $\theta$ , the moment $M$ , and the shear force $Q$ . Then, the state vector at both ends of the beam element can be expressed as

$$\begin{array}{c}{S}_{0(i)}={\left\{\begin{array}{cccc}{\omega }_0\hfill & {\theta }_0\hfill & M_0\hfill & \begin{array}{cc}{Q}_0\hfill & 1\hfill \end{array}\hfill \end{array}\right\}}_{(i)}^T\hfill \\ S_{l(i)}={\left\{\begin{array}{cccc}{\omega }_l\hfill & {\theta }_l\hfill & M_l\hfill & \begin{array}{cc}{Q}_l\hfill & 1\hfill \end{array}\hfill \end{array}\right\}}_{(i)}^T\hfill \end{array}$$

Figure 5.

Free-body diagram of an infinitesimal element of beam under axial force.

According to the different boundary conditions at both ends of the laminated beam, including the common methods of free end, simple support end, and fixed end, the state quantity of the boundary condition selects different combinations of values, as shown in Table 1.

Table 1.

State quantity of different boundary conditions.

Free end	Simple supported end	Fixed end
$\begin{array}{c}M=0,\theta \ne 0\hfill \\ Q=0,w\ne 0\hfill \end{array}$	$\begin{array}{c}w=0,\theta \ne 0\hfill \\ M=0,Q\ne 0\hfill \end{array}$	$\begin{array}{c}\theta =0,M\ne 0\hfill \\ w=0,Q\ne 0\hfill \end{array}$

Ignoring the effect of lateral motion on longitudinal motion, ²⁴ the following differential equations can be obtained according to the stress equilibrium of the infinitesimal beam element

$$\{\begin{array}{c}\mathrm{EI}\frac{d^4{\omega }_{l(i)}}{d{l_{(i)}}^4}+N\frac{d^2{\omega }_{l(i)}}{d{l_{(i)}}^2}=0\hfill \\ \multicolumn{1}{c}{{\theta }_{l_{(i)}}=-\frac{d{\omega }_{l(i)}}{d{l}_{(i)}}}\\ \multicolumn{1}{c}{M_{l_{(i)}}=-\mathrm{EI}\frac{d^2{\omega }_{l(i)}}{d{l_{(i)}}^2}}\\ \multicolumn{1}{c}{Q_{l_{(i)}}=-\frac{\mathrm{dM}(x)}{d{l}_{(i)}}=-\mathrm{EI}\frac{d^3{\omega }_{l(i)}}{d{l_{(i)}}^3}}\end{array}$$

(26)

Let ${\lambda }^2=\mathrm{EI}/N$ , then equation (26) can be transformed to following expression

$${\lambda }^2{\omega }_{l_{(i)}}^{(4)}+{\omega }_{l_{(i)}}^{(2)}=0$$

(27)

Laplace transform is applied to equation (27). Then, the relationship between right-end deflection and left-end deflection of the infinitesimal beam element can be calculated

$$\begin{array}{cc}\multicolumn{1}{c}{{\omega }_{l_{(i)}}=}& {\omega }_{0(i)}+l_{(i)}\omega {\prime }_{0(i)}+(1-\cos \frac{l_{(i)}}{\lambda }){\lambda }^2{\omega ″}_{0(i)}\hfill \\ \multicolumn{1}{c}{}& +(l_{(i)}-\lambda \sin \frac{l_{(i)}}{\lambda }){\lambda }^2\omega ″{\prime }_{0(i)}\hfill \end{array}$$

(28)

Combining equations (26) and (28), the deduced relationship between the vertical deflection and other state variables can be derived from equation (29)

$$\begin{array}{cc}\multicolumn{1}{c}{{\omega }_{l_{(i)}}=}& {\omega }_{0(i)}-l_{(i)}{\theta }_{0(i)}-(1-\cos \frac{l_{(i)}}{\lambda })\frac{M_{0(i)}}{N}\hfill \\ \multicolumn{1}{c}{}& -(l_{(i)}-\lambda \sin \frac{l_{(i)}}{\lambda })\frac{Q_{0(i)}}{N}\hfill \end{array}$$

(29)

Similarly, for the angle of rotation, the moment and shear force can be derived from equation (30).

$$\{\begin{array}{c}{\theta }_{l_{(i)}}={\theta }_{0(i)}+\frac{M_{0(i)}}{\mathrm{EI}\lambda }\sin \frac{l_{(i)}}{\lambda }+(1-\cos \frac{l_{(i)}}{\lambda })\frac{Q_{0(i)}}{N}\hfill \\ M_{l(i)}=\cos \frac{l_{(i)}}{\lambda }{M}_{0(i)}+\lambda \sin \frac{l_{(i)}}{\lambda }{Q}_{0(i)}\hfill \\ Q_{l(i)}=-\frac{1}{\lambda }\sin \frac{l_{(i)}}{\lambda }{M}_{0(i)}+\cos \frac{l_{(i)}}{\lambda }{Q}_{0(i)}\hfill \end{array}$$

(30)

And then, the relationship between the state vectors at both ends of the beam element can be expressed as a matrix form with the transform matrix

$$\begin{array}{c}{S}_{l_{(i)}}={\left\{\begin{array}{c}{\omega }_l\hfill \\ \multicolumn{1}{c}{{\theta }_l}\\ \multicolumn{1}{c}{M_l}\\ \multicolumn{1}{c}{Q_l}\\ \multicolumn{1}{c}{1}\end{array}\right\}}_{(i)}={\left\{\begin{array}{ccccc}1\hfill & -l\hfill & U_{13}\hfill & U_{14}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & U_{23}\hfill & U_{24}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & U_{33}\hfill & U_{34}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & U_{43}\hfill & U_{44}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right\}}_{(i)}\hfill \\ {\left\{\begin{array}{c}{\omega }_0\hfill \\ \multicolumn{1}{c}{{\theta }_0}\\ \multicolumn{1}{c}{M_0}\\ \multicolumn{1}{c}{Q_0}\\ \multicolumn{1}{c}{1}\end{array}\right\}}_{(i)}=U_i{S}_{0_{(i)}}\hfill \end{array}$$

(31)

where

$$\begin{array}{c}{U}_{13}=-(1-\cos \frac{l_{(i)}}{\lambda })\frac{1}{N};U_{14}=-(l_{(i)}-\lambda \sin \frac{l_{(i)}}{\lambda })\frac{1}{N}\hfill \\ U_{23}=\frac{1}{\mathrm{EI}\lambda }\sin \frac{l_{(i)}}{\lambda };U_{24}=(1-\cos \frac{l_{(i)}}{\lambda })\frac{1}{N}\hfill \\ U_{33}=U_{44}=\cos \frac{l_{(i)}}{\lambda };U_{34}=\lambda \sin \frac{l_{(i)}}{\lambda };\hfill \\ U_{43}=-\frac{1}{\lambda }\sin \frac{l_{(i)}}{\lambda }\hfill \end{array}$$

Assuming that the whole beam structure is divided into n units, the transfer matrix of the left and right ends of the structure is written as follows

$$S_{l(n)}=\underset{i=1}{\overset{n}{\Pi }}{U}_i{S}_{0(1)}$$

(32)

Combining the boundary condition which is shown in Table 1, the basic state vector at the boundary node is obtained by solving the residual matrix. Then, the state vectors of all elements in the beam structure are obtained by the recursion method through the transfer matrix relations of each unit.

Program design of solving slippage length based on recursive algorithm

The parameters of transfer matrix are related to the state of slippage between layers, while the slippage between layers gradually unfolds with the increasing of outer load. Incremental iterative algorithm can be used for calculation. The loading process is divided into multiple load increment steps, the increment of state vector under each load increment step is calculated separately, and the load increment of the current step is added to the calculation result of the previous load step. The algorithm for solving the slippage length is shown in Figure 6.

Figure 6.

The main program flow chart and the transfer matrix subroutine flow chart.

Numerical simulation

Take a laminated beam with a span of 5 m as the analysis object. This laminated beam is composed of 11 layers of single-layer beams with a thickness of 7 mm, as shown in Figure 7. The material of the laminated beam of this numerical example is steel, Poisson’s ratio is 0.3, and the elastic modulus is 195 GPa. The beam is fixed at both ends. The initial tension in the beam is 106 N, and a gradually increasing concentration force is applied at the midspan position.

Figure 7.

The schematic drawing of the laminated beam.

The width of the laminated beam layers is assumed to be symmetrical and shown in Table 2.

Table 2.

Equivalent beam width of the PWS-91 parallel steel wire cable.

Layer number N	Width of the beam $b_i(\mathrm{mm})$
1, 11	25
2, 10	29
3, 9	33
4, 8	37
5, 7	41
6	45

The ANSYS FEM of the laminated beam is simulated by three-dimensional beam element BEAM188 with 6 degrees of freedom (DOF) in the space, and each laminated beam layer is divided into 500 elements in total. The rigid arm unit is established between the units of every two adjacent laminated beam layer elements, and the interlayer contact is simulated through the DOF-coupling of the common node of the rigid arms. Boundary conditions of the laminated beam model are fixed constraints.

Simulate the initial strain by applying axial tension of 106 N. The initial external load is 120 N, and the load increment for each load step is 5 N. The FEM is shown in Figure 8.

Figure 8.

Finite element model of the laminated beam.

The moment of inertia of the laminated beam is related to the slippage state, as shown in Figure 9.

Figure 9.

Relationship between the slippage state and the cross-section moment of inertia.

According to the presented equations (22)–(24) about slippage, it can be deduced that when the external load reaches 120 N, the laminated beam gets the critical slippage state. To verify the accuracy of the proposed algorithm, the results at different span locations were calculated by ANSYS and compared as shown in Figure 10.

Figure 10.

Comparison of the displacements and the shear forces calculated by ANSYS and proposed algorithm: (a) comparison of the displacement results calculated by FEM and the presented algorithm and (b) comparison of the shear force results calculated by FEM and the presented algorithm.

Figure 10 shows that the calculation results of the algorithm are very close to FEM simulation. Due to the rigid girder element in FEM, the shear locking phenomenon occurs and results in the simulation shear force which is smaller than the proposed algorithm, but the numerical errors are very small, as shown in Figure 10(b).

The shear force at the coupling node of the rigid girder element in FEM under the critical slippage state is extracted, that is, the shear force value of the contact layer in the laminated beams can be simulated, and the distribution rule is drawn as shown in Figure 11.

Figure 11.

The shear forces of the contact layers in FE model.

Figure 11 shows that the shear force of the laminated beam contact layer is quadratic in the transverse direction and antisymmetric quadratic distribution in the longitudinal direction with respect to the external load application position. The shear force distribution law solved by ANSYS is in accordance with the deduced results of the proposed algorithm.

According to previous equations, the ultimate interlayer friction resistance of this example is 12.3 N. Suppose that the shear force of the rigid girder element in the FEM is the integral of the unit shear stress on the contact layer area, the calculated shear force of the contact layer in Figure 11 shows that the maximum value of the shear force in rigid girder nodes is 11.9 N. Based on the law of interlayer slippage deduced in the previous sections, the first slippage occurs at both the fixed ends and the midspan at the same time with the superimposition of load substeps.

Under the current substep, if the shear force of the contact layer of a finite element unit is greater than the ultimate frictional resistance, then the longitudinal restraint of the coupling node of the rigid girder connected to this element is released and the shear force of the current contact layer is reversely applied to the node for simulating the friction between layers.

The number of slippage elements near the left fixed end in each substep is extracted and described as shown in Figure 12. The result describes the contact layer slippage development process in each of the four higher shear force contour areas. The first all-section slip state of the laminated beam happens in 10th substep and the length of the interlayer slippage at this state is 15 elements.

Figure 12.

Slippage development process of contact layers in the laminated beam.

Displacement comparison in all substeps between FEM and the proposed algorithm is shown in Figure 13. It indicates that the proposed algorithm is still accurate and efficient in analyzing the mechanical characteristics of the laminated beam when considering interlayer slippage. FEM needs to adjust the iterative process manually. The proposed algorithm based on incremental superposition method automatically completes iterative calculation and seems to be simple and efficient.

Figure 13.

Comparison of the displacement results in all substeps between FEM and presented program.

Figure 9 shows that the interlayer slippage reduces the stiffness of the laminated beam and increases the vertical displacement. The higher the proportion of slippage spread, the more significant the vertical displacement of the laminated beam is. Figure 14 shows the vertical displacement of midspan nodes with or without considering slippage in different substeps.

Figure 14.

Variation curve of displacements for the middle node.

It can be seen from Figure 14 that for a beam without considering the slippage, the displacement and force are in a linear relationship because the axial force is assumed to be constant during the deformation, whereas for the beam taking the slippage into consideration, the displacement tends to increase non-linearly.

Conclusion

This article provides a mechanical analysis of the interlayer slippage in the contact layer of the frictional laminated beam based on the section analysis method. A new algorithm for calculating the slippage effect of laminated beam is proposed based on the incremental iterative recursive method and TMM. Comparisons against FEM numerical analysis results indicate that the present algorithm is accurate:

1. The interlayer slippage develops from the neutral axis to the outer layer symmetrically in the cross-section direction, and the interlayer slippage area in the longitudinal direction is related to the external loads applying position and the boundaries.

2. Interlayer slippage in the laminated beam will destroy the integrity of the section, thereby reducing the bending stiffness of the section and increasing the maximum bending stress. Therefore, even in the case of components mainly subjected to axial forces, the possibility of the components being subjected to lateral loads should be considered in structural strength design or fatigue calculation. If the components are subjected to lateral forces, the slippage effect needs to be considered according to the actual slippage effect of bending stress.

Future work

The continuation of this research will include the following:

1. Mechanical analysis of local bending behavior of the laminated beam by interlayer slip effect in different axial force conditions.

2. Better understanding of the effect of interlayer slip on the vibration characteristics of the frictional laminated beam.

3. Usage of the derivatives of eigen parameters method to identify the interlayer slip damage of the laminated beam.