A peridynamic formulation for planar arbitrarily curved beams with Timoshenko

Abstract

In this study, a peridynamic formulation is proposed for analysis of planar arbitrarily curved beams. Assumptions of Timoshenko–Ehrenfest beam model are considered such that linear kinematic relations are described by displacements of the beam axis and rotation of the cross-section. Equilibrium conditions and boundary conditions are derived by means of the principle of virtual work. Then, peridynamic functions are constructed, and incorporated with equilibrium conditions to develop a peridynamic formulation. Several examples are exhibited to examine the performance of the proposed formulation in some regards, i.e., numerical accuracy, convergence properties, and locking effects. It is delineated that the proposed formulation provides a good level of precision for the analysis of relatively thick beams. However, membrane and shear locking effects are present in the analysis of slender beams, and peridynamic functions of higher polynomial degree can be used to mitigate these unfavorable effects.

Keywords

Planar arbitrarily curved beams Timoshenko–Ehrenfest beam model peridynamic beam formulation peridynamic differential operators locking effects

1. Introduction

In classical continuum mechanics, equations of motion are written as partial differential equations, which are not well-defined with the presence of discontinuities. Consequently, peridynamic theory [1 –5] has been introduced as a reformulation of classical continuum mechanics to eliminate their roadblocks in simulation of damage of solid bodies. Indeed, equations of motion in the peridynamic theory are expressed as integro-differential equations, which are valid on discontinuities. Furthermore, the peridynamic theory can be classified as a nonlocal one, meaning that a material particle interacts with other counterparts within a finite distance. These two distinct features make the peridynamic theory suitable for solving problems with evolving discontinuities and significant nonlocality [6 –15].

In the implementation of the peridynamic theory, each material particle requires the construction of its own equations of motion, which involves its own unknown variables as well as those of other particles it interacts with. Consequently, the obtained system of equations is normally large and non-symmetric, and computational costs of the peridynamic theory become prohibitively expensive for the analysis of three-dimensional objects. Indeed, for the analysis of slender structures such as beams or plates, several layouts of particles along the direction of thickness are essential to ensure precise analysis [16] if these structures are simulated as fully three-dimensional objects. In this scenario, it is beneficial to develop formulations by incorporating dimension-reduced structural models with the peridynamic theory for the simulation of slender three-dimensional objects.

Regarding beam models, Timoshenko–Ehrenfest considers shear deformation by decoupling the deformation mechanisms of the beam axis and the cross-section. More precisely, the rigid cross-section is not necessarily always orthogonal to the beam axis during deformation. Hence, the Timoshenko–Ehrenfest beam model is widely used since it is capable of simulating both thick and thin beams. Recently, several beam formulations [17 –26] have been developed in the context of peridynamic theory using kinematic assumptions of Timoshenko–Ehrenfest beam model. The common approach to derive the peridynamic equations of motion is first to construct the nonlocal expressions of the kinetic and total potential energy, and then use the Euler–Lagrange equation to obtain the equations of motion for each particle. This approach has been performed in previous works [17,18,21,24,25] for the geometrically linear analysis and in Nguyen and Oterkus [22] for the geometrically nonlinear analysis. The incorporation of peridynamic beam formulations with the finite element framework is carried out by Jafari et al. [19] and Yang et al. [20]. An alternative approach is to use the existing governing equations derived from the classical continuum mechanics with peridynamic differential operator [26], which converts the local derivatives into the nonlocal integration.

It is well known that the numerical precision of Timoshenko–Ehrenfest beam formulations can be degraded in analysis of slender beams owing to membrane and shear locking effects. These unfavorable effects have been observed with some other numerical methods, e.g., finite element methods [27 –29], isogeometric analysis [30 –36], and mesh-free methods [37]. For Timoshenko–Ehrenfest beam formulations with peridynamic theory, the presence of membrane and shear locking effects is also confirmed [24,25].

With the publications reviewed above, it is recognized that a locking-free peridynamic Timoshenko–Ehrenfest beam formulation for analysis of beams having arbitrary geometry has not yet been introduced. In this study, a peridynamic Timoshenko–Ehrenfest beam formulation is presented for the analysis of planar arbitrarily curved beams, and locking effects on numerical performance of the proposed formulation are examined. More precisely, with assumptions of Timoshenko–Ehrenfest beam model, displacements of the beam axis and the rotation of cross-section are used to describe linear kinematic relations. Equilibrium conditions together with boundary conditions are derived by means of the principle of virtual work. Then, a peridynamic formulation is developed by incorporating peridynamic functions with derived equilibrium conditions.

The rest of this paper is structured as follows. Section 2 involves the derivation of equilibrium conditions and boundary conditions of planar arbitrarily curved beams using kinematic assumptions of Timoshenko–Ehrenfest beam model. Section 3 presents the detailed construction of peridynamic functions using the concept of peridynamic differential operators [38,39], and the proposed peridynamic formulation is detailed. In section 4, three examples with several rigorous tests are designed to elucidate the performance of the proposed formulation. Conclusions are given in section 5.

2. Equilibrium conditions and boundary conditions of planar curved beams with Timoshenko–Ehrenfest beam model

In this section, linear kinematic relations of planar curved beams obeying assumptions of Timoshenko–Ehrenfest beam model are briefly detailed. Then, generalized strain and stress measures are presented for use of the principle of virtual work to derive equilibrium conditions and boundary conditions of planar arbitrarily curved beams.

2.1. Linear kinematic relations

A representative planar curved beam is pictorially portrayed in Figure 1. In the following derivations, the beam axis is assumed to be a regular curve, and parameterized by the arc length parameter $S$ , i.e., $R (S) : [0, L] \to R^{2}$ , with $L$ representing the length of the beam. At a material point located on the beam axis, the unit tangent and normal vectors are identified as follows:

A_{1} = \frac{d R (S)}{d S} = R^{'} (S), A_{2} = Λ A_{1}

(1)

where $Λ = [\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}]$ poses a rotation of $90 °$ in the counterclockwise direction. Throughout this study, primes denote derivatives with respect to the arc length parameter. Since the vectors $A_{1}$ and $A_{2}$ are orthonormal, they are used as basic vectors of a local Cartesian coordinate system moving along the beam axis.

Figure 1.

