Numerical Analysis of Cracked Double-Beam Systems

1. Introduction

The literature regarding the mechanical behaviour of beams is very rich. In the majority of the earliest papers, the solutions presented in the literature can be divided into theoretical and numerical ones. In this context, there exist a lot of studies on evaluating the bending, buckling, post-buckling and vibration behaviours of beams using Euler- Bernoulli and Timoshenko beams models [1–10].

It is well-known that the structural behaviour of beams is sensitive to presence of cracks. In a beam, their occurrence, introduces local flexibility, a localised increase in the bending flexibility, which may lead to excessive deflections and unexpected failures. Because of their practical relevance, many studies are performed to explore the static deflection and vibration response of the beams with different boundary conditions and resting on various elastic foundations [11–17]. Some works are cited herein. In [11] Biondi and Caddemi have studied the problem of the integration of the static governing equations for uniform Euler−Bernoulli beam with two types discontinuities and have presented a closed form solutions of governing differential equations. Cicirello and Palmeri [12] dealt with the static analysis of pre-damaged Euler−Bernoulli beams with any number of unilateral cracks and subjected to tensile or compression forces combined with arbitrary transverse loads. Khaji et al. in [13] have developed an analytical approach for crack identification procedure in uniform Timoshenko beams with an open edge crack, based on bending vibration measurements. Ghannadiasl and Khodapanah in [14] have presented an analytical solution of dynamic analysis of cracked Euler-Bernoulli beam on the elastic foundation subjected to the concentrated load. Also, the effects of depth and location of the crack on natural frequency and deflection of the cracked beam on an elastic foundation have been evaluated. In his PhD thesis [15], Batihan dealt with the transverse vibration of a cracked beam on an elastic foundation and the effect of crack and foundation parameters on transverse vibration natural frequencies has been presented. Yang et al. have studied the bending deformation of the Timoshenko beam with switching cracks and the influence of the beam slenderness ratio, the crack depth, and the external load on the crack state and bending performances of the cracked beam have been assessed [16]. Finally, Alijani et al. have investigated the static behavior of cracked Euler−Bernoulli beams resting on an elastic foundation implementing analytical, approximate and numerical approaches [17].

Although the single beam models under various kind of loadings and boundary conditions are the most studied structural solutions, they cannot be employed in many engineering applications, such as sandwich or composites beams, nanostructures, adhesively bonded joints, floating-slab tracks, continuous dynamic vibration absorbers. Starting from these assumptions, double-beam systems have attracted much attention from researchers and engineers and they play an important role in many fields of structural and foundation engineering.

Double Beam System (BS) models are structural models consisting of two parallel beams structure interconnected by a uniformly distributed elastic layer which is generally considered as Winkler elastic medium. Beam theories are assumed to govern the beams and elastic foundation models are assumed to represent the elastic layers. Thanks to their remarkable structural properties, such as better vibration absorption than a single beam, lighter weight, and higher strength and stiffness, led to an explosion of interest within the scientific community. As a result, progressive research activities regarding BS have been ongoing in recent years, and they have been used widely in many engineering fields. For instance, special configurations of multiple pipeline systems, auxiliary underground structures (i.e. passageway and drainage systems) can be represented by double Beam System (BS) models. In civil engineering field, this system helps to reduce the energy from the earthquake and based on this peculiar characteristic, many vibration absorbers have been developed. Although static and vibration analysis of beams resting on elastic foundation is a widely studied topic, few works for static analysis of elastically connected system can be found in the literature. In this sense, applying theoretical and numerical methods, some significant results have been obtained in the study of beams and nanotubes [18–23].

The present paper deals with the static analysis of a uniform cracked double beam system in presence of a Winkler medium. The structure is also supposed to be constrained at both ends by elastically flexible springs, with transverse stiffness and rotational stiffness. Energy principle is formulated for the static analysis of the double beam system and the goverming equations is solved analytically and numerically. Cell-Discretization Method (CDM) is employed to discretize and solve the governing equations and boundary conditions. This numerical method has already been used by the authors [23–26], by Raithel and Franciosi [27] and Franciosi and Franciosi [28] for different structural problems. The analysis is performed reducing the two beams to a set of rigid bars linked together by means of elastic sections (elastic cells), in which the stiffness of the beam is properly concentrated. In this way, the structure is reduced to a system with finite number of degrees of freedom, and the global stiffness matrix can be easily calculated.

Numerical examples are then given to demonstrate the reliability and effectiveness of the current model. Also, in order to outline the efficiency of the proposed computational model, the Authors have considered the derivation of a beam finite element and have presented the derivation of stiffness matrix for cracked double-beam system. The results obtained by closed-form solution are compared by Cell-discretization method (CDM) and Finite element method (FEM). Accuracy of the results has been evaluated by making comparisons with the results in literature and reported in bibliography. Also, the effects of different parameters, such as crack depth and position, boundary conditions, elastic medium stiffness, slenderness on the static behaviour of the structure are investigated. It is demonstrated that the proposed algorithm provides a simple and powerful tool in dealing with the static analysis of a double-beam system. Finally, some concluding remarks are made.

2. Problem formulation

Let us consider the system which is composed of two parallel beams, having the same length L and translational and rotational elastic constraints at their ends, as shown in Figure 1. Assuming that the material and geometrical properties of the two beams are the same, let be E and I the Young’s modulus and moment of inertia, respectively. It is assumed that the two beams are joined by a Winkler-type medium with modulus $k_m$. The upper beam is subjected to a uniformly distributed load, q.

Based on the Euler-Bernoulli theory and applying the variational formulation, the total potential energy of the system under consideration assumes the following form:

$$\begin{array}{cc}\hfill E_1+\mathrm{P}=& \frac{1}{2}\sum _{i=1}^2{\int }_0^{\mathrm{L}}\left[{\mathrm{EI}}_j{\left(\frac{{\partial }^2{\mathrm{v}}_j\left(\mathrm{z}\right)}{\partial {\mathrm{z}}^2}\right)}^2\right]\mathrm{dz}+\frac{1}{2}{k}_{\mathrm{jRL}}{\mathrm{v}}_j^{{ }^{\prime }2}(\mathrm{z}=0)+\frac{1}{2}{k}_{\mathrm{jTL}}{\mathrm{v}}_j^2(\mathrm{z}=0)+\hfill \\ \hfill & \frac{1}{2}{k}_{\mathrm{jRR}}{\mathrm{v}}_j^{{ }^{\prime }2}(\mathrm{z}=\mathrm{L})+\frac{1}{2}{k}_{\mathrm{jTR}}{\mathrm{v}}_j^2(\mathrm{z}=\mathrm{L})-{\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{v}}_1\left(z\right)\mathrm{dz}+\frac{1}{2}{k}_{\mathrm{jds}}{\left(\Delta v_j^{{ }^{\prime }}\left(L_{\mathrm{jc}}\right)\right)}^2\hfill \end{array}$$

(1)

$$E_2=\frac{1}{2}{\int }_0^{\mathrm{L}}\left[k_m{\left({\mathrm{v}}_2\left(\mathrm{z}\right)-{\mathrm{v}}_1\left(\mathrm{z}\right)\right)}^2\right]\mathrm{dz}$$

(2)

where $E_1$ is strain energy of two beams, P is potential energy of applied load and $E_2$ the elastic energy which includes the contribution of the elastic medium. Let v(z) the transverse displacement, z is the spatial coordinate, $k_{\mathrm{jRL}}$ and $k_{\mathrm{jTL}}$ are rotational and translational stiffness at z=0 and $k_{\mathrm{jRR}}$ e $k_{\mathrm{jTR}}$ rotational and translational stiffnesses at z=L, respectively, with j=1,2. The index j=1,2 refers to the order of the beams: the upper beam is indicated by j=1 and the lower beam by j=2, respectively. Finally, $k_{\mathrm{jds}}$ is equivalent rotational springs stiffness of the two corresponding sections of two beams and $L_{\mathrm{jc}}$ is their abscissa.

Applying the principle of stationary potential energy, the following dimensionless form of the static governing equations of the Euler−Bernoulli beam model is obtained:

$$\begin{array}{c}\hfill {\mathrm{v}}_1{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_1\left(\zeta \right)-{\mathrm{v}}_2\left(\zeta \right)\right)=\frac{\mathrm{q}{\mathrm{L}}^4}{\mathrm{EI}}\\ \hfill {\mathrm{v}}_2{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_2\left(\zeta \right)-{\mathrm{v}}_1\left(\zeta \right)\right)=0\end{array}$$

(3)

in which $\zeta =\frac{z}{\mathrm{L}}$. In this way, we shift the domain [0,L] to the domain [0,1] .

As it is known, the presence of a crack or concentrated force introduces a discontinuity and therefore the displacement function of two beams is not a smooth function of z. Consequently, in order to calculate the displacement function of the two beams, a separate differential equations system must be written for each beam segment between the discontinuities. Let n be the number of discontinuities of the upper section and s those of the lower section, $n+s+1=m$ differential equations systems are written. Therefore the number of systems of equations is the sum of the discontinuities between the upper and lower sections of the two beams.

For the double beam system under consideration, one gets:

$$\begin{array}{c}\hfill {\mathrm{v}}_1{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_1\left(\zeta \right)-{\mathrm{v}}_{m+1}\left(\zeta \right)\right)=\frac{\mathrm{q}{\mathrm{L}}^4}{\mathrm{EI}}0<\zeta <{\gamma }_1\\ \hfill {\mathrm{v}}_{m+1}{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_{m+1}\left(\zeta \right)-{\mathrm{v}}_1\left(\zeta \right)\right)=00<\zeta <{\gamma }_1\end{array}$$

(4)

$$\begin{array}{c}\hfill {\mathrm{v}}_i{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_i\left(\zeta \right)-{\mathrm{v}}_{m+i}\left(\zeta \right)\right)=\frac{\mathrm{q}{\mathrm{L}}^4}{\mathrm{EI}}{\gamma }_{i-1}<\zeta <{\gamma }_{i}i=2,...,m-1\\ \hfill {\mathrm{v}}_{m+i}{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_{m+i}\left(\zeta \right)-{\mathrm{v}}_i\left(\zeta \right)\right)=0{\gamma }_{i-1}<\zeta <{\gamma }_{i}i=2,...,m-1\end{array}$$

(5)

$$\begin{array}{c}\hfill {\mathrm{v}}_m{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_m\left(\zeta \right)-{\mathrm{v}}_{2m}\left(\zeta \right)\right)=\frac{\mathrm{q}{\mathrm{L}}^4}{\mathrm{EI}}{\gamma }_{m-1}<\zeta <1\\ \hfill {\mathrm{v}}_{2m}{ }^{\prime \prime \prime \prime }\left(\zeta \right)+\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}\left({\mathrm{v}}_{2m}\left(\zeta \right)-{\mathrm{v}}_m\left(\zeta \right)\right)=0{\gamma }_{m-1}<\zeta <1\end{array}$$

(6)

From the first equations of (4)–(6) the displacements ${\mathrm{v}}_{m+i}$($\zeta$) and ${\mathrm{v}}_{2m}\left(\zeta \right)$ are obtained, respectively, which substituted in the second equation of (4)–(6), lead to:

$${\mathrm{v}}_{m+1}\left(\zeta \right)=\frac{{\mathrm{v}}_1{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_1\left(\zeta \right)-\frac{p}{{\alpha }^4}$$

(7)

$${\mathrm{v}}_{m+i}\left(\zeta \right)=\frac{{\mathrm{v}}_i{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_i\left(\zeta \right)-\frac{p}{{\alpha }^4}i=2,...,m-1$$

(8)

$${\mathrm{v}}_{2m}\left(\zeta \right)=\frac{{\mathrm{v}}_m{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_m\left(\zeta \right)-\frac{p}{{\alpha }^4}$$

(9)

and the following system of equations yields:

$${\mathrm{v}}_1{ }^{\prime \prime \prime \prime \prime \prime \prime \prime }\left(\zeta \right)+2{\alpha }^4{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime }}\left(z\right)=p{\alpha }^{4}0<\zeta <{\gamma }_1$$

(10)

$${\mathrm{v}}_i{ }^{\prime \prime \prime \prime \prime \prime \prime \prime }\left(\zeta \right)+2{\alpha }^4{\mathrm{v}}_i^{{ }^{\prime \prime \prime \prime }}\left(z\right)=p{\alpha }^4{\gamma }_{i-1}<\zeta <{\gamma }_{i}i=2,...,m-1$$

(11)

$${\mathrm{v}}_m{ }^{\prime \prime \prime \prime \prime \prime \prime \prime }\left(\zeta \right)+2{\alpha }^4{\mathrm{v}}_m^{{ }^{\prime \prime \prime \prime }}\left(z\right)=p{\alpha }^4{\gamma }_{m-1}<\zeta <1$$

(12)

setting

$$\alpha =\sqrt[4]{\frac{k_m{\mathrm{L}}^4}{\mathrm{EI}}};p=\frac{\mathrm{q}{\mathrm{L}}^4}{\mathrm{EI}}$$

(13)

The general solutions are:

$$\begin{array}{cc}\hfill {\mathrm{v}}_i\left(\zeta \right)=& a_{1+(i-1)8}+a_{2+(i-1)8}\zeta +a_{3+(i-1)8}{\zeta }^2+a_{4+(i-1)8}{\zeta }^3+\hfill \\ & a_{5+(i-1)8}\mathrm{Cos}\mathrm{h}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]\mathrm{Sin}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]+a_{6+(i-1)8}\mathrm{Cos}\mathrm{h}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]\mathrm{Cos}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]+\hfill \\ & a_{7+(i-1)8}\mathrm{Sin}\mathrm{h}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]\mathrm{Sin}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]+a_{8+(i-1)8}\mathrm{Sin}\mathrm{h}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]\mathrm{Cos}\left[\frac{\alpha }{\sqrt[4]{2}}\zeta \right]+\hfill \\ & \frac{p{\zeta }^4}{48}i=1,m\hfill \end{array}$$

(14)

where $a_j$ are eight constants that are determined, for each segment, by imposing the boundary conditions. At the n discontinuity point for the upper beam with traslational and rotational constraints, the boundary conditions are defined as:

$$\begin{array}{c}\hfill {\mathrm{v}}_1^{{ }^{\prime \prime }}(\zeta =0)-K_{1\mathrm{RL}}{\mathrm{v}}_1^{{ }^{\prime }}(\zeta =0)=0\\ \hfill {\mathrm{v}}_1^{{ }^{\prime \prime \prime }}(\zeta =0)+K_{1\mathrm{TL}}{\mathrm{v}}_1(\zeta =0)=0\end{array}$$

(15)

$$\begin{array}{c}\hfill {\mathrm{v}}_i\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{i+1}\left(\zeta ={\gamma }_i\right)=0i=2,m-1\\ \hfill {\mathrm{v}}_1^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{i+1}^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)+\frac{\Psi }{\mathrm{L}}{\mathrm{v}}_1^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)=0\\ \hfill {\mathrm{v}}_i^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{i+1}^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)=0\\ \hfill {\mathrm{v}}_i^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{i+1}^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)+F_t=0\end{array}$$