A representative planar arbitrarily curved beam.

In the context of Timoshenko–Ehrenfest beam model, the rigid cross-section does not necessarily remain orthogonal to the beam axis during deformation. Hence, linear kinematic relations of planar curved beams are described by the displacement vector $u$ of the beam axis and the rotation $ϕ$ of the cross-section. Here, the displacement vector $u$ and the cross-sectional rotation $ϕ$ are decoupled. In the local Cartesian coordinate system, the displacement vector $u$ is expressed as follows:

u (S) = u_{1} A_{1} + u_{2} A_{2}

(2)

where $u_{1}$ and $u_{2}$ can be interpreted as the displacement components of the beam axis along the vectors $A_{1}$ and $A_{2}$ , respectively.

The deformation state of the beam is characterized by the membrane strain $ε$ , the bending curvature $χ$ , and the shear strain $γ$ as [40] follows:

ε = u_{1}' - κ u_{2}, χ = - ϕ', γ = κ u_{1} + u'_{2} - ϕ

(3)

with $κ$ being the signed curvature of the beam axis. Here, the membrane strain $ε$ describes the axial deformation, the bending curvature $χ$ measures the change in the curvature of the beam axis, and the shear strain $γ$ represents the shear deformation.

In this study, only linear elastic isotropic materials are considered, and consequently, the following constitutive relations are used as below:

N = EA ε, M = EI χ, T = GA γ

(4)

where $N$ , $M$ , and $T$ are, respectively, the axial force, the bending moment, and the shear force. For material properties, $E$ and $G$ , respectively, denote Young’s modulus and the shear modulus. In addition, the cross-sectional area and moment of inertia are, respectively, indicated as $A$ and $I$ .

2.2. Equilibrium conditions and boundary conditions

For planar curved beams, the internal virtual work $δ W_{int}$ and the external virtual work $δ W_{ext}$ can be written as follows:

δ W_{i n t} = \int_{0}^{L} (N δ ε + M δ χ + T δ γ) d S

(5)

δ W_{e x t} = \int_{0}^{L} (f_{1} δ u_{1} + f_{2} δ u_{2} + m δ ϕ) d S + {(F_{1} δ u_{1} + F_{2} δ u_{2} + Q δ ϕ) |}_{S = 0}^{S = L} .

(6)

Here, $f = f_{1} A_{1} + f_{2} A_{2}$ is the distributed force, whereas $m$ stands for the distributed moment. In addition, $F = F_{1} A_{1} + F_{2} A_{2}$ is the concentrated force, and $Q$ is the concentrated moment. These concentrated loads are prescribed at the ends of the beam.

The principle of virtual work is now expressed as follows:

δ W_{int} - δ W_{ext} = 0 .

(7)

With the aid of Equations (3)–(7), one can arrive at the equilibrium conditions together with boundary conditions by employing the integration by parts and fundamental theorem of variational calculus, as follows.

Equilibrium conditions:

- N' + κ T = f_{1}

(8)