(16)

$$\begin{array}{c}\hfill {\mathrm{v}}_m^{{ }^{\prime \prime }}(\zeta =1)+K_{1\mathrm{RR}}{\mathrm{v}}_m^{{ }^{\prime }}(\zeta =1)=0\\ \hfill {\mathrm{v}}_m^{{ }^{\prime \prime \prime }}(\zeta =1)-K_{1\mathrm{TR}}{\mathrm{v}}_m(\zeta =1)=0\end{array}$$

(17)

For the lower beam, having translational and rotational elastic constraints at the ends and subjected to a concentrated force, the boundary conditions are:

$$\begin{array}{c}\hfill {\mathrm{v}}_{m+1}^{{ }^{\prime \prime }}(\zeta =0)-K_{2\mathrm{RL}}{\mathrm{v}}_{m+1}^{{ }^{\prime }}(\zeta =0)=0\\ \hfill {\mathrm{v}}_{m+1}^{{ }^{\prime \prime \prime }}(\zeta =0)+K_{2\mathrm{TL}}{\mathrm{v}}_{m+1}(\zeta =0)=0\end{array}$$

(18)

$$\begin{array}{c}\hfill {\mathrm{v}}_{m-1+i}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{m+i}\left(\zeta ={\gamma }_i\right)=0i=2,m\\ \hfill {\mathrm{v}}_{m-1+i}^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{m+i}^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)+\frac{\Psi }{\mathrm{L}}{\mathrm{v}}_{m-1+i}^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)=0\\ \hfill {\mathrm{v}}_{m-1+i}^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{m+i}^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)=0\\ \hfill {\mathrm{v}}_{m-1+i}^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)-{\mathrm{v}}_{m+i}^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)+F_t=0\end{array}$$

(19)

$$\begin{array}{c}\hfill {\mathrm{v}}_{2m}^{{ }^{\prime \prime }}(\zeta =1)+K_{2\mathrm{RR}}{\mathrm{v}}_{2m}^{{ }^{\prime }}(\zeta =1)=0\\ \hfill {\mathrm{v}}_{2m}^{{ }^{\prime \prime \prime }}(\zeta =1)-K_{2\mathrm{TR}}{\mathrm{v}}_{2m}(\zeta =1)=0\end{array}$$

(20)

where

$$\begin{array}{c}\hfill F_t=\frac{F{\mathrm{L}}^3}{\mathrm{EI}}{K}_{j\mathrm{RL}}=\frac{k_{\mathrm{jRL}}\mathrm{L}}{\mathrm{EI}};K_{j\mathrm{RR}}=\frac{k_{\mathrm{jRR}}\mathrm{L}}{\mathrm{EI}};\\ \hfill K_{\mathrm{jTL}}=\frac{k_{\mathrm{jTL}}{\mathrm{L}}^3}{\mathrm{EI}};K_{\mathrm{jTR}}=\frac{k_{\mathrm{jTR}}{\mathrm{L}}^3}{\mathrm{EI}};j=1,2\end{array}$$

(21)

denote dimensionless parameters of rotational and translational stiffnesses at two ends, for $\zeta$=0 and $\zeta$=1, respectively and taking into account the presence of concentrated non-dimensional force.

By substituting the following equations

$${\mathrm{v}}_{m+1}(\zeta =0)=\frac{{\mathrm{v}}_1{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_1\left(\zeta \right)-\frac{p}{{\alpha }^4}$$

(22)

$${\mathrm{v}}_{m+i}\left(\zeta ={\gamma }_i\right)=\frac{{\mathrm{v}}_i{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_i\left(\zeta \right)-\frac{p}{{\alpha }^4}$$

(23)

$${\mathrm{v}}_{2m}\left(\zeta \right)=\frac{{\mathrm{v}}_n{ }^{\prime \prime \prime \prime }\left(\zeta \right)}{{\alpha }^4}+{\mathrm{v}}_m\left(\zeta \right)-\frac{p}{{\alpha }^4}$$

(24)

into the equations (18-20), one obtains the following system of equations:

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime \prime \prime }}(\zeta =0)}{{\alpha }^4}+{\mathrm{v}}_1^{{ }^{\prime \prime }}(\zeta =0)-K_{2\mathrm{RL}}\left(\frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime \prime }}(\zeta =0)}{{\alpha }^4}+{\mathrm{v}}_1^{{ }^{\prime }}(\zeta =0)\right)=0\\ \hfill \frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime \prime \prime \prime }}(\zeta =0)}{{\alpha }^4}+{\mathrm{v}}_1^{{ }^{\prime \prime \prime }}(\zeta =0)+K_{2\mathrm{TL}}\left(\frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime }}(\zeta =0)}{{\alpha }^4}+{\mathrm{v}}_1(\zeta =0)-\frac{p}{{\alpha }^4}\right)=0\end{array}$$

(25)

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_i^{{ }^{\prime \prime \prime \prime }}\left(\left(\zeta ={\gamma }_i\right)\right.}{{\alpha }^4}+{\mathrm{v}}_i\left(\zeta ={\gamma }_i\right)-\frac{{\mathrm{v}}_{i+1}^{{ }^{\prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}-{\mathrm{v}}_{i+1}\left(\zeta ={\gamma }_i\right)=0i=2,m-1\end{array}$$

(26)

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_i^{{ }^{\prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}+{\mathrm{v}}_i^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)-\frac{{\mathrm{v}}_2^{{ }^{\prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}+\\ & -{\mathrm{v}}_2^{{ }^{\prime }}\left(\zeta ={\gamma }_i\right)+\frac{\Psi }{\mathrm{L}}\left(\frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}+{\mathrm{v}}_1^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)\right)=0\hfill \end{array}$$

(27)

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_i^{{ }^{\prime \prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}+{\mathrm{v}}_i^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)-\frac{{\mathrm{v}}_2^{{ }^{\prime \prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}-{\mathrm{v}}_2^{{ }^{\prime \prime }}\left(\zeta ={\gamma }_i\right)=0\end{array}$$

(28)

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_1^{{ }^{\prime \prime \prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}+{\mathrm{v}}_i^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)-\frac{{\mathrm{v}}_2^{{ }^{\prime \prime \prime \prime \prime \prime \prime }}\left(\zeta ={\gamma }_i\right)}{{\alpha }^4}-{\mathrm{v}}_2^{{ }^{\prime \prime \prime }}\left(\zeta ={\gamma }_i\right)+F_t=0\end{array}$$

(29)

$$\begin{array}{c}\hfill \frac{{\mathrm{v}}_m^{{ }^{\prime \prime \prime \prime \prime \prime }}(\zeta =1)}{{\alpha }^4}+{\mathrm{v}}_m^{{ }^{\prime \prime }}(\zeta =1)+K_{2\mathrm{RR}}\left(\frac{{\mathrm{v}}_m^{{ }^{\prime \prime \prime \prime \prime }}(\zeta =1)}{{\alpha }^4}+{\mathrm{v}}_m^{{ }^{\prime }}(\zeta =1)\right)=0\\ \hfill \frac{{\mathrm{v}}_m^{{ }^{\prime \prime \prime \prime \prime \prime \prime }}(\zeta =1)}{{\alpha }^4}+{\mathrm{v}}_m^{{ }^{\prime \prime \prime }}(\zeta =1)-K_{2\mathrm{TR}}\left(\frac{{\mathrm{v}}_m^{{ }^{\prime \prime \prime \prime }}(\zeta =1)}{{\alpha }^4}+{\mathrm{v}}_m(\zeta =1)-\frac{p}{{\alpha }^4}\right)=0\end{array}$$

(30)

3. Cracks: modelling and method of solution

3.1. Modelling of cracks: an overview of the discrete spring model

The identification of the location and depth of a crack in beam type structures is an important topic of structural health monitoring and has been the subject of a significant body of research. The literature concerning the modelling of cracks is very rich and various models, such as local stiffness reduction, discrete spring models, and complex models in two or three dimensions, have been proposed in the technical literature. In this topic, the state of the art can be found in a review works by Alijani et al. [17] and by Palmeri and Cicirello [29]. In both papers the authors report a coherent yet concise review of as many of these publications as possible and the main themes treated are the modelling and simulation of static behaviour of cracked beams.

According to the classification by Friswell and Penny [30], the proposed approach falls in the broad category of “discrete spring models”, being equivalent to an internal hinge coupled with a linear elastic spring. This model is widely adopted when the focus of the structural analysis is on the global performance of frames and trusses, rather than on phenomena of crack initiation and crack propagation. Although very simple, the discrete spring model proves to be very efficient for static problems. The idea of treating cracked beams with equivalent linear springs at the cracks position is based on the partition of each member into undamaged pieces between two consecutive cracks. One of the main advantages of this model is the effective representation of the crack in terms of position and severity. Different attempts can be found in the literature to provide the values for the rotational spring stiffness using crack parameters as depth or geometry. Among them, Palmeri and Cicirello [29] have investigated the static analysis of multi-cracked Euler–Bernoulli and Timoshenko beams using this relationship. Afterward, Okamura et al. [31] have showed the buckling behavior of cracked columns with a rectangular cross-section. Furthermore, Ricci and Viola [32] have extended the method of Kienzler and Herrmann [33] to evaluate stress intensity factors of cracked beams and bars and have derived two relationships between stress intensity factor and rotational spring stiffness. In the following analysis, two relationships between stress intensity factor and rotational spring stiffness have been considered, so as introduced by Alijani et al. in [17].

Figure 2 shows a beam with a crack in the non-dimensional position ${\gamma }_i$ and depth a. Also, it is assumed that the shear effect as well as axial and torsional loads are neglected. According to this assumption, the stress intensity factors may be evaluated through the following equations:

$$K_I=\frac{6M}{b{h}^2}\sqrt{\pi a}{F}_m\left(\xi \right)0\le \xi \le 0.6,\xi =\frac{a}{h}$$

(31)

$$K_I=\frac{3.99M}{bh\sqrt{h}\sqrt{{\left(1-{\xi }^2\right)}^3}}0.6\le \xi \le 1.0,\xi =\frac{a}{h}$$

(32)

where $F_m\left(\xi \right)$ is given from the following geometric function

$$F_M\left(\xi \right)=\sqrt{\left(\frac{2}{\pi \xi }\right)tan\frac{\pi \xi }{2}}\frac{0.923+0.199{\left(1-sin\left(\frac{\pi \xi }{2}\right)\right)}^4}{cos\left(\frac{\pi \xi }{2}\right)}$$

(33)

Finally, the following relationship is defined:

$$\frac{1}{k_{\mathrm{ds}}}=\frac{2b\left(1-{\vee }^2\right)}{E}{\int }_0^a{\left(\frac{K_I}{M}\right)}^{2}da$$

(34)

which is used to determine a rotational spring stiffness factor equivalently in terms of the geometric and material parameters of the crack and where $\Psi$ of Equation (27) is equal to $\frac{\mathrm{EI}}{k_{\mathrm{ds}}}$ and $k_{\mathrm{ds}}$ is the stiffness of the i-th equivalent spring.

In the present paper the rotational spring stiffness factor is assumed to be different for each beam. In particular ${\Psi }_1$ and ${\Psi }_2$ indicate the rotational spring stiffness factor for upper and lower beams, respectively.

3.2. Method of Solution: Cell-discretization method (CDM)

The Cell-discretization method (CDM) is an efficient numerical method for the solution of linear partial differential equations. It has been becoming an important tool in field of the structural engineering, thanks to its approximation abilities and easiness to be implemented. The advent of sophisticated and totally generalized discretization tools, such as Finite Element method (FEM), Boundary Element Method (BEM), allowed to simulate the behaviour of structures taking into account several variables due to the removal of as many as simplified hypotheses, but on the other hand such procedures may induce to lose the physical sense of the real behaviour of the structures that should be always at the basis of engineering studies. In this sense, the CDM may be regarded as a technique able to address such issues. Since the beginning of the 20th century, this method found several applications, for instance: the dynamics and stability of arches; masonry arches; static and dynamical analysis of Euler-Bernoulli beams under several load and boundary conditions; static and dynamical analysis of Timoshenko and Rayleigh beams; static analysis of plates under several load and boundary conditions [23–28]. More recently, some of the present authors applied the method to the dynamical analysis of single- and double-walled carbon nanotubes, by taking into account non - local effects [23,26] and they obtained results showing that the method is able to describe the nanostructure behaviour satisfactorily with a little computational effort. In several papers, the procedure is shown to be very versatile and able to work in any case on a finite number of Lagrangian parameters by bringing the solution into the alveo of the usual numerical analysis methods.

In the present paper, the two beams are reduced to a set of t rigid bars with the same length $\mathit{l}$, linked together by n=t+1 elastic cells (see Figure 3). The moment of inertia $I_j$, with j=1,2, will be evaluated at the cells abscissa, obtaining the concentrated stiffness $k_{1i}=\frac{{\mathrm{EI}}_{1i}}{\mathit{l}}$ and $k_{2i}=\frac{{\mathrm{EI}}_{2i}}{\mathit{l}}$, for the upper and lower beams, respectively. Both these quantities can be organized into the so-called unassembled stiffness diagonal matrix ${\mathbf{k}}_j$ with dimension (n x n), j=1,2, for each two beams.

If the cross-section is non-uniform, for each rigid section of length $\mathit{l}$ the average of the inertia across the elastic cell can be considered.

In this way, the structure is reduced to a classical holonomic system, with 2n degrees of freedom, in particular, n vertical displacements ${\mathrm{v}}_{1i}$ for the upper beam, and n vertical displacements ${\mathrm{v}}_{2i}$, for the lower beam, at the cells abscissa will be conveniently assumed as Lagrangian coordinates and will be organized into the 2n-dimensional vector v. Moreover, for the upper and lower beams the n-1 rotations of the rigid bars can be calculated as a function of the Lagrangian coordinates as follows:

$$\begin{array}{c}\hfill {\varphi }_{1,i}=\frac{{\mathrm{v}}_{1,i+1}-{\mathrm{v}}_{1,i}}{l}\\ \hfill {\varphi }_{2,i}=\frac{{\mathrm{v}}_{2,i+1}-{\mathrm{v}}_{2,i}}{l}\end{array}$$

(35)

or, in matrix form ${\varphi }_{\mathbf{1}}=$V${\mathrm{v}}_1$ and ${\varphi }_{\mathbf{2}}=$V${\mathrm{v}}_2$ where $\mathrm{V}$ is a rectangular transfer matrix with n-1 rows and n columns.

The relative rotations between the two faces of the elastic cells are given by:

$${\psi }_{j,1}={\varphi }_{j,1},{\psi }_{j,i}={\varphi }_{j,i}-{\varphi }_{j,i-1},{\psi }_{j,n}=-{\varphi }_{j,n-1}$$

(36)

or, in matrix form ${\psi }_{\mathbf{1}}$ =$\Delta$${\varphi }_{\mathbf{1}}$ for upper rigid bar and ${\psi }_{\mathbf{2}}$ =$\Delta$${\varphi }_{\mathbf{2}}$, for lower rigid bar, where $\Delta$ is another rectangular transfer matrix with n rows and n-1 columns.

The strain energies $L_{\mathrm{je}}$, with j=1,2, (the first two terms of Equation (1)) are given by:

$$\begin{array}{c}\hfill L_{1e}=\frac{1}{2}\sum _{i=1}^n{k}_{1,\mathrm{ii}}{\psi }_{1,i}^2\\ \hfill L_{2e}=\frac{1}{2}\sum _{i=1}^n{k}_{2,\mathrm{ii}}{\psi }_{2,i}^2\end{array}$$

(37)

and they are concentrated at the cells of the upper and lower beams, respectively.

The strain energies should be expressed as functions of the Lagrangian coordinates as follows:

$$\begin{array}{c}\hfill L_{1e}=\frac{1}{2}{\psi }_1^T{\mathbf{k}}_1{\psi }_1=\frac{1}{2}{\varphi }_1^T{\Delta }^T{\mathbf{k}}_1\Delta {\varphi }_1=\frac{1}{2}{{\mathrm{v}}_{\mathbf{1}}}^T\left({\mathbf{V}}^T{\Delta }^T{\mathbf{k}}_1\Delta \mathbf{V}\right){\mathrm{v}}_{\mathbf{1}}\\ \hfill L_{2e}=\frac{1}{2}{\psi }_2^T{\mathbf{k}}_2{\psi }_2=\frac{1}{2}{\varphi }_2^T{\Delta }^T{\mathbf{k}}_2\Delta {\varphi }_2=\frac{1}{2}{{\mathrm{v}}_{\mathbf{2}}}^T\left({\mathbf{V}}^T{\Delta }^T{\mathbf{k}}_2\Delta \mathbf{V}\right){\mathrm{v}}_{\mathbf{2}}\end{array}$$

(38)

so that, the total strain energy can be expressed as:

$$L_e=\frac{1}{2}{\mathrm{v}}^T\left(\begin{array}{cc}{\mathbf{K}}_{\mathbf{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_{\mathbf{2}}\end{array}\right)\mathrm{v}$$

(39)

where ${\mathbf{K}}_{\mathbf{1}}$=$\left({\mathbf{V}}^T{\Delta }^T{\mathbf{k}}_1\Delta \mathbf{V}\right)$ and${\mathbf{K}}_{\mathbf{2}}$=$\left({\mathbf{V}}^T{\Delta }^T{\mathbf{k}}_2\Delta \mathbf{V}\right)$. The global assembled stiffness matrix K,with 2n rows and 2n columns, assumes the following form:

$$\mathbf{K}=\left(\begin{array}{cc}{\mathbf{K}}_{\mathbf{1}}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_{\mathbf{2}}\end{array}\right)$$

(40)

The potential energy, as function of the Lagrangian coordinates, is given by:

$$P_1=\sum _{i=1}^n{q}_i{\mathrm{v}}_{1,i}$$

(41)

Being the double beam system subjected to a uniformly distributed load q

$$\mathrm{q}=\left(q_1,...,q_i,...,q_n\right)$$

(42)

with

$$\begin{array}{c}\hfill q_1=ql/2;q_n=ql/2;\\ \hfill q_i=qli=2,...,n-1\end{array}$$

(43)

the global assembled load matrix Q with 2n rows and 2n columns, assumes the following form:

$$\mathbf{Q}=\left(\begin{array}{cc}\mathbf{q}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right)$$

(44)

The strain energy due to the elastic medium, Equation (2), can be expressed as:

$$E_2=\frac{1}{2}{\mathrm{v}}^T\mathrm{C}{\mathrm{v}}_1+\frac{1}{2}{\mathrm{v}}_2^T\mathrm{C}{\mathrm{v}}_2-{\mathrm{v}}_1^T\mathrm{C}{\mathrm{v}}_2$$

(45)

The terms of matrix $\mathrm{C}$ are given by:

$$\begin{array}{c}\hfill C_{i,i}=2\frac{l}{3}{k}_m,i=3,n-3\\ \hfill C_{i+1,i}=C_{ii,+1}=\frac{l}{3}{k}_m,i=2,n-2\\ \hfill C_{1,1}=\frac{l}{6}{k}_m,C_{2,2}=\frac{l}{3}{k}_m,C_{1,2}=C_{2,1}=\frac{l}{12}{k}_m\end{array}$$

(46)

The matrix $\mathrm{C}$ with n rows and n columns and have half-band widths equal to 2 takes the following form:

$${\mathrm{C}}_t=\left(\begin{array}{cc}\mathrm{C}& -\mathrm{C}\\ -\mathrm{C}& \mathrm{C}\end{array}\right)$$

(47)

Finally, the static governing equation can be written as:

$${\mathbf{K}}_t\mathrm{v}=\mathrm{Q}$$

(48)

with ${\mathbf{K}}_t$ is the global assembled stiffness matrix:

$${\mathbf{K}}_t=\mathrm{K}+{\mathrm{C}}_t$$

(49)

3.2.1. Boundary conditions in presence of a crack

Finally, starting to the strain energy terms of the flexible constraints at the ends in the eq (1), the assembled stiffness matrix K must be modified as follows:

$$\begin{array}{c}\hfill K[1,1]=K[1,1]+k_{1\mathrm{TL}};\\ \hfill K[n,n]=K[n,n]+k_{1\mathrm{TR}};\\ \hfill K[n+1,n+1]=K[n+1,n+1]+k_{2\mathrm{TL}};\\ \hfill K[2n,2n]=K[2n,2n]+k_{2\mathrm{TR}};\end{array}$$

(50)

The rotational stiffnesses of the constraints of each beam can be taken into account by summing up the corresponding flexibilities of the rigid bars and one gets:

$$\begin{array}{c}\hfill k_1[1,1]=\frac{k_1[1,1]k_{1\mathrm{RL}}}{k_{1\mathrm{RL}}+k_1[1,1]}\\ \hfill k_1[n,n]=\frac{k[n,n]k_{1\mathrm{RR}}}{k_{1\mathrm{RR}}+k_1[n,n]}\\ \hfill k_2[1,1]=\frac{k_2[1,1]k_{2\mathrm{RL}}}{k_{2\mathrm{RL}}+k_2[1,1]}\\ \hfill k_2[n,n]=\frac{k_2[n,n]k_{2\mathrm{RR}}}{k_{2\mathrm{RR}}+k_2[n,n]}\end{array}$$

(51)

These terms will be organized in two matrices ${\mathbf{k}}_1$ and ${\mathbf{k}}_2$ furnished in (40). In presence of a crack at the $j^{\mathrm{ma}}$ cell abscissa, the local stiffness of the upper beam is given by:

$$k_1[j+1,j+1]=\frac{k_1[j+1,j+1]k_{1\mathrm{ds}}}{k_{1\mathrm{ds}}+k_1[j+1,j+1]}$$

(52)

whereas at the $s^{\mathrm{ma}}$ cell abscissa relative to the lower beam, the local stiffness is:

$$k_2[s+1,s+1]=\frac{k_2[s+1,s+1]k_{2\mathrm{ds}}}{k_{2\mathrm{ds}}+k_2[s+1,s+1]}$$

(53)

if the height of crack is different in the two beams. The values $k_{1\mathrm{ds}}$ and $k_{2\mathrm{ds}}$ are derived from Equation (34).

3.3. Method of solution: Finite Element Method (FEM)

In order to outline the efficiency of the proposed computational model (CDM), the Authors considered the derivation of a beam finite element for static analysis of a cracked double-beam system. In particular, a finite element of a cracked double-beam system has been developed on the basis of a variational approach and the shape functions for rotational and translational displacements have been used to develop the stiffness matrix in presence and absence of cracks.

For the structure under consideration (see Figure 4), the total potential energy is given by:

$$\begin{array}{cc}\hfill E_t=& \frac{1}{2}{\int }_0^{\mathrm{L}}{\mathrm{EI}}_1{\left(\frac{{\partial }^2{\mathrm{v}}_1\left(\mathrm{z}\right)}{\partial {\mathrm{z}}^2}\right)}^2\mathrm{dz}+{\int }_0^{\mathrm{L}}{\mathrm{EI}}_2{\left(\frac{\partial \mathrm{v}\left(\mathrm{z}\right)}{\partial {\mathrm{z}}^2}\right)}^2\mathrm{dz}-{\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{v}}_1\left(\mathrm{z}\right)\mathrm{dz}+\hfill \\ \hfill & \frac{1}{2}{\int }_0^{\mathrm{L}}{k}_m{\left({\mathrm{v}}_2\left(\mathrm{z}\right)-{\mathrm{v}}_1\left(\mathrm{z}\right)\right)}^2\mathrm{dz}\hfill \end{array}$$

(54)

Using the cubic polynomial functions for the transverse displacements, the following shape functions are derived:

$${\mathrm{N}}_i=\left(\begin{array}{c}1-\frac{3{\mathrm{z}}^2}{{\mathrm{L}}^2}+\frac{2{\mathrm{z}}^3}{{\mathrm{L}}^3}\\ -\mathrm{z}+\frac{2{\mathrm{z}}^2}{\mathrm{L}}-\frac{{\mathrm{z}}^3}{{\mathrm{L}}^2}\\ \frac{3{\mathrm{z}}^2}{{\mathrm{L}}^2}-\frac{2{\mathrm{z}}^3}{{\mathrm{L}}^3}\\ \frac{{\mathrm{z}}^2}{\mathrm{L}}-\frac{{\mathrm{z}}^3}{{\mathrm{L}}^2}\end{array}\right);$$

(55)

where i=1,2 refers to the upper and lower beams, respectively. The transverse displacements are thus given as:

$${\mathrm{v}}_i\left(\mathrm{z}\right)={\mathrm{N}}_i^T{\mathrm{d}}_i$$

(56)

where ${\mathrm{d}}_i$ denotes the vector of nodal displacements and is given by:

$${\mathrm{d}}_i={\left({\mathrm{V}}_1^{\left(i\right)},{\Theta }_1^{\left(i\right)},{\mathrm{V}}_2^{\left(i\right)},{\Theta }_2^{\left(i\right)}\right)}^T$$

(57)

Deriving the nodal displacements respect to spatial coordinate z, the following expressions are obtained:

$$\frac{{\mathrm{dv}}_i\left(\mathrm{z}\right)}{\mathrm{dz}}={\mathrm{N}}_i^{{ }^{\prime }T}{\mathrm{d}}_i;\frac{{\mathrm{d}}^2{\mathrm{v}}_i\left(\mathrm{z}\right)}{{\mathrm{dz}}^2}={\mathrm{N}}_i^{’’T}{\mathrm{d}}_i$$

(58)

which substituted into Equation (54) leads to:

$$\begin{array}{cc}\hfill E_t=& \frac{1}{2}{\mathrm{d}}_1^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1^{’’}\mathrm{EI}{\mathrm{N}}_1^{’’T}\mathrm{dz}\right){\mathrm{d}}_1+\frac{1}{2}{\mathrm{d}}_2^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2^{’’}\mathrm{EI}{\mathrm{N}}_2^{’’T}\mathrm{dz}\right){\mathrm{d}}_2-{\mathrm{d}}_1^T{\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{N}}_1\mathrm{dz}+\hfill \\ \hfill & \frac{1}{2}{\mathrm{d}}_1^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1+\frac{1}{2}{\mathrm{d}}_2^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2+\hfill \\ \hfill & -{\mathrm{d}}_1^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2-{\mathrm{d}}_2^T\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1\hfill \end{array}$$

(59)

Applying the principle of stationary potential energy, one gets:

$$\begin{array}{cc}\hfill & \left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1^{’’}\mathrm{EI}{\mathrm{N}}_1^{’’T}\mathrm{dz}\right){\mathrm{d}}_1+\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1-\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2={\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{N}}_1\mathrm{dz}\hfill \\ \hfill & \left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2^{’’}\mathrm{EI}{\mathrm{N}}_2^{’’T}\mathrm{dz}\right){\mathrm{d}}_2+\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2-\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1=\mathbf{0}\hfill \end{array}$$

(60)

rearrangement the terms in ${\mathbf{d}}_1$ and ${\mathbf{d}}_2$, one gets:

$$\begin{array}{cc}\hfill & \left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1^{’’}\mathrm{EI}{\mathrm{N}}_1^{’’T}\mathrm{dz}+{\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1-\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2={\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{N}}_1\mathrm{dz}\hfill \\ \hfill & {\int }_0^{\mathrm{L}}{\mathrm{N}}_2^{’’}\mathrm{EI}{\mathrm{N}}_2^{’’T}\mathrm{dz}+\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_2^T\mathrm{dz}\right){\mathrm{d}}_2-\left({\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_1^T\mathrm{dz}\right){\mathrm{d}}_1=\mathbf{0}\hfill \end{array}$$

(61)

or

$$\left(\begin{array}{cc}{\mathrm{K}}_{11}& -{\mathrm{K}}_{12}\\ -{\mathrm{K}}_{21}& {\mathrm{K}}_{22}\end{array}\right)\left(\begin{array}{c}{\mathrm{d}}_1\\ {\mathrm{d}}_2\end{array}\right)=\left(\begin{array}{c}\mathrm{Q}\\ \mathbf{0}\end{array}\right)$$

(62)

being

$$\begin{array}{c}\hfill {\mathrm{K}}_{11}={\int }_0^{\mathrm{L}}{\mathrm{N}}_1^{’’}\mathrm{EI}{\mathrm{N}}_1^{’’T}\mathrm{dz}+{\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_1^T\mathrm{dz};\\ \hfill {\mathrm{K}}_{12}={\int }_0^{\mathrm{L}}{\mathrm{N}}_1{k}_m{\mathrm{N}}_2^T\mathrm{dz};\\ \hfill {\mathrm{K}}_{21}={\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_1^T\mathrm{dz};\\ \hfill {\mathrm{K}}_{22}={\int }_0^{\mathrm{L}}{\mathrm{N}}_2^{’’}\mathrm{EI}{\mathrm{N}}_2^{’’T}\mathrm{dz}+{\int }_0^{\mathrm{L}}{\mathrm{N}}_2{k}_m{\mathrm{N}}_2^T\mathrm{dz};\\ \hfill \mathrm{Q}={\int }_0^{\mathrm{L}}\mathrm{q}{\mathrm{N}}_1\mathrm{dz}.\end{array}$$

(63)

Equation (62) can be rewritten as:

$$\begin{array}{c}\hfill {\mathrm{K}}_{11}{\mathrm{d}}_1-{\mathrm{K}}_{12}{\mathrm{d}}_2=\mathrm{Q}\\ \hfill -{\mathrm{K}}_{21}{\mathrm{d}}_1+{\mathrm{K}}_{22}{\mathrm{d}}_2=\mathbf{0}\end{array}$$

(64)

From the second Equation (64) we obtain ${\mathbf{d}}_1$ which substituted in the first equation leads to:

$$\begin{array}{c}\hfill {\mathrm{d}}_2={\mathrm{K}}_{22}^{-1}{\mathrm{K}}_{21}{\mathrm{d}}_1;\\ \hfill \left({\mathrm{K}}_{11}-{\mathrm{K}}_{12}{\mathrm{K}}_{22}^{-1}{\mathrm{K}}_{21}\right){\mathrm{d}}_1=\mathrm{Q}\end{array}$$

(65)

or

$$\begin{array}{c}\hfill {\mathrm{d}}_1={\left({\mathrm{K}}_{11}-{\mathrm{K}}_{22}^{-1}{\mathrm{K}}_{12}{\mathrm{K}}_{21}\right)}^{-1}\mathrm{Q};\\ \hfill {\mathrm{d}}_2={\mathrm{K}}_{22}^{-1}{\mathrm{K}}_{12}{\mathrm{d}}_1\end{array}$$

(66)

In presence of a crack the stiffness matrix of finite element with a crack must be modified. Several investigators have been derived the stiffness matrix for a cracked beam element and the crack effect on the stiffness matrix has been investigated (see for example [34,35]). In the present analysis the structure under consideration has been analyzed implementing beam finite elements that have not taken the effect of shear forces into account. Therefore the stiffness matrix, which represents the relationship between nodal force vector and nodal displacements vector, can be expressed as:

$$\begin{array}{c}\hfill {\mathbf{K}}_{\mathbf{c}}=\frac{\mathrm{EI}}{\delta }\left(\begin{array}{cccc}12a_{11}& -6a_{12}& -12a_{13}& -6a_{14}\\ -6a_{21}& a_{22}& 6a_{23}& a_{24}\\ -12a_{31}& 6a_{32}& 12a_{33}& 6a_{34}\\ -6a_{41}& a_{42}& 6a_{43}& a_{44}\end{array}\right)\end{array}$$

(67)

so as deduced in [35] (see formula 22), where $k_r=\frac{\mathrm{EI}}{\Psi }$, $\delta$ is given by:

$$\delta =L\left(k_r{\mathrm{L}}^3+\mathrm{EI}\left(12{\mathrm{L}}_1^2-12\mathrm{L}{\mathrm{L}}_1+{\mathrm{L}}^{2}4\right)\right)$$

(68)

and where

$$\begin{array}{c}\hfill a_{11}=a_{33}=\left(\mathrm{EI}+k_r\mathrm{L}\right);\\ \hfill a_{22}=\left(12\mathrm{EI}{\mathrm{L}}_1^2+k_r{\mathrm{L}}^{3}4\right);\\ \hfill a_{44}=\left(12\mathrm{EI}{\left(\mathrm{L}-{\mathrm{L}}_1\right)}^2+4k_r{\mathrm{L}}^3\right);\\ \hfill a_{12}=a_{21}=\left(k_r{\mathrm{L}}^2+2\mathrm{EI}{\mathrm{L}}_1\right);\\ \hfill a_{13}=a_{31}=\left(\mathrm{EI}+k_r\mathrm{L}\right);\\ \hfill a_{14}=a_{41}=\left(k_r{\mathrm{L}}^2+2\mathrm{EI}\left(\mathrm{L}-{\mathrm{L}}_1\right)\right);\\ \hfill a_{23}=a_{32}=\left(k_r{\mathrm{L}}^2+2\mathrm{EI}{\mathrm{L}}_1\right);\\ \hfill a_{24}=a_{42}=\left(12\mathrm{EI}\left(\mathrm{L}-{\mathrm{L}}_1\right){\mathrm{L}}_1+2k_r{\mathrm{L}}^3\right);\\ \hfill a_{34}=a_{43}=\left(k_r{\mathrm{L}}^2+2\mathrm{EI}\left(\mathrm{L}-{\mathrm{L}}_1\right)\right).\end{array}$$

(69)

This matrix will be translated into ${\mathrm{K}}_{22}$ at the appropriate abscissa.

4. Numerical examples: results and discussion

In order to evaluate the effects of crack, taper ratio, elastic medium parameter, boundary conditions and rotational and traslational stiffness on the two beams deflection, several numerical examples are presented, using a general code developed in Mathematica [36]. The solutions in terms of deflection function are derived by using analytical and numerical approaches. The results obtained using the Euler-Bernoulli beam model are also presented to show the validity and reliability of the proposed approach (CDM).

In all numerical examples, the structure under consideration is composed of two parallel beam having the same length L equal to 2.5 m and the following material properties: Young’s modulus is 30 GPa and Poisson’s ratio $\nu$=0.3. The cross-section of two beams is a rectangular with base and height dimensions of b=0.10 m and height 0.18 m, respectively. Also, let be $k_m$=1 MPa the media-type foundation modulus and q=100 kN/m a uniformly distributed load applied to the upper beam.

4.1. Case 1: Simply supported-simply supported double-beam system in presence of a crack, on the lower beam, and a uniformly distributed load applied to the upper beam

In the first numerical example, static deflections of a simply supported-simply supported (SS-SS) double-beam system, in presence of a crack and a uniformly distributed load, has been investigated, by using a large number of cells, i.e. n = 300. The crack is introduced at the dimensionless abscissa z=1.25 m from the left end of the lower beam. The non-dimensional crack height is assumed to be $\xi$=0.5. The upper beam is loaded with a uniformly distributed load q=100 kN/m.

Figure 5 displays the dimensionless transverse displacements of two beams. Particularly, the dashed line refers to the vertical displacements of the upper beam, whereas the solid line shows the displacements of the lower beam. As can be seen the maximum value of transverse displacement ${\mathrm{v}}_1$ is equal to 0.0291243 m, at the dimensionless abscissa z=1.25 m, whereas the transverse displacement ${\mathrm{v}}_2$ at the crack position (z=1.25 m) is equal to 0.00890037. In order to check the correctness of the numerical calculations of CDM, a numerical comparison with the results given by closed-form solutions is proposed. Applying the exact procedure the tranverse displacements values of two beams are ${\mathrm{v}}_1$(z=1.25 m)=0.0291244 m and ${\mathrm{v}}_2$ (z=1.25 m)=0.00890037 m, respectively. As can be seen, the deflections of two beams show that there is an excellent agreement among the results obtained by the theoretical and numerical procedures. Also, in Figure 4 is shown the influence of the crack position on the deflection of SS−SS two beams. This figure emphasizes that the deflection will be maximized if the position of crack approaches the midspan of the beam.

In order to show the robustness,the accuracy and reliability of the proposed numerical method, the structure under consideration has also been analyzed implementing beam finite element.

In a first case the structure has been analyzed implementing beam finite element that takes the effect of a crack and distributed load into account. The crack is introduced in the midspan of the lower beam ($L_1=0.5L$) and a uniformly distributed load is applied to the upper beam only. For different values of Winkler-type medium stiffness $k_m$, the maximum displacements of the beam system have been calculated by means Finite element method (FEM), closed-form solution and CD method and the results obtained are quoted in Table 1. As can be easily seen from the table, using the proposed analytical approach and numerical methods, the estimated values of maximum displacements of double beam system are in excellent agreement with numerical calculations based on the FE method. Also, the results obtained show that the displacements decrease as $k_m$ increases.

In the second case, the structure has been analyzed implementing beam finite element that have not taken the effect of crack into account. Particularly, a simply supported-simply supported (SS-SS) double-beam system submitted to a uniformly distributed load, acting on the upper beam only, has been considered. For different values of $k_m$ the maximum displacements of the beam system have been calculated and the results obtained from the closed-form procedure and numerical methods based on the Cell-Discretization Method (CDM) and Finite Element method (FEM) are compared and quoted in Table 2.

As shown from the table, the estimated values of maximum displacements of double beam system are in excellent agreement with numerical calculations based on the FE method. Also, the results obtained show that the displacements decrease as $k_m$ increases.

4.2. Case 2: Simply supported-simply supported double-beam system in presence of a uniformly ditributed load, on the upper beam, and of a crack, on the upper and lower beams

In the second example, a simply supported-simply supported (SS-SS) double-beam system with two transverse cracks has been considered. The transverse cracks are applied to the upper and lower beams at the non-dimensional distances of z=0.5 m and z=1.75 m, respectively. The upper beam is loaded with a uniformly distributed load q=100 kN/m.

In Figure 6 the deflections of two beams have been reported. Particularly, the vertical displacements of the upper beam are represented with a dashed line, whereas the solid line shows the displacements of the lower beam. As can be noted, for the upper beam the maximum value of transverse displacement ${\mathrm{v}}_1$ is equal to 0.0321605 m, at the dimensionless abscissa z=1.18 m, whereas, for the lower beam, the transverse displacement ${\mathrm{v}}_2$ at the crack position (z=1.2375 m) is equal to 0.00847839. The closed-form solution, on the other hand, provides the following results: the maximum value of deflection ${\mathrm{v}}_1$(z=1.18 m)=0.0321452 m, for the upper beam, and the maximum value of deflection ${\mathrm{v}}_2$(z=1.2375 m)=0.0084943 m, for lower beam. It can be seen that the predicted results of CDM are in excellent agreement the closed-form solutions.

4.3. Case 3: Effect of elastic constraints and crack on the static behaviour of double-beam system

This example considers a double-beam system with the same geometrical and material properties of the previous numerical examples. The beams are elastically restrained against translational and rotational at either ends, with the following fixed values $k_{1\mathrm{TL}}$$=k_{2\mathrm{TL}}$$=k_{1\mathrm{TR}}$= $k_{2\mathrm{TR}}={10}^{10}$ and varying the non-dimensional rotational stiffness in the range [0-$\left.{10}^4\right]$. The transverse crack is located at the mid-span of the lower beam. A uniformly distributed load q=100 kN/m is applied to the upper beam. Table 3 presents the deflections first value for simply supported double beam system, for different values of the non-dimensional stiffnesses. In Table 3, results involving respective maximum transverse displacements ${\mathrm{v}}_{1max}$ and ${\mathrm{v}}_{2max}$, at the midspan for the upper and lower beams, are presented. As can be noted, as the stiffnesses values increase, deflections decrease throughout the two beams, especially near the crack. An analogy between increasing the stiffness and decreasing the deflections is illustrated in Figure 7 which represents the clamped-clamped beam case.

4.4. Case 4: Effect of taper ratio coefficient and crack on the static behaviour of double-beam system

This numerical example considers the influence of taper ratio and crack on the static behaviour of two non-uniform double-beam system, having the following law of moment of inertia:

$$\mathrm{I}\left(z\right)={\mathrm{I}}_0{\left((\beta -1)\frac{z}{\mathrm{L}}+1\right)}^4$$

(70)

with $\beta <$1 cross-section decreases with the abscissa, on the contrary if $\beta >$1 cross-section increaes with the abscissa. A uniformly distributed load, q=100 kN/m, is applied to the upper beam and a crack is introduced at the midspan of the lower beam. For varying $\beta$ in the range [0.5-2], the effect of taper ratio is evaluated as obtained by using CDM. In Figure 8 and Figure 9, the displacements diagrams are depicted, for $\beta$=0.5, 1, respectively. The upper beam is simply supported at the left end and clamped at the right end, while the lower beam is simply supported at either ends.

For $\beta$=1.5,2, in Figure 10 and Figure 11 the displacements diagrams are plotted. The upper beam is simply-supported at either ends, while the lower beam is simply-supported, at the left end, and clamped, at the right end.

For varying $\beta$=0.5-2, in Table 4 the maximum deflections values for two beams are quoted. It is interesting to note that the maximum deflections value decreases when the taper ratio coefficient increases.

4.5. Case 5: Effect of a concentrated force and crack located at the midspan of upper and lower beams, respectively

In this numerical example, the effect of concentrated force is evaluated. Consider a simply supported-simply supported (SS-SS) double-beam system in presence of a crack. The transverse crack is located at the midspan of the lower beam and a concentrated force F is located at the midspan of the upper beam. Table 5 shows the maximum deflections values for different values of the concentrated force and in absence of a distributed load. As can be noted, as the value of the concentrated force F increases, the values of displacements at the midspan of the two beams increase.

For the same numerical example, fixed the force value F=${10}^4$ N and varying the nondimensional height of crack $\xi =\frac{a}{h}$, located at the mid-distance of the lower beam, the maximum deflections are calculated and listed in Table 6. It is interesting to note that the maximum deflections value increases as the nondimensional height of crack increases.

4.6. Case 6: Effect of the slenderness $\lambda$ on the deflections for both upper and lower beams in presence of a crack

In this numerical example, to study the effects of the slenderness $\lambda =\frac{L}{h}$ on displacements and rotations, considers a simply-supported double beam system having the same geometrical and material properties of the previous numerical examples. The beam system is submitted to a uniformly distributed load q=100 kN/m acting on the upper beam only. A crack is introduced in the midspan of the lower beam.

Setting $\Gamma =\frac{h}{b}=1.8$ and being $\lambda =\frac{L}{h}$, the expressions of the parameters of Eq. (13) became:

$$\alpha =\lambda \sqrt[4]{\frac{12k_m\Gamma }{\mathrm{E}}};p=\frac{12\mathrm{q}\Gamma {\lambda }^4}{\mathrm{E}}$$

(71)

Eq. (71) denotes the relationship between the deflection and slenderness.

For different values of slenderness $\lambda$ of the two beams, the maximum displacements and the rotations at the left and right ends (for $\zeta$=0 and $\zeta$=1), respectively, are calculated. The results obtained from the closed-form solution and numerical method based on the Cell-Discretization Method (CDM) are compared and listed in Table 7. As can be easily seen from the table, the displacements values are in excellent agreement. Also, the results obtained show that the displacements increase as the slenderness $\lambda =\frac{L}{h}$ of the two beams increases.

5. Concluding remarks

In the present paper, the static behaviour of a double-beam system, in the presence of Winkler medium, carrying a crack, at generic position, on the lower beam, and subjected to a distributed load, on the upper beam, have been studied. The double-beam system is also supposed to be constrained at the ends by elastically flexible springs, with transverse stiffness and rotational stiffness. According to the Euler-Bernoulli beam theory, the static governing equations are derived using an variational formulation and have been solved through implementing analytical and numerical approaches. Among numerical approach, the Cell-Discretization Method (CDM) is employed to solve the governing equations, in which the beam is reduced to a set of rigid bars linked together by means of elastic constraints. A comparative analysis has been performed in order to verify the accuracy and validity of the proposed numerical method. The effects of crack, taper ratio, elastic medium parameter, boundary conditions and rotational and traslational stiffness on the two beams deflections through theoretical and numerical approaches have been presented. Through the obtained results, the following observations apply:

-

if the transverse stiffness value increase the maximum deflections value decreases;

-

if the taper ratio coefficient increases the maximum deflections value decreases;

-

the crack plays a key role for static behaviour of double-beam system: the maximum deflections value increases when the crack increases.

Finally, the numerical examples demonstrate that the results determined by CDM perfectly match with the solutions from the governing equations of the double-beam model and are in good agreement with the results obtained by FEM. It is shown that the CD method has a good and rapid convergence regardless of the beam theory, crack and elastic foundation parameters.