- κ N - T^{'} = f_{2}

(9)

M' - T = m .

(10)

Boundary conditions:

N = F_{1} or u_{1} = {\bar{u}}_{1}

(11)

T = F_{2} or u_{2} = {\bar{u}}_{2}

(12)

- M = Q or ϕ = \bar{ϕ} .

(13)

The overlines signify prescribed values of the displacement components and cross-sectional rotation at the boundary.

The primal forms of the equilibrium conditions and boundary conditions are obtained by substituting Equations (3) and (4) into Equations (8)–(13) as below.

Primal forms of the equilibrium conditions:

- {EAu ″}_{1} + κ^{2} GA u_{1} + κ (EA + GA) u'_{2} + κ' EA u_{2} - κ GA ϕ = f_{1}

(14)

- κ (E A + G A) {u^{'}}_{1} - κ^{'} G A u_{1} - G A {u^{″}}_{2}^{+} κ^{2} E A u_{2} + G A ϕ^{'} = f_{2}

(15)

- κ GA u_{1} - GAu'_{2} - EI ϕ ″ + GA ϕ = m .

(16)

Primal forms of the boundary conditions:

EA (u'_{1} - κ u_{2}) = F_{1} or u_{1} = {\bar{u}}_{1}

(17)

GA (κ u_{1} + u'_{2} - ϕ) = F_{2} or u_{2} = {\bar{u}}_{2}

(18)

EI ϕ^{'} = Q or ϕ = \bar{ϕ} .

(19)

3. Peridynamic formulation for planar arbitrarily curved beams

This section discusses the construction of univariate peridynamic functions. Then, a peridynamic formulation for the analysis of planar curved beams is derived by means of the peridynamic functions.

3.1. Discretization and numerical integration scheme

Prior to detailed construction of the peridynamic functions, a quadrature scheme for numerical integration is discussed. Figure 2 exhibits a specific discretization of a beam axis with $N_{P}$ material particles, each with specific entities. For a particle $k$ , its location on the beam axis is identified by an arc length coordinate $S_{(k)}$ , and the length of the beam occupied by this material particle is denoted by $L_{(k)}$ . The displacement components of the material particle $k$ along the vectors $A_{1}$ and $A_{2}$ are, respectively, denoted by $u_{1 (k)}$ and $u_{2 (k)}$ . In addition, the cross-sectional rotation at the particle $k$ is represented by $ϕ_{(k)}$ .

Figure 2.

A peridynamic discretization.

A particle $j$ is called a family member of particle $k$ if the following condition holds as below:

| S_{(j)} - S_{(k)} | \leq δ_{(k)}

(20)

where $δ_{(k)}$ is a positive scalar referred to as the horizon of particle $k$ . All family members encompassed by the horizon $δ_{(k)}$ form the family $H_{(k)}$ , and the number of family members in $H_{(k)}$ is denoted by $N_{P (k)}$ . The global index of the $i$ th particle in $H_{(k)}$ can be extracted by the operator $H_{(k)} [i]$ . The length of the beam occupied by all family members in $H_{(k)}$ defines the interaction domain $D_{(k)}$ of the particle $k$ .

Given a scalar-valued function $f (S)$ , integration over the interaction domain $D_{(k)}$ can be performed by summing the entities of the members in family $H_{(k)}$ , as follows:

\int_{D_{(k)}} [f (S_{(k)} + Δ S) - f (S_{(k)})] d S = \sum_{i = 1}^{N_{P (k)}} [f (S_{(H_{(k)} [i])}) - f (S_{(k)})] L_{(H_{(k)} [i])} .

(21)

It is worth noting that the quadrature scheme mentioned above is not limited to uniform distributions of material particles. In practice, it is beneficial to employ discretization with nonuniform distributions of material particles for the analysis of beams having arbitrary geometry. Furthermore, the horizon is not necessarily a constant quantity, meaning that different values of horizon can be assigned to different particles. As illustrated in Figure 2, the particles $k$ and $j$ are, respectively, highlighted in red and magenta. These particles are assigned different values of horizon, i.e., $δ_{(k)} \neq δ_{(j)}$ . The family members of particle $k$ are encompassed by dashed red lines, whereas those of particle $j$ are surrounded by dashed–dotted magenta lines. One can easily recognize that particle $j$ is a family member of particle $k$ , but not vice versa.

3.2. Peridynamic functions

In the vicinity of a material particle $k$ , the Taylor series expansion of a univariate function $f (S)$ is written as follows:

f (S_{(k)} + Δ S) = f (S_{(k)}) + \sum_{p = 1}^{P} \frac{1}{p!} {\frac{d^{p} f (S)}{d S^{p}} |}_{S_{(k)}} {(Δ S)}^{p} + R [{(Δ S)}^{P}]

(22)

where $R [{(Δ S)}^{P}]$ is the remainder, which is considered negligibly small.

At particle $k$ , a set of peridynamic functions $g_{P (k)}^{q} (Δ S)$ is constructed posing the orthogonality property as follows:

\frac{1}{n!} \int_{D_{(k)}} {(Δ S)}^{n} g_{P (k)}^{q} (Δ S) d S = δ_{nq}, (n = 1, 2, \dots, P)

(23)

with $δ_{nq}$ representing Kronecker delta.

With the aid of Taylor series expansion and peridynamic functions, the derivatives of any order can be evaluated at particle $k$ as follows:

{\frac{d^{q} f (S)}{d S^{q}} |}_{S_{(k)}} = \int_{D_{(k)}} [f (S_{(k)} + Δ S) - f (S_{(k)})] g_{P (k)}^{q} (Δ S) d S .

(24)

The remaining issue is the explicit forms of peridynamic functions. Although not a limitation, the peridynamic functions are commonly constructed as polynomials, i.e.,

g_{P (k)}^{q} (Δ S) = \sum_{p = 1}^{P} a_{p}^{q} w_{p} (| Δ S |) {(Δ S)}^{p}

(25)

where $a_{p}^{q}$ is an unknown coefficient, and $w_{p} (| Δ S |)$ is a weight function associated with the term ${(Δ S)}^{p}$ in the polynomial expansion. For simplicity, an identical weight function can be selected for all polynomial terms. By substituting Equation (25) into Equation (23), a system of linear equations with respect to unknown coefficients $a_{p}^{q}$ is achieved as follows:

\sum_{p = 1}^{P} (\int_{D_{(k)}} w_{p} (| Δ S |) {(Δ S)}^{p + n} dS) a_{p}^{q} = n! δ_{nq}, (n = 1, 2, \dots, P) .

(26)

The determination of the unknown coefficients $a_{p}^{q}$ completes the construction of the peridynamic functions.

3.3. Peridynamic formulation for planar curved beams

By incorporating the peridynamic functions into Equations (14)–(16), the equilibrium conditions of a particle $k$ can be obtained as follows:

\begin{matrix} - EA \int_{D_{(k)}} [u_{1} (S_{(k)} + Δ S) - u_{1} (S_{(k)})] g_{P (k)}^{2} (Δ S) d S + κ^{2} (S_{(k)}) GA u_{1} (S_{(k)}) \\ + κ (S_{(k)}) (EA + GA) \int_{D_{(k)}} [u_{2} (S_{(k)} + Δ S) - u_{2} (S_{(k)})] g_{P (k)}^{1} (Δ S) d S \\ + κ^{'} (S_{(k)}) EA u_{2} (S_{(k)}) - κ (S_{(k)}) GA ϕ (S_{(k)}) = f_{1} (S_{(k)}) \end{matrix}

(27)

\begin{array}{l} - κ (S_{(k)}) (E A + G A) \int_{D_{(k)}} [u_{1} (S_{(k)} + Δ S) - u_{1} (S_{(k)})] g_{P (k)}^{1} (Δ S) d S - κ^{'} (S_{(k)}) G A u_{1} (S_{(k)}) \\ - G A \int_{D_{(k)}} [u_{2} (S_{(k)} + Δ S) - u_{2} (S_{(k)})] g_{P (k)}^{2} (Δ S) d S + κ^{2} (S_{(k)}) E A u_{2} (S_{(k)}) \\ + G A \int_{D_{(k)}} [ϕ (S_{(k)} + Δ S) - ϕ (S_{(k)})] g_{P (k)}^{1} (Δ S) d S = f_{2} (S_{(k)}) \end{array}

(28)

\begin{matrix} - κ (S_{(k)}) G A u_{1} (S_{(k)}) - G A \int_{D_{(k)}} [u_{2} (S_{(k)} + Δ S) - u_{2} (S_{(k)})] g_{P (k)}^{1} (Δ S) d S \\ - EI \int_{D_{(k)}} [ϕ (S_{(k)} + Δ S) - ϕ (S_{(k)})] g_{P (k)}^{2} (Δ S) d S + GA ϕ (S_{(k)}) = m (S_{(k)}) . \end{matrix}

(29)

Furthermore, in case the loads are prescribed, the boundary conditions in Equation (17)–(19) can be alternatively written as follows:

EA {\int_{D_{(a)}} [u_{1} (S_{(a)} + Δ S) - u_{1} (S_{(a)})] g_{P (a)}^{1} (Δ S) dS - κ (S_{(a)}) u_{2} (S_{(a)})} = F_{1} (S_{(a)})

(30)

GA {κ (S_{(a)}) u_{1} (S_{(a)}) + \int_{D_{(a)}} [u_{2} (S_{(a)} + Δ S) - u_{2} (S_{(a)})] g_{P (a)}^{1} (Δ S) dS - ϕ (S_{(a)})} = F_{2} (S_{(a)})

(31)

EI \int_{D_{(a)}} [ϕ (S_{(a)} + Δ S) - ϕ (S_{(a)})] g_{P (a)}^{1} (Δ S) d S = Q (S_{(a)})

(32)

with $a$ being either $1$ or $N_{P}$ for the particles at the ends of the beam.

For numerical implementation, the integration involved in Equations (27)–(32) is performed by means of quadrature scheme presented in section 3.1. Consequently, the equilibrium conditions at particle $k$ can be recast as follows:

\begin{matrix} - EA \sum_{i = 1}^{N_{P (k)}} [u_{1 (H_{(k)} [i])} - u_{1 (k)}] g_{P (k)}^{2} (S_{(ki)}) L_{(H_{(k)} [i])} + κ_{(k)}^{2} GA u_{1 (k)} \\ + κ_{(k)} (EA + GA) \sum_{i = 1}^{N_{P (k)}} [u_{2 (H_{(k)} [i])} - u_{2 (k)}] g_{P (k)}^{1} (S_{(ki)}) L_{(H_{(k)} [i])} \\ + κ_{(k)}^{'} EA u_{2 (k)} - κ_{(k)} GA ϕ_{(k)} = f_{1 (k)} \end{matrix}

(33)

\begin{matrix} - κ_{(k)} (EA + GA) \sum_{i = 1}^{N_{P (k)}} [u_{1 (H_{(k)} [i])} - u_{1 (k)}] g_{P (k)}^{1} (S_{(ki)}) L_{(H_{(k)} [i])} - κ_{(k)}^{'} GA u_{1 (k)} \\ - GA \sum_{i = 1}^{N_{P (k)}} [u_{2 (H_{(k)} [i])} - u_{2 (k)}] g_{P (k)}^{2} (S_{(ki)}) L_{(H_{(k)} [i])} + κ_{(k)}^{2} EA u_{2 (k)} \\ + GA \sum_{i = 1}^{N_{P (k)}} [ϕ_{(H_{(k)} [i])} - ϕ_{(k)}] g_{P (k)}^{1} (S_{(ki)}) L_{(H_{(k)} [i])} = f_{2 (k)} \end{matrix}

(34)

\begin{matrix} - κ_{(k)} GA u_{1 (k)} - GA \sum_{i = 1}^{N_{P (k)}} [u_{2 (H_{(k)} [i])} - u_{2 (k)}] g_{P (k)}^{1} (S_{(ki)}) L_{(H_{(k)} [i])} \\ - EI \sum_{i = 1}^{N_{P (k)}} [ϕ_{(H_{(k)} [i])} - ϕ_{(k)}] g_{P (k)}^{2} (S_{(ki)}) L_{(H_{(k)} [i])} + GA ϕ_{(k)} = m_{(k)} . \end{matrix}

(35)

Here, $κ_{(k)}$ is the signed curvature of the beam axis measured at particle $k$ . In addition, $f_{(k)} = f_{1 (k)} A_{1 (k)} + f_{2 (k)} A_{2 (k)}$ and $m_{(k)}$ , respectively, denote the distributed force vector and moment applied at particle $k$ . The difference in the arc length coordinates of particle $k$ and its $i$ th family member is defined as $S_{(ki)} = S_{(H_{(k)} [i])} - S_{(k)}$ .

In an analogous fashion, the boundary conditions in Equations (30)–(32) are alternatively expressed as follows:

EA {\sum_{i = 1}^{N_{P (a)}} [u_{1 (H_{(a)} [i])} - u_{1 (a)}] g_{P (a)}^{1} (S_{(ai)}) L_{(H_{(a)} [i])} - κ_{(a)} u_{2 (a)}} = F_{1 (a)}

(36)

GA {κ_{(a)} u_{1 (a)} + \sum_{i = 1}^{N_{P (a)}} [u_{2 (H_{(a)} [i])} - u_{2 (a)}] g_{P (a)}^{1} (S_{(ai)}) L_{(H_{(a)} [i])} - ϕ_{(a)}} = F_{2 (a)}

(37)

EI \sum_{i = 1}^{N_{P (a)}} [ϕ_{(H_{(a)} [i])} - ϕ_{(a)}] g_{P (a)}^{1} (S_{(ai)}) L_{(H_{(a)} [i])} = Q_{(a)} .

(38)

4. Illustrative examples

This section presents three examples to assess the performance of the proposed peridynamic formulation for the analysis of planar arbitrarily curved beams. Some remarks regarding numerical implementation are given before detailed discussions. In all tests, Young’s modulus $E = 2 \times 10^{11}$ and Poisson’s ratio $ν = 0.3$ . The shear modulus is calculated as $G = E / [2 (1 + ν)]$ . For the family of each particle, $P + 1$ nearest particles on each side of the particle of interest are included, with $P$ being the highest order of derivatives retained in the Taylor series expansion. This strategy is suggested [39] to balance the accuracy and efficiency in constructing the peridynamic functions. In addition, the associated weight functions of polynomial terms of the peridynamic functions are identical, and specifically chosen as follows:

w_{p} (| S_{(ki)} |) = e^{{(- \frac{4 | S_{(ki)} |}{(P + 1) δ_{m}})}^{2}}

(39)

where $δ_{m} = L / (N_{p} - 1)$ is the average spacing between particles.

The enforcement of the equilibrium conditions at all particles leads to a system of linear equations as follows:

Ku = F

(40)

in which $K$ is the coefficient matrix arisen from the equilibrium conditions, $u$ contains displacement components and cross-sectional rotation of all particles, and $F$ is the loads applied at the particles. Regarding boundary conditions, the constraints of kinematic quantities are imposed directly on the particles located at the ends of the beam. The imposition of boundary conditions can be compactly expressed as follows:

Gu = C

(41)

where $G$ is the coefficient matrix occurred due to the boundary conditions, and $C$ includes the specified values of kinematic quantities. Eventually, the solutions of systems in Equations (40) and (41) can be sought by solving the following system as [38,39] below:

[\begin{matrix} K & G^{T} \\ G & 0 \end{matrix}] [\begin{matrix} u \\ λ \end{matrix}] = [\begin{matrix} F \\ C \end{matrix}] .

(42)

Here, $λ$ is the vector of Lagrange multipliers.

Regarding cross-sectional geometry, square cross-sections with the side length of $h$ is considered.

4.1. A straight beam subjected to sinusoidal force

The first example deals with the analysis of a straight beam. As visualized in Figure 3, a straight beam is fixed at both ends, and subjected to a sinusoidal force described as follows:

f_{2} = q \sin (\frac{2 π S}{L}) .

(43)

with $q = 45 EI / L^{2}$ . These settings allow straightforward derivations of analytical solutions for all kinematic quantities, including the displacement component $u_{2}$ , the cross-sectional rotation $ϕ$ , the bending moment $M$ , and the shear force $T$ , i.e.,

\begin{matrix} u_{2}^{a} = \frac{q L^{2}}{16 π^{4} EIGA (12 EI + GA L^{2})} [48 {(π EI)}^{2} \sin (\frac{2 π S}{L}) + EIGA L^{2} (12 + 4 π^{2}) \sin (\frac{2 π S}{L}) \\ + {(GA L^{2})}^{2} \sin (\frac{2 π S}{L}) - 4 π {(GA)}^{2} L S^{3} + 6 π {(GA)}^{2} L^{2} S^{2} - 2 π {(GA)}^{2} L^{3} S] \end{matrix}

(44)

ϕ^{a} = \frac{q L^{3}}{8 π^{3} EI (12 EI + GA L^{2})} [- (12 EI + GA L^{2}) \cos (\frac{2 π S}{L}) + 6 GA S^{2} - 6 GALS + 12 EI + GA L^{2}]

(45)

M^{a} = \frac{q L^{2}}{4 π^{3} (12 EI + GA L^{2})} [(12 π EI + π GA L^{2}) \sin (\frac{2 π S}{L}) + 6 GALS - 3 GA L^{2}]

(46)

T^{a} = \frac{qL}{2 π^{3} (12 EI + GA L^{2})} [(12 π^{2} EI + π^{2} GA L^{2}) \cos (\frac{2 π S}{L}) + 3 GA L^{2}] .

(47)

Figure 3.

A straight beam subjected to sinusoidal force. (a) Problem setup. (b) A peridynamic discretization.

In this example, discretization with uniform distributions of material particles is employed. In the first test, the side length of the cross-section is set as $L / h = 10$ . Figure 4 shows the results obtained by the proposed formulation using a discretization of $31$ particles and $P = 4$ . Instantly, one can recognize the good agreement between the numerical results and analytical solutions for all considered kinematic quantities.

Figure 4.

A straight beam subjected to sinusoidal force—numerical results with $L / h = 10$ . (a) Displacement component $u_{2}$ . (b) Cross-sectional rotation. (c) Bending moment. (d) Shear force.

The convergence properties of the proposed formulation are investigated with respect to different values of $P$ . It should be recalled that the highest order of derivatives occurred in the equilibrium conditions is $2$ , which is the possible minimum value of $P$ for the convergence tests. Hence, the convergence tests are carried out with $P = 2, 3, 4, 5, 6$ , and $7$ . The relative error of a specific kinematic quantity (■) is calculated by means of discrete $L 2$ norm and $H 1$ norm as follows:

L 2 = \frac{\sqrt{\sum_{k = 1}^{N_{P}} {[{(■)}_{(k)} - {(■)}_{(k)}^{r}]}^{2} L_{(k)}}}{\sqrt{\sum_{k = 1}^{N_{P}} {[{(■)}_{(k)}^{r}]}^{2} L_{(k)}}}

(48)

H 1 = \frac{\sqrt{\sum_{k = 1}^{N_{P}} {[{(■)}_{(k)} - {(■)}_{(k)}^{r}]}^{2} L_{(k)} + \sum_{k = 1}^{N_{P}} {[{(■^{'})}_{(k)} - {(■^{'})}_{(k)}^{r}]}^{2} L_{(k)}}}{\sqrt{\sum_{k = 1}^{N_{P}} {[{(■)}_{(k)}^{r}]}^{2} L_{(k)} + \sum_{k = 1}^{N_{P}} {[{(■^{'})}_{(k)}^{r}]}^{2} L_{(k)}}} .

(49)

Here, ${(■)}_{(k)}$ and ${(■)}_{(k)}^{r}$ , respectively, denote the numerical and reference values of the quantity (■) at particle $k$ . Analogously, ${(■^{'})}_{(k)}$ and ${(■^{'})}_{(k)}^{r}$ , respectively, represent the numerical and reference values for the derivatives of the quantity (■) at particle $k$ . Since analytical solutions of this example are available, they are treated as reference ones. Throughout this section, the convergence behavior in the displacement of the beam axis and the cross-sectional rotation is examined with $L 2$ norm, whereas the convergence property regarding the axial force, shear force, and bending moment is inspected with $H 1$ norm.

In the convergence tests, the relative errors are considered as the functions of the number of particles $N_{P}$ , i.e., $RelativeError = c {(N_{P})}^{- d}$ . Here, $c$ is a constant, whereas $d$ is referred to as the order of convergence. In the log–log scale, $d$ represents the slope of the convergence line. Convergence results with $L / h = 10$ are exhibited in Figures 5(a)–(d), respectively, for displacement of the beam axis, cross-sectional rotation, bending moment, and shear force. For all examined quantities, two common observations can be highlighted, i.e., the relative errors are reduced for higher values of $P$ , and the order of convergence is improved along with increasing $P$ . Here, it should be recalled that $P$ not only indicates the highest order of derivatives retained in the Taylor series expansion, but also the highest polynomial degree of peridynamic functions. This feature implies that both the accuracy and order of convergence are enhanced by using peridynamic functions of pairs of higher polynomial degree. Furthermore, it is found that $P$ and $P - 1$ are, respectively, the order of convergence for even and odd values of $P$ used.

Figure 5.

A straight beam subjected to sinusoidal force—convergence results with $L / h = 10 .$ (a) Displacement of the beam axis. (b) Cross-sectional rotation. (c) Bending moment. (d) Shear force.

To scrutinize the shear locking effect on the numerical performance of the proposed formulation, the convergence tests are performed with $L / h = 10^{3}$ , which represents extremely thin beams. Convergence results are presented in Figure 6. Again, six values of $P$ are considered, i.e., $2, 3, 4, 5, 6$ , and $7$ . For all examined kinematic quantities, a loss in order of convergence is recognized for $P = 2, 3, 4$ , and $5$ . For $P = 6$ and $7$ , the asymptotic order of convergence is only achieved with a sufficiently large value of $N_{P}$ . This behavior can be attributed to the presence of the shear locking effect.

Figure 6.

A straight beam subjected to sinusoidal force—convergence results with $L / h = 10^{3}$ . (a) Displacement of the beam axis. (b) Cross-sectional rotation. (c) Bending moment. (d) Shear force.

4.2. A semi-circular beam

The second example is devoted to the analysis of a semi-circular beam of radius $R$ . The geometry of the beam together with loading and boundary conditions are portrayed in Figure 7. Both ends of the beam are fixed, whereas uniform tangential and radial forces are applied on the beam axis as follows:

f_{1} = 50 \frac{EI}{R L^{2}}, f_{2} = 0.5 \frac{EA}{L} .

(50)

Figure 7.

A semi-circular beam. (a) Problem setup. (b) A peridynamic discretization.

Discretization with uniform distributions of material particles is used for all tests of this example. The analytical solutions can be attained with the procedure detailed in the Appendix.

The slenderness of the semi-circular beam is defined as $R / h$ . The first test of this example further validates the accuracy of the proposed formulation. A discretization with $N_{p} = 31$ and $P = 4$ is used. Furthermore, the slenderness is set as $R / h = 10$ , representing relatively thick beams. The comparisons between the numerical results and analytical solutions are presented in Figure 8(a)–(d), respectively, for the displacement component $u_{1}$ , cross-sectional rotation, axial force, and bending moment. Excellent matches between numerical results and analytical solutions are clearly observed for all quantities of interest. The deformed configuration is plotted in Figure 9 with the color scale indicating the total displacement of material particles.

Figure 8.

A semi-circular beam—numerical results with $R / h = 10 .$ (a) Displacement component $u_{1}$ . (b) Cross-sectional rotation. (c) Axial force. (d) Bending moment.

Figure 9.

A semi-circular beam—Deformed configuration with $R / h = 10$ .

Convergence properties of the proposed formulation are further inspected with respect to the relative errors in the displacement of the beam axis, cross-sectional rotation, axial force, and bending moment. Since the analytical solutions of this example exist, they are used as reference ones for the calculation of relative errors. Convergence results with $R / h = 10$ are exhibited in Figure 10. Analogous to the previous example, the utilization of peridynamic functions of higher polynomial degree not only improves numerical precision, but also enhances the order of convergence. Furthermore, the order of convergence for even and odd values of $P$ is still reported at $P$ and $P - 1$ .

Figure 10.

A semi-circular beam—convergence results with $R / h = 10$ . (a) Displacement of the beam axis. (b) Cross-sectional rotation. (c) Axial force. (d) Bending moment.

Figure 11 portrays the deformed configuration with $R / h = 10^{3}$ obtained with $N_{P} = 31$ and $P = 4$ . By juxtaposing Figures 9 and 11, one can easily recognize the discrepancy in the deformed configurations obtained with $R / h = 10$ and $10^{3} .$ Here, it is worth mentioning that the applied forces are appropriately normalized by the membrane and bending rigidities to ensure that the deformed configurations are insensitive with respect to the slenderness. This feature implies that, with the used discretization, the proposed formulation provides accurate results for $R / h = 10$ , but fails to properly predict the response of the beam for $R / h = 10^{3}$ . A closer examination of the deformed configuration in Figure 11 reveals that the beam axis nearly deforms into a circular configuration of a smaller radius, meaning that the axial deformation is more predominant than the bending and shear counterparts. This suggests that, besides the shear locking effect, the proposed formulation also suffers from the membrane locking effect.

Figure 11.

A semi-circular beam—deformed configuration with $R / h = 10^{3}$ .

The performance of the proposed formulation with $R / h = 10^{3}$ is further elucidated by means of convergence tests. As shown in Figure 12, only the convergence properties of the axial force are in good correspondence with those for the slenderness of 10, i.e., Figures 10(c) and 12(c). On the contrary, for other kinematic quantities, the order of convergence obtained with $P = 2, 3, 4,$ and $5$ decreases significantly as the improvement in relative errors is marginal with the increasing number of particles. For $P = 6$ and $7$ , their asymptotic orders of convergence are only obtained with a sufficient number of particles.

Figure 12.

A semi-circular beam—convergence results with $\frac{R}{h} = 10^{3} .$ (a) Displacement of the beam axis. (b) Cross-sectional rotation. (c) Axial force. (d) Bending moment.

A closer inspection reveals the deterioration in the convergence of the axial force with $P = 6$ and $7$ , i.e., Figures 10(c)and 12(c). More precisely, the order of convergence does not approach the asymptotic value of $P$ and $P - 1$ . However, before the deterioration happens the asymptotic values are already achieved with the small number of particles used. The reason for this phenomenon is still not clear. Since this occurs for both the values of slenderness, the locking effects can be excluded. Furthermore, the relative errors are only $10^{- 7}$ , which is still far from machine precision. On the contrary, different schemes have been applied to improve the condition of the system of linear equations, but the deterioration is not yet successfully solved.

4.3. A quintic beam

The two previous examples involve only simple beams having vanishing or constant curvature. In this last example, a quintic beam is considered to represent beams having arbitrary geometry. The geometrical configuration along with loading and boundary conditions are illustrated in Figure 13(a). The beam axis is described as follows:

x = 60 η^{5} - 150 η^{4} + 140 η^{3} - 60 η^{2} + 20 η y = 100 η^{4} - 200 η^{3} + 100 η^{2}

(51)

with $η \in [0, 1]$ being the parameter. A discretization of the beam is visualized in Figure 13(b). Here, it is worth mentioning that the distribution of material particles is not uniform, i.e., the larger the magnitude of curvature, the smaller the spacing between particles.

Figure 13.

A quintic beam. (a) Problem setup. (b) A peridynamic discretization.

For beams having beam axis with varying curvature, the equilibrium conditions contain non-constant coefficients, i.e., see Equations (14)–(16). Hence, it is challenging to derive analytical solutions for displacements and cross-sectional rotation if one prescribes boundary conditions and external loads. In fact, it is more straightforward to assign prior functions to displacements and cross-sectional rotation, and boundary conditions together with external loads are ascertained accordingly. These boundary conditions and external loads are then used for implementation of the proposed formulation, and the specified functions of the displacements and cross-sectional rotation are viewed as analytical solutions. In the coming tests, the displacements and cross-sectional rotation are specified as follows:

u_{1} = u_{2} = ϕ = \sin (\frac{2 π S}{L}) .

(52)

With the above functions, the beam is fixed at both ends, and the external loads are determined by substituting Equation (51) into Equations (14)–(16). Distributions of external loads are shown in Figure 14.

Figure 14.

A quintic beam—external forces and moment with $d / h = 10 .$ (a) Distributed forces. (b) Distributed moment.

The slenderness of the quintic beam is characterized by $d / h$ , with $d$ being the horizontal span of the beam (see Figure 13). Figure 15 exhibits the comparison between the numerical results and analytical solutions for the slenderness $d / h = 10$ . For this comparison, a discretization with $N_{P} = 31$ and $P = 6$ is used. For all inspected kinematic quantities, excellent matches between the numerical results and analytical solutions are recognized, validating the accuracy of the proposed formulation. The deformed configuration is visualized in Figure 16 with the color scale indicating the total displacement of material particles.

Figure 15.

A quintic beam—numerical results with $d / h = 10$ . (a) Displacement component $u_{1}$ . (b) Cross-sectional rotation. (c) Axial force. (d) Bending moment.

Figure 16.

A quintic beam—deformed configuration with $d / h = 10$ .

To close this example, the examination on the sensitivity of the relative errors in the displacement of the beam axis with respect to the slenderness $d / h$ is performed. In this test, four values of slenderness are considered, i.e., $d / h = 10, 10^{2}, 10^{3}$ , and $10^{4}$ , to represent from relatively thick to extremely thin beams. Furthermore, a discretization with $N_{P} = 101$ is used. As plotted in Figure 17, for $P = 2$ and $3$ , the relative errors quickly increase as the slenderness changes from $10$ to $10^{2}$ , then remain flat for the slenderness of $10^{3}$ and $10^{4}$ . Regarding $P = 4$ and $5$ , although the relative errors also increase with higher values of slenderness, the percentage difference is less than $2.5 %$ , which is practically acceptable. On the contrary, a great insensitivity with respect to the slenderness is noticed with $P = 6$ and $7$ .

Figure 17.

A quintic beam—relative error in displacement of the beam axis with different values of slenderness.

As a conclusive remark, membrane and shear locking effects are present in the numerical performance of the proposed formulation, but they can be significantly mitigated by utilizing peridynamic functions of higher polynomial degree.

5. Conclusion

This study presents a peridynamic formulation for the analysis of planar arbitrarily curved beams considering kinematic assumptions of Timoshenko–Ehrenfest beam model. Displacement components of the beam axis and rotation of the cross-section are used to describe linear kinematic relations. Equilibrium conditions and boundary conditions are derived by means of the principle of virtual work. Then, peridynamic functions are incorporated into the primal forms of equilibrium conditions to develop a peridynamic formulation.

Several rigorous tests are designed to assess the numerical performance of the proposed formulation. For relatively thick beams, the proposed formulation provides a good level of accuracy regardless of the highest polynomial degree used for the construction of peridynamic functions. However, the membrane and shear locking effects are found in the analysis of slender beams. In this scenario, it is suggested to utilize peridynamic functions of higher polynomial degree for the analysis of slender beams to ensure numerical precision. Indeed, the membrane and shear locking effects are significantly alleviated with the use of peridynamic functions of higher polynomial degree. Regarding the convergence properties, the order of convergence for even and odd values of $P$ is, respectively, determined as $P$ and $P - 1$ . Here, $P$ represents the highest polynomial degree of peridynamic functions.

Footnotes

Appendix

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This Research is funded by Thailand Science research and Innovation Fund Chulalongkorn University (awarded to J.R.). This research project is supported by the Second Century Fund (C2F), Chulalongkorn University, Thailand (awarded to D.Vo). P.S. acknowledges partial support from Chiang Mai University, Thailand.

ORCID iDs

Duy Vo

Jaroon Rungamornrat

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. Data will be made available on request.

References

Silling

. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000; 48(1): 175–209.

Silling

Epton

Weckner

, et al. Peridynamic states and constitutive modeling. J Elasticity 2007; 88(2): 151–184.

Silling

Lehoucq

. Peridynamic theory of solid mechanics. In: Aref

Van Der Giessen

(eds) Advances in applied mechanics. Amsterdam: Elsevier, 2010, pp. 73–168.

Madenci

Oterkus

. Peridynamic theory. In: Madenci

Oterkus

(eds) Peridynamic theory and its applications. New York: Springer, 2014, pp. 19–43.

Madenci

Roy

Behera

. Advances in peridynamics. Cham: Springer International Publishing, 2022.

Breitenfeld

Geubelle

Weckner

, et al. Non-ordinary state-based peridynamic analysis of stationary crack problems. Comput Methods Appl Mech Eng 2014; 272: 233–250.

Amani

Oterkus

Areias

, et al. A non-ordinary state-based peridynamics formulation for thermoplastic fracture. Int J Impact Eng 2016; 87: 83–94.

Ren

Askari

. A 3D discontinuous Galerkin finite element method with the bond-based peridynamics model for dynamic brittle failure analysis. Int J Impact Eng 2017; 99: 14–25.

Javili

Morasata

Oterkus

, et al. Peridynamics review. Math Mech Solids 2018; 24(11): 3714–3739.

10.

Hattori

Trevelyan

Coombs

. A non-ordinary state-based peridynamics framework for anisotropic materials. Comput Methods Appl Mech Eng 2018; 339: 416–442.

11.

Imachi

Tanaka

Bui

, et al. A computational approach based on ordinary state-based peridynamics with new transition bond for dynamic fracture analysis. Eng Fract Mech 2019; 206: 359–374.

12.

Diehl

Prudhomme

Lévesque

. A review of benchmark experiments for the validation of peridynamics models. J Peridyn Nonlocal Modeling 2019; 1(1): 14–35.

13.

Imachi

Tanaka

Ozdemir

, et al. Dynamic crack arrest analysis by ordinary state-based peridynamics. Int J Fract 2020; 221(2): 155–169.

14.

Jafaraghaei

Bui

. Peridynamics simulation of impact failure in glass plates. Theor Appl Fract Mech 2022; 121: 103424.

15.

Wang

Tanaka

Oterkus

, et al. Study on two-dimensional mixed-mode fatigue crack growth employing ordinary state-based peridynamics. Theor Appl Fract Mech 2023; 124: 103761.

16.

Naumenko

Pander

Würkner

. Damage patterns in float glass plates: experiments and peridynamics analysis. Theor Appl Fract Mech 2022; 118: 103264.

17.

Diyaroglu

Oterkus

, et al. Peridynamics for bending of beams and plates with transverse shear deformation. Int J Solids Struct 2015; 69-70: 152–168.

18.

Nguyen

Oterkus

. Peridynamics formulation for beam structures to predict damage in offshore structures. Ocean Eng 2019; 173: 244–267.

19.

Jafari

Ezzati

Atai

. Static and free vibration analysis of Timoshenko beam based on combined peridynamic-classical theory besides FEM formulation. Comput Struct 2019; 213: 72–81.

20.

Yang

Oterkus

Nguyen

, et al. Implementation of peridynamic beam and plate formulations in finite element framework. Contin Mech Thermodyn 2019; 31(1): 301–315.

21.

Yang

Oterkus

. Peridynamic formulation for Timoshenko beam. Procedia Struct Integr 2020; 28: 464–471.

22.

Nguyen

Oterkus

. Peridynamics for geometrically nonlinear analysis of three-dimensional beam structures. Eng Anal Bound Elem 2021; 126: 68–92.

23.

Shen

Xia

, et al. Modeling of peridynamic beams and shells with transverse shear effect via interpolation method. Comput Methods Appl Mech Eng 2021; 378: 113716.

24.

Bai

Naceur

Liang

, et al. Alleviation of shear locking in the peridynamic Timoshenko beam model using the developed mixed formulation method. Comput Part Mech 2023; 10(3): 627–643.

25.

Bai

Naceur

Zhao

, et al. Improved numerical integration for locking treatment in the peridynamic Timoshenko beam model. Eng Comput 2023; 40(9/10): 2225–2247.

26.

Roy

Behera

Madenci

. Peridynamic modeling of elastic instability and failure in lattice beam structures. Comput Methods Appl Mech Eng 2023; 415: 116210.

27.

Prathap

Bhashyam

. Reduced integration and the shear-flexible beam element. Int J Numer Methods Eng 1982; 18(2): 195–210.

28.

Raveendranath

Singh

Pradhan

. A two-noded locking-free shear flexible curved beam element. Int J Numer Methods Eng 1999; 44(2): 265–280.

29.

Lee

P-G

Sin

H-C

. Locking-free curved beam element based on curvature. Int J Numer Methods Eng 1994; 37(6): 989–1007.

30.

Bouclier

Elguedj

Combescure

. Locking free isogeometric formulations of curved thick beams. Comput Methods Appl Mech Eng 2012; 245–246: 144–162.

31.

Nanakorn

. A total Lagrangian Timoshenko beam formulation for geometrically nonlinear isogeometric analysis of planar curved beams. Acta Mech 2020; 231(7): 2827–2847.

32.

Nanakorn

. Geometrically nonlinear multi-patch isogeometric analysis of planar curved Euler–Bernoulli beams. Comput Methods Appl Mech Eng 2020; 366: 113078.

33.

Nanakorn

Bui

. A total Lagrangian Timoshenko beam formulation for geometrically nonlinear isogeometric analysis of spatial beam structures. Acta Mech 2020; 231(9): 3673–3701.

34.

Nanakorn

Bui

. Geometrically nonlinear multi-patch isogeometric analysis of spatial Euler–Bernoulli beam structures. Comput Methods Appl Mech Eng 2021; 380: 113808.

35.

Nanakorn

Bui

, et al. On invariance of spatial isogeometric Timoshenko–Ehrenfest beam formulations for static analysis. Comput Methods Appl Mech Eng 2022; 394: 114883.

36.

Nanakorn

Bui

, et al. Locking-free isogeometric Timoshenko–Ehrenfest beam formulations for geometrically nonlinear analysis of planar beam structures. Mech Adv Mater Struct 2022; 31(2): 356–377.

37.

Wang

Chen

J-S

. A locking-free meshfree curved beam formulation with the stabilized conforming nodal integration. Comput Mech 2006; 39(1): 83–90.

38.

Madenci

Barut

Futch

. Peridynamic differential operator and its applications. Comput Methods Appl Mech Eng 2016; 304: 408–451.

39.

Madenci

Barut

Dorduncu

. Peridynamic differential operator for numerical analysis. Cham: Springer International Publishing, 2019.

40.

Van Nguyen

Nanakorn

. Effects of discretization schemes on free vibration analysis of planar beam structures using isogeometric Timoshenko-Ehrenfest beam formulations. In: Geng

Qian

Poh

L-H

, et al. (eds) Proceedings of the 17Th East Asian-pacific conference on structural engineering and construction. Singapore: Springer Nature Singapore, 2023, pp. 829–836.

A peridynamic formulation for planar arbitrarily curved beams with Timoshenko–Ehrenfest beam model

Abstract

Keywords

1. Introduction

2. Equilibrium conditions and boundary conditions of planar curved beams with Timoshenko–Ehrenfest beam model

2.1. Linear kinematic relations

2.2. Equilibrium conditions and boundary conditions

3. Peridynamic formulation for planar arbitrarily curved beams

3.1. Discretization and numerical integration scheme

3.2. Peridynamic functions

3.3. Peridynamic formulation for planar curved beams

4. Illustrative examples

4.1. A straight beam subjected to sinusoidal force

4.2. A semi-circular beam

4.3. A quintic beam

5. Conclusion

Footnotes

Appendix

Declaration of conflicting interests

Funding

ORCID iDs

Data availability

References