Numerical analysis of large elasto-plastic deflection of constant curvature beam under follower load

International Journal of Non-Linear Mechanics 84 (2016) 46–55

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Numerical analysis of large elasto-plastic deflection of constant curvature beam under follower load D. Pandit n, Sivakumar M. Srinivasan Department of Applied Mechanics, IIT Madras, Chennai 600036, India

art ic l e i nf o

a b s t r a c t

Article history: Received 6 January 2016 Received in revised form 31 March 2016 Accepted 26 April 2016 Available online 27 April 2016

This paper describes a method to analyze the elasto-plastic large deflection of a curved beam subjected to a tip concentrated follower load. The load is made to act at an arbitrary inclination with the tip tangent. A moment-curvature based constitutive law is obtained from linearly hardening model. The deflection governing equation obtained is highly non-linear owing to both kinematics and material nonlinearity. It is linearized to obtain the incremental differential equation. This in turn is solved using the classical Runge–Kutta 4th order explicit solver, thereby avoiding iterations. Elastic results are validated with published literature and the new results pertaining to elasto-plastic cases are presented in suitable non-dimensional form. The load to end angle response of the structure is studied by varying normalized material and kinematic parameters. It is found that the response curves overlap at small deflection corresponding to elastic deformation and diverge for difference in plastic property. The divergent response curves intersect with each other at higher deflection. The results presented also show that the approach may be used to obtain desired non-uniformly curved beam by suitably loading a uniform curvature beam. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Large deflection Linear hardening Curved beam Incremental formulation Non-linear differential equation Moment-curvature

1. Introduction Large deflection of slender elastic beam under the influence of a terminal force a.k.a elastica has been of interest since the time of Euler [1]. With the invention of smart materials, engineering applications involving large beam deflection have expanded considerably over the last few decades [2]. Many of such applications involve curved beams with terminally acting follower forces inducing inelastic large deflections. In the present work an attempt is made to introduce a simple approach to predict and subsequently study these problems. The elastica problem in its simplest form is that of a horizontal slender cantilever under the action of a conservative vertical terminal load. The large deflection in elastica is essentially a quasistatic phenomenon due to finite rotation under small strain condition applied to an Euler–Bernoulli beam. Assuming this condition to hold true, the elastica problem has evolved to include material non-linearity, non uniform cross section, and non-conservative forces etc. The governing equation of the elastica is a two point boundary value problem and is non-linear due to geometric non-linearity. Analytical solutions to such problems are limited to evaluation of n

Corresponding author. E-mail address: [email protected] (S.M. Srinivasan).

http://dx.doi.org/10.1016/j.ijnonlinmec.2016.04.013 0020-7462/& 2016 Elsevier Ltd. All rights reserved.

elliptic or Jacobian integrals. For a recent review on analytical methods, the reader may find [3] useful. The analytical solutions are implicit in nature for load or displacement and are generally expressed in terms of the end slope. To obtain explicit load–displacement responses required in engineering applications, various semi-analytical techniques are devised. Homotopy perturbation and Adomian decomposition methods are two of the most popular semi-analytical techniques; some of the relevant ones are [4–8] and many of the references therein. The semi analytical techniques generally deal with simpler variations of elastica dealing with elastic problems and render solutions in the form of long expressions. To solve more complicated problems of elastica like non-linear elasticity, various numerical approaches are adopted. Though FEM appears to be the most popularly adopted approach, many complicated elastica problems can be solved by easier or more computationally economic approaches. These non-FEM approaches may be broadly looked into as an algorithm which involve a numerical integration scheme coupled with a root finding iterative technique [9,10]. Yu and Johnson [11] coined the terminology ‘plastica’ indicating an extension of the closed form elastica theory to incorporate plasticity. They solved the problem a cantilever under conservative compressive force using the perturbation technique coupled with numerical integration, considering an elastic-perfectly plastic materials model. Refs. [12,13] considered bi-linear elasto-plastic moment–curvature based constitutive law

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

in solving large deflection problems of slender structures. Non-conservative terminal force which follow the deformed beam profile, induce additional complexity for analysis as it poses the question of adequacy of static approaches. It was Beck [14] who first estimated the buckling load of an elastic column under a tangentially acting compressive follower load by dynamic analysis. Rao and Rao [15] found that static approaches are applicable to sub-tangential follower load systems, provided the load is lesser than the critical flutter (dynamic instability) causing force. Though an extensive literature exist for problems of elastic straight beam under follower forces [16–19], relatively less literature deals with curved beam under follower forces. In the seminal work of Argyris and Symeonidis [20], elastic curved beam under various follower loads is studied in depth by employing FEM. Spric and Saje [21] used finite difference method to study the large deflection of curved elastica under tip concentrated follower load. Nallathambi et al. [22] solved the curved cantilever beam under a tip concentrated follower load by numerical integration starting from free end coupled with one parameter shooting method. Shvartstman [23] solved the same by direct integration of an initial value problem obtained by change of variable from the two point boundary value problem. Additionally he showed that a curved elastic beam under follower load can become unstable only by flutter i.e. to say divergence instability does not exist. Very recently, Ghuku and Saha [24] obtained closed form solution of a curved cantilever elastica under dead load. So far, only planar cases have been discussed. In a recent work by [25], a comprehensive analysis has been carried out on both in-plane and out-ofplane response using an intrinsic formulation and shooting. The motivation of the present work is driven by a relative dearth of literature dealing in large elasto-plastic deflection of curved beam under follower load. Evidently additional complexity in terms of material non-linearity is brought in by considering plasticity into the existing geometrically non-linear problem. The aim of the paper is to present a simple explicit methodology of solving such problems and analyze the numerical results. Such problems when studied in detail may find application in fields like lumbar spine prosthetic, smart structures, flexible robotic arms, and nonuniform curvature hook manufacturing. The methodology adopted here consists of linearizing the governing nonlinear equation about current time to obtain an incremental form of the differential equation for the current step which in turn is solved by employing Runge–Kutta 4th order method. The problem is solved using the static method hence the current “time” actually refers to pseudo time instants. The mathematical formulation for a general elasto-plastic curved beam is presented in Section 2. In Section 3 the solution methodology is explained. In Section 4, results pertaining to various curves and loading are presented and discussed. An example application of the results is also presented herein.

Fig. 1. Schematic of a general curved cantilever under a general follower force at its tip.

2.1. Assumption and scope The beam is assumed to follow Euler–Bernoulli hypothesis and hence axial and shear deformations are neglected. It is assumed to be of uniform curvature in the un-deformed state. It is made up of homogeneous isotropic elastic-linearly hardening rate independent material. A quasi static load acts in the plane of symmetry of the beam following the deformation of the beam maintaining a constant angle (α) with the tip tangent (see Fig. 1). The beam is assumed to undergo large curvature change within small strain framework. In this work, static analysis methodology is adopted, hence the solutions are limited to statically stable configurations. 2.2. Governing relations In Fig. 1, the un-deformed configuration is shown in which a Cartesian coordinate system XOY is chosen wherein OX axis is directed along the right horizontal while OY points vertically down. The arc length is measured from the fixed end along the deformed beam and is denoted by s. The tangent at any point on the beam axis makes an angle ϕ (s ) with OX axis. For convenience in notation, the end angle is denoted by ψ i.e. ϕ (s = smax = l ) = ψ . The load P is assumed to follow the beam such that it maintains α orientation with respect to the beam tangent at its tip. To derive the governing equation of deformation from statics, we note from kinematics, the curvature is given by:

κ (s , t ) = 2. Formulation A constant curvature cantilever of length l and included angle γ completely describes its geometry in un-deformed configuration, as shown in Fig. 1. With tip follower load P oriented at an angle α with respect to the tangent at its tip deforms the beam to produce an end angle of ψ with horizontal. Clearly, in the un-deformed configuration, ψ ¼ γ when P¼ 0. The assumptions, scope and the equation of the curved cantilever are described in detail in this section. The development of the incremental constitutive law used in the present analysis is also presented in this section.

47

∂ϕ ∂s

(1)

Considering kinetics at any s along the deformed beam, and differentiating the hogging bending moment M with respect to s, we get:

∂M = − P sin (α − ψ + ϕ) ∂s

(2)

A general rate independent material model relating the moment and curvature, may be conveniently expressed as:

D=

dM dκ

(3)

In Eq. (3), D denotes the elasto-plastic flexural rigidity of the beam.

48

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

In the linear elastic case it is a constant and is equal to EI. In the elasto-plastic case, for plastic loading, it is the slope of moment– curvature non linear curve and in all other cases of linear elastic loading or unloading, it is equal to EI. In industrial practice of FEM analysis of steel structures, a bi-linear material model is often used for elasto-plastic analysis. The corresponding moment–curvature relationship for a beam of rectangular cross section of width b and thickness 2a is given by: (under monotonic load1)

M⁎ = κ ⁎, κ ⁎ ≤ 1 ⎛ 1⎛ 1 − m⎞ 3⎞ ⎟ + m ⎜ κ ⁎ − ⎟, M⁎ = ⎜ 3 − ⎝ 2⎝ 2⎠ κ ⁎2 ⎠

κ⁎ > 1

Here, initial yield curvature is κy0 =

σy Ea

, initial yield stress in uni-

axial tension is sy, the initial yield moment is My0 = normalized moment is M⁎ =

κ⁎ =

M My0

(4)

2ba2σy 3

, yield

and yield normalized curvature is

κ . κ y0

In Eq. (4), m is a dimensionless parameter referred here as postyield parameter. It is a quantity such that, after yielding of the material, mE is the tangent modulus. Clearly, m can have a value from 0 to 1, where m ¼0 represents elastic-perfectly plastic material and m ¼ 1 represents completely elastic material model. To obtain the governing deformation equation, the kinematic condition of Eq. (1), kinetics Eq. (2) and the constitutive law given in Eq. (3), are combined. Subsequently, comparing the coefficient of ds, the governing non-linear differential equation is obtained as:

D

∂ 2ϕ + P sin (α − ψ + ϕ) = 0, ∂s2

0≤s≤l

(5)

The boundary conditions are given by:

∂ϕ γ |s = l = ∂s l

ϕ|s = 0 = 0,

(6)

2.3. Incremental constitutive law

=

dκ ⁎

−

dκ e⁎

(7)

Since, increment in moment is assumed to depend on increment of elastic curvature only, the moment–curvature relationship, in view of Eq. (7) reads:

dM⁎ = dκ ⁎ − dκp⁎

(8)

Under monotonic loading, when the current yield moment M⁎y evolves, the slope of moment–curvature curve also changes in view of Eq. (4). For a given M⁎y , considering an isotropic hardening rule, this slope remains the same at the start of reverse plastic dM⁎ deformation. In order to determine this slope (i.e. D⁎ = dκ ⁎ ), the κ ⁎ corresponding to the known M⁎y , is required to be found from the following cubic equation: (obtained from Eq. (4))

f

(κ ⁎)

=

M⁎y

⎛ 1⎛ 1 − m⎞ 3⎞ ⎟ − m ⎜ κ⁎ − ⎟ = 0 − ⎜3 − ⎝ 2⎝ 2⎠ κ ⁎2 ⎠

(9)

In the present work, Eq. (9) is solved analytically by employing Vieta's procedure. From the obtained κ ⁎, the slope and its 1

The derivation is provided in Appendix A.1.

derivative may subsequently be computed from:

1−m κ ⁎3 ∂D⁎ −3 (1 − m) = ∂κ ⁎ κ ⁎4

D⁎ = m +

(10)

With the relevant details furnished thus far, the incremental form of the complete moment–curvature relationship along with the hardening rule is presented in Algorithm (1). The law is also pictorially depicted in Fig. 2 obtained by employing a return mapping algorithm as given below: Algorithm 1. Incremental moment–curvature relation (superscript “n” implied). 1:

To obtain an incremental moment–curvature law, it is noted that the current yield moment M⁎y will evolve according to Eq. (4) when plastic deformation takes place. Plastic deformation is quantified by introduction of a quantity called plastic curvature κp⁎, which is the difference between total curvature κ ⁎ and elastic curvature κe⁎, defined in incremental form as follows:

dκp⁎

Fig. 2. Moment-curvature relationship for a beam with rectangular cross section.

f (κ ) = My −

1 2

(3 − ) − m (κ − ) = 0 ! solve for κ, given 1 −m κ2

3 2

My 2:

if M2 − My2 = 0, & Mdκ > 0 then ! plastic loading step

3:

D=m+

4:

( )=

5:

dMy = D|dκ |

6:

dM = sign (M ) dMy

∂D ∂κ

1 −m κ3

−3 (1 −m) κ4

7:

dκp = dκ − dM 8: else! Elastic loading step 9: D ¼1 ∂D 10: =0

( ) ∂κ

dMy ¼ 0 dM = dκ dκp = 0 14: end if

11: 12: 13:

3. Methodology for numerical solution The governing Eq. (5) in its present form is a two point boundary value problem. By adopting a change of variable technique as suggested in [23], the equation is first transformed into an initial value problem. Since the problem involves elasto-plastic non-linearity which is history dependent, an incremental formulation is devised. Thus the equation is linearized about the

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

current time and solved in steps. A variable θ is introduced to transform the boundary value problem into an initial value problem, defined as:

θ (s ) = α − ψ + ϕ (s )

(11)

With subsequent normalization of pertinent quantities, Eq. (5) reads:

D¯ θ ″ + P¯ sin θ = 0,

(12)

0 ≤ s¯ ≤ 1

And the corresponding initial conditions:

θ| s¯= 1 = α, where: D¯ =

(13)

θ ′| s¯= 1 = γ D , EI

P¯ =

Pl2 , EI

θ′ =

∂θ , ∂¯s

θ″ =

at the fixed end reads:

∂ 2θ ∂¯s 2

s

and s¯ = l . Also, the condition

(14)

θ| s¯= 0 = α − ψ

Since the deformation process is assumed to be quasi-static in nature with end angle control loading, the quantities at any time t can be equivalently considered to be at ψ (t ) instead. To devise an incremental formulation, an incremental operator is introduced for small time increment Δt , which reads:

Δ (·) =

∂(·) ∂(·) Δt = ∂t t ∂ψ

Δψ (15)

ψ (t )

For the ease of readability we introduce an additional operator:

(·)′ =

(16)

It can be inferred from Eqs. (15) and (16) that: Δ ((·)′) = (Δ (·))′ and Δ ((·)″) = (Δ (·))″. To obtain the linearized incremental form of the governing equation, and initial conditions, Δ(·) is applied to Eqs. (12) and (13) respectively. The incremental governing equation now reads:

⎛ ∂D¯ ⎞ D¯ (Δθ )″ + ⎜ θ ″⎟ (Δθ )′ + P¯ cos θ Δθ = − sin θ ΔP¯ ⎝ ∂κ¯ ⎠ In which κ¯ =

∂ϕ ∂¯s

(17)

is the curvature normalized with respect to l

(i.e. κ¯ = κl ). Similar to the governing equation, the incremental form of the initial conditions reads:

Δθ ′| s¯= 1 = 0

(18)

Clearly, Eq. (18) is an initial value problem with dependent variable Δθ and independent variable s¯ in which the known conditions are specified at s¯ = 1. It may be noted that the incremental quantity Δθ is defined over 0 ≤ s¯ ≤ 1 for the step from end angle ψ to ψ + Δψ occurring over time instants t to t + Δt . The end angle, in turn is defined over: γ ≤ ψ (t ) ≤ ψf in which ψ (tmax ) = ψf is the final end angle required at the end of all the time steps. Now since ΔP¯ is an unknown2 for a ψ controlled problem, another complementary initial value problem is formed to aid the solution, defined as follows:

⎛ ¯ ⎞ ¯ ″ + ⎜ ∂D θ ″⎟ v′ + P¯ cos θ v = − sin θ , Dv ⎝ ∂κ¯ ⎠

v = v (s¯ ),

0 ≤ s¯ ≤ 1

(19)

With initial conditions:

v| s¯= 1 = 0,

Δθ = vΔP¯

(21)

From Eqs. (6), (11) and (21) the condition to determine ΔP¯ is obtained, as given by:

ΔP¯ = −

Δψ v| s¯= 0

(22)

As it appears in Eq. (19), state variables that are required to ¯ obtain the solution are: D¯ , ∂D , P¯ and θ. Assuming these state ∂κ¯ variables to be known for a given intermediate ψ, the following steps are followed to obtain solution to Eq. (13) till the end of pseudo time steps when the final end angle ψf is achieved: 1. Knowing the state variables for a given ψ, Eq. (19) is solved for v (s¯ ). 2. ΔP¯ is determined from Eq. (22) for a given Δψ and subsequently P¯ is updated. 3. Δθ is computed from Eq. (21) . 4. Increment in ϕ is computed from: Δϕ = Δψ + Δθ , subsequently ϕ (s¯ ) is updated which renders the profile of the beam at ψ + Δψ . ∂Δϕ 5. Increment in curvature is computed from: Δκ¯ = ∂¯s . ∂D¯ ¯ 6. D and at ψ + Δψ , required for the next step are obtained from ∂κ¯ the material module employing return mapping algorithm.3 In the present problem the initial un-deformed profile of the beam is given by: ϕ (s¯ , t = 0) = s¯γ . For evaluation of the state ¯ variables D¯ and ∂D , a constitutive material module as depicted in ∂κ¯

∂(·) ∂s¯

Δθ| s¯= 1 = 0,

49

v′| s¯= 1 = 0

(20)

Evidently, the introduced variable v is related to Δθ by: 2 For load control case, Eq. (17) can be directly solved by Runge–Kutta 4th method.

Fig. 2 is invoked at each material point along the beam for the dM⁎ incremental input Δκ¯ . It can be shown that D¯ = D⁎ = dκ ⁎ . The output of the material module are D⁎,

∂D⁎ , ∂κ ⁎

updated yield bending

moment ( M⁎y ), bending moment (Mn) and the internal variable called plastic curvature κp⁎. A non-dimensional parameter, which depends on both material and geometry, referred here as elastica parameter ζ = κy l naturally appeared while seeking the relationship between same quantities but normalized differently. Moments and curvatures are normalized differently for ease of implementation of the structural and constitutive module coherently. The elastica parameter appears in the following pertinent relationships:

κ¯ = κ ⁎ζ ∂D¯ 1 ∂D⁎ = ζ ∂κ ⁎ ∂κ¯ ¯ = Ml = M⁎ζ M EI

(23)

The role of ζ and the post yield parameter m in the deformation response of a curved beam is discussed in the next section.

4. Results and discussion In this section first the validity of the method is presented. This is followed by results pertaining to elasto-plastic deformation of both straight and curved beams under various loading conditions. The results presented in this section have been verified to have converged within a maximum of 0.1% change for an increase in the number of time steps by at least 20% keeping the space discretization value to a maximum of 0.0025 of the length of the beam. In the current work, the load follows the deformed beam. Hence, for a convenient depiction of response of the structure, load versus relative end angle i.e. ψ − γ relationship is chosen, which may give a fair intuition of both x and y deflections of the 3

Pseudo-code in Appendix A.2.

50

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

Fig. 3. Elastic response of straight beam (γ ¼0) under normal follower load, present method compared with [19].

tip, simultaneously. In this section the elasto-plastic deformation response of both straight and curved beams are presented with emphasis on the functional dependence on the elastica parameter ζ and the post yield parameter m. 4.1. Validation 4.1.1. Elastic validation Using the presented method, the elastic results pertaining to straight beam are validated with [19] and depicted by Shv2007 legend in the figure. The curved beam results are validated with [23] and shown with Shv2013 legend in the figures. In Fig. 3, the end angle and the variation of tip coordinates for a straight beam with respect to a normal follower load are presented. The validation of the response of a semicircular curved beam under normal and tangential follower forces are presented in Figs. 4 and 5 respectively. 4.1.2. Elasto-plastic validation In the absence of available literature on elasto-plastic response

Fig. 5. Elastic response of a semicircular beam ( γ = π ) under compressive tangential load, present method compared with [23].

for curved beams under follower force, a finite difference method (FDM) based approach is developed.4 For comparison of deflected profiles under elasto-plastic deformation as predicted by present method and FDM, it is required that the fixed end moment M⁎ > 1. Setting M⁎ = 1.2 and m = 0.5, ζ = 2 it is found from the present method that P¯ = 2.95 is required to be applied normally to an initially straight beam. Similarly for the semicircular beam, to obtain M⁎ = 1.2 with m = 0.5, ζ = 2, it is found from the present method that a normal follower load of P¯ = 9.13 is required to be applied. When these loads are used in the FDM based approach, the deflected profiles matched very closely with present method prediction. The deflected profiles obtained from FDM and the present method are compared in Fig. 6. 4.2. Straight beam The straight beam results are obtained when γ ¼0 is used for the condition specified in Eq. (13). For a perpendicularly acting π follower load, α = 2 is used in Eq. (11). In Fig. 7, the response of a

beam for α = 2 and ζ ¼1 is shown. In industrial practice of FEM analysis for steel components, it is common to consider the post 1 yield parameter to be 200 of that of Young's modulus. Hence m ¼0.005 is chosen here. Since m ¼ 1 represents elastic condition, an intermediate value of m ¼ 0.25 is also considered for the purpose of comparison. It can be seen that with decreasing m, the beam yields larger end angle. It is found that the beam with π

moderately high m always turns back (beyond

dP¯ dψ

= ∞ point) with

a greater stiffness than in the forward path. But this turning back phenomenon occurs at higher angles for lower m. In the m ¼0.005 dP¯

case, a near flat ( dψ ≈ 0) response is obtained beyond certain critical end angle ( ψ ≈ 0.6), without any tendency of turning back. It is at this critical angle the elasto-plastic deformation starts and thus the curves corresponding to m < 1 begin to separate out from m ¼1 curve. When m is kept lower than 1 and ζ is varied, it is observed that the response curves do not separate out from a single critical angle, as can be seen in Fig. 8. The curves show non-linear behavior with mutual crossings. It may be emphasized here that though the curves intersect with one another, they do not correspond to the exactly same bend profile, as the material and structural Fig. 4. Elastic response of a semicircular beam ( γ = π ) under normal follower load, present method compared with [23].

4

Discussed briefly in Appendix B.

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

Fig. 6. Deflected profiles as predicted by FDM and present method under normal follower load producing elasto-plastic deformation, with m = 0.5, ζ = 2.

51

Fig. 9. Bent profiles of a semicircular curved beam under normalized force P¯ = 10 acting tangentially and perpendicularly.

properties are different. This inference can deduced from the fact that the solution of the initial value problem defined in Eqs. (12) and (13) is unique. 4.3. Curved beam

Fig. 7. Response of straight beam under normal follower load and ζ = 1, γ = 0 .

To obtain curved beam response, a semicircular beam is choπ sen, i.e. γ = π is set. Both tangential (α ¼ 0) and normal (α = 2 ) load cases are studied. The bending stiffness of a curved beam under normal and tangential follower forces are significantly different. This can be observed from the loaded bent profiles of an elastic semicircular cantilever beam under a normalized load of 10 units, as shown in Fig. 9. Similar to what is observed in the straight beam case, for the semicircular beam also the elasto-plastic deformation curves separate out from the elastic response curve at a critical value of the relative end angle ( ψ − γ ). This phenomenon is observed both for perpendicularly and tangentially loaded cases as shown in Figs. 10 and 11 respectively. However depending on the value of m, the curvature nature of the elasto-plastic response curves are seen to

Fig. 8. Response of straight beam under normal follower load and m = 0.25, γ = 0 .

Fig. 10. Response of semi circular beam under normal follower load and ζ = 1, γ = π .

52

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

Fig. 13. Bent profiles of semi-circular beams under normal follower load, m ¼0.25 and P¯ = 10.67.

Fig. 11. Response of semi-circular beam under tangential compressive follower load and ζ = 1, γ = π .

Fig. 14. Response of semi-circular beam under tangential compressive load and m = 0.25, γ = π .

Fig. 12. Response of semi-circular beam under normal follower load and m = 0.25, γ = π . π

vary for α = 2 case which is not seen to influence much for tangential case. The role of ζ is important only for elasto plastic deformation, hence to study its influence, m ¼0.25 is considered. For a normally loaded semi-circular cantilever, the response observed is similar to that of a straight beam response in development of very high

dP¯ dψ

and intersection of curves. It is intuitive and also observed in Fig. 12 that the stiffness increases with increase in ζ for a given m. In Fig. 12 the curves corresponding to ζ ¼0.5 and ζ ¼1 are seen to intersect for a load of 10.67 and relative angle of 0.609 radian. The corresponding bent profiles of the beams are shown in Fig. 13. As evident from the figure, the profiles are considerably different though they have the same end angle and load. The response curves in tangential loading case are seen to have opposite sense as compared to their perpendicularly loaded counterpart as is shown in Fig. 14. A comparison of response curves for a straight and a semi-circular cantilever under elasto-

Fig. 15. Response of a straight and a semi-circular beam under normal follower load and ζ ¼1, with γ¼ 0 for straight and γ = π for semicircular beam.

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

53

obtain an incremental linear governing differential equation. The resulting equation is solved by employing classical Runge–Kutta 4th order method. The method is explicit and does not require any iterations that is needed in shooting based iterative techniques. New results are obtained for elasto-plastic deformation of a curved beam in which a strong influence of both material and geometry on the response of the beam is observed. A brief account of the structural response on normalized quantities m and ζ are discussed. In all of the response curves presented here, it may be broadly said that for relatively small deformation, the response curves initially overlap before diverging away locally. On the other hand, at considerable high deflection, intersection of the response curves are observed. As pointed out earlier that overlap or intersection generally do not correspond to the same bent profile. But here the initial overlap region corresponds to elastic response so the bent profiles are identical. It is also observed that for relatively high m, the perpendicularly loaded response curves have two loads for a given end angle, due to occurrence of

Fig. 16. Unloading response of a semi-circular beam under normal follower load for m = 0.25, ζ = 0.5, γ = π .

plastic loading is presented in Fig. 15. From the figure it may be observed that a semi-circular beam is stiffer than a straight beam. 4.3.1. Unloading response and application Plastic deformations are a unique kind of nonlinear analysis which cannot be modeled as non-linear elastic problems due to the requirement of incorporation of flow rule in plastically deforming structures. Modeling unloading and reverse loading phenomena are difficult using analytical formulation. It is precisely here that the strength of an incremental formulation like the one presented here is better appreciated. The unloading response of a curved beam with (m, ζ ) = (0.25, 0.5) subjected to perpendicular follower load is compared with its elastic counterpart in Fig. 16. It can be seen that for relatively higher peak loads, the unloaded relative end angles were higher than that of the corresponding peak load relative end angles. When the peak load was 10.67, the loading and unloading path were seen to intersect at a load of 7.39. It may be noted that the intersection point corresponds to different bend profiles, as the stiffness are different along the loading and unloading paths. Interestingly, when the beam was made to a peak load of 7.39, the unloaded relative end angle obtained was 0.838 as compared to 1.303 obtained for peak load 10.67. This phenomenon of different unloading paths, reemphasizes the history dependent nature of plastic deformation. It can be appreciated that a constant curvature beam is easier to manufacture than a non-uniform curvature curved beam. The results presented here may aid in design of such a non-uniformly curved beam. For example let a beam with end angle ψ = π + 0.6 = 3.742 radians is required to be manufactured from a semicircular curved beam of length 131.3 mm and cross section 25  0.625 mm made of ferric material modeled with E¼210 GPa and m¼ 0.25 and yield stress σy = 250 MPa . Such a beam can be obtained by applying a perpendicularly acting follower load of magnitude 37.2 N to the tip of the semicircular beam. The value of the load is computed from the norEI malized load of 6 units of Fig. 16 and using the formula: P = 6 2 . l

5. Concluding remarks An incremental static analysis procedure is presented here to analyze a constant curvature beam made of linearly hardening elastoplastic material under tip concentrated follower load. In this method, the governing equation is linearized about a pseudo time step to

dP¯ dψ

= 0.

Further studies by varying the material and kinematic properties (i.e. m, ζ , α , γ ) of the structure may pave the way for more optimized design and application of flexible elasto-plastic curved beams. It is also shown with an example problem, how a nonuniformly curved beam with a desired end angle may be obtained from a semicircular beam. This observation may also lead to further research in the direction of inverse analysis. In the presented approach the material model is independent to the structure and is “called in” to solve the structural problem, and hence the approach is generic enough to include user defined material models. The current method may easily be extended to solve non prismatic and non-uniformly curved beams, under follower loads. But to obtain solution of conservative load problems, the method should be modified to solve two initial value problems instead of one, within a given load step.

Acknowledgment The authors would like to thank Ms. Bhakti Patel for the valuable discussion.

Appendix A. Moment-curvature relationship A.1. Derivation of monotonic moment–curvature law In Fig. A1, 5 the strain and stress profiles for a rectangular beam cross section of thickness 2a and width b under elasto-plastic state of deformation is depicted. The portions bounded by −c ≤ y ≤ c correspond to linear elastic region. Here the stress–strain relation is given by:

σ = E ϵ;

|ϵ| ≤

σy = ϵy E

(A.1)

where E is the Young's modulus and ϵy is the maximum elastic strain of the material under uni-axial tension. In the regions defined by: c ≤ y ≤ a and −c ≥ y ≥ − a , the deformation which is elasto-plastic is governed by the following linear hardening stress– strain law:

σy ⎞ ⎞ ⎛ ⎛ σ = ⎜ σy + Et ⎜ |ϵ| − ⎟ ⎟ sign (ϵ) ⎝ ⎝ E ⎠⎠ where Et is the tangent modulus and is given by: 5

To increase readability.

(A.2)

54

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

normalized bending moment to be M⁎ =

M My0

and subsequently

simplifying, we get:

M⁎ =

⎛ 1⎛ 1 − m⎞ 3⎞ ⎜3 − ⎟ + m ⎜ κ ⁎ − ⎟, ⎝ 2⎝ 2⎠ κ ⁎2 ⎠

κ⁎ ≥ 1

(A.12)

A.2. Return mapping algorithm In order to solve Eq. (19) for the step ψn to ψn þ 1 ( n = 1, 2, 3…); ∂D⁎ and ( ∂κ ⁎ )n + 1 are required at each node. These are supplied by the constitutive module when invoked at each node. The pseudocode which solves the incremental constitutive law presented in Algorithm 1 by the return mapping algorithm is presented as follows:

D⁎ n + 1

Algorithm 2. Return mapping algorithm (superscript “n” implied). 1:

Fig. A1. Schematic of the strain (ϵ) and stress (s) profiles of a beam cross section under elasto-plastic deformation, −a ≤ y ≤ a .

Et = mE;

0≤m≤1

(A.3)

2a2bE ϵ y My0 = 3

(A.5)

ϵy a

(A.6)

For the deformation state shown in Fig. A1, we have:

ϵ y = κc

(A.7)

Now, noting the symmetry of deformation and stress pattern about the neutral axis of the cross section, the total bending moment (M) of the cross section is given by:

∫0

3:

Output: M n + 1, M yn + 1, κpn + 1, Dn + 1,

4: 5:

! trial yield function evaluation

a

yσb dy

(A.8)

Integrating the integrand of Eq. (A.8) by considering Eqs. (A.1)– (A.3), (A.5)–(A.7), we get:

Eϵy 2 mE ϵ y 2 M Eκc3 mEκ 3 = + (a − c2) + (a − c3) − (a − c2) 2b 3 2 3 2

y0

to read: E ϵ3y ⎛ 1 1 ⎞⎞ 1 ⎞ m⎛ ⁎ 1 ⎞ m⎛ 1⎛ M ⎜κ − 2⎟ − ⎜ 1 − 2 ⎟⎟ = 2 ⎜ + ⎜1 − 2⎟ + 3⎝ 2⎝ 2b 2⎝ κ⁎ ⎠ κ⁎ ⎠ κ⁎ ⎠⎠ κ y0 ⎝ 3κ ⁎2

(A.10)

∂D n + 1 ∂κ

( )

ftr = (M n + dκ )2 − M yn2 ! solve for procedure

κ from Eq. (9) for Myn following Vieta's

κ ¼funVieta(Mny)

8:

if ftr ≥ 0 then ! plastic loading step

(

dMy = m +

1 −m κ3

)|dκ| ! increment in yield moment

(M n )

10: 11:

if

12: 13:

else dM = − (dMy + M yn + M n) ! increment in moment

14: 15: 16: 17: 18: 19: 20: 21: 22:

≥ 0 then dM = dMy + M yn − M n ! increment in moment

end if dκp = dκ − dM ! increment in plastic curvature else ! elastic loading step dM = dκ ! increment in moment dMy ¼0 ! no change in yield moment dκp = 0 ! no change in plastic curvature end if M n + 1 = M n + dM ! update moment

M yn + 1 = M yn + dMy ! update yield moment

23:

κpn + 1 = κpn + dκp ! update plastic curvature

24:

if M yn + 1 − M yn ≤ 0 then ! no change in yield curvature

25: 26:

(A.9)

By noting that the normalized curvature is defined to be: κ κ ⁎ = κ and considering Eqs. (A.6) and (A.7), Eq. (A.9) is simplified

∂D n + 1⎞ ⎟ ∂κ ⎠

( )

7:

(A.4)

and

M=2

Input: M n, M yn, κpn, dκ

9:

ϵ = κy

κy0 =

2:

6:

From Euler–Bernoulli linear elastic beam theory, we know that the maximum elastic bending moment (My0), strain (ϵ)-curvature (κ) law and maximum elastic curvature (κy0) are given by:

procedure CONSTITUTIVE MODULE ⎛ n ⎜M , M n, κ n, dκ , M n + 1, M n + 1, κ n + 1, Dn + 1, y p y p ⎝

27: 28: 29:

implies elastic step Dn + 1 = 1 ! updated tangent modulus ∂D n + 1 ∂κ

( )

=0

else ! plastic loading step

Dn + 1 = m + ∂D n + 1 ∂κ

( )

=

1 −m κ3

! updated tangent modulus

−3 (1 −m) κ4

30: end if 31: end procedure

Eq. (A.4) is simplified in view of Eqs. (A.6) and (A.7) to read:

E ϵ3y 3My0 = 2 2b κy0

(A.11)

Now dividing Eq. (A.10) by (A.11) and noting the definition of

Appendix B. Finite difference method (FDM) It is required to solve Eq. (5) along with the boundary

D. Pandit, S.M. Srinivasan / International Journal of Non-Linear Mechanics 84 (2016) 46–55

conditions defined by Eq. (6) by central difference scheme of FDM. The length of the beam is first discretized into n nodes in which the fixed end node reads 1 and the free end node number is n. To any known trial solution or intermediate solution ϕi0 where i = 1, 2, … , n, denote node number we have error function g defined as:

gi = g (ϕi0− 1, ϕi0, ϕi0+ 1, ϕn0 ) ⎛ ϕ 0 − 2ϕ 0 + ϕ 0 ⎞ i i−1⎟ 0 0 = Di ⎜⎜ i + 1 ⎟ + P sin (α − ϕn + ϕi ); δ2 ⎝ ⎠

(B.1)

, …, n − 1; l

where: δ = n −1 , Di = D (κi ) ; κi = The BC at fixed end reads:

i = 2, 3

ϕi + 1 −ϕi − 1 2δ

.

g1 = ϕ1

(B.2)

The BC at the free end employing ghost point formulation and assuming that the free end always remains in elastic state, reads:

gn =

δγ ⎞ 2EI ⎛⎜ 0 ϕ − ϕn0 + ⎟ + P sin α δ2 ⎝ n − 1 l ⎠

(B.3)

The updating or correction term Δϕi to the known solution ϕ0i for i = 1, 2, … , n, is determined by solving the following system of n linear equations:

∂g Δϕ = − g (ϕ0) ∂ϕ0

(B.4)

The process of updating is:

ϕ0 = ϕ0 + Δϕ

(B.5)

This process of solving the linear system of equations followed by updating is continued until convergence is achieved. It may be noted that like in the incremental method, the same constitutive module is invoked here to evaluate D and ∂D at each node. ∂κ

References [1] V.G.A. Goss, The history of the planar elastica: insights into mechanics and scientific method, Sci. Edu. 18 (2009) 1057–1082. [2] G. Wang, M. Shahinpoor, Design, prototyping and computer simulations of a novel large bending actuator made with a shape memory alloy contractile wire, Smart Mater. Struct. 6 (1997) 214. [3] M. Batista, Analytical treatment of equilibrium configurations of cantilever under terminal loads using Jacobi elliptical functions, Int. J. Solids Struct. 51 (2014) 2308–2326.

55

[4] J. Wang, J.-K. Chen, S. Liao, An explicit solution of the large deformation of a cantilever beam under point load at the free tip, J. Comput. Appl. Math. 212 (2008) 320–330. [5] L. Dong, A. Alotaibi, S. Mohiuddine, S. Atluri, Computational methods in engineering: a variety of primal & mixed methods, with global & local interpolations, for well-posed or ill-posed bcs, CMES: Comput. Model. Eng. Sci. 99 (2014) 1–85. [6] T. Elgohary, L. Dong, J. Junkins, S. Atluri, Solution of post-buckling & limit load problems, without inverting the tangent stiffness matrix & without using arclength methods, CMES: Comput. Model. Eng. Sci. 98 (2014) 543–563. [7] A. Banerjee, B. Bhattacharya, A. Mallik, Large deflection of cantilever beams with geometric non-linearity: analytical and numerical approaches, Int. J. Non-Linear Mech. 43 (2008) 366–376. [8] S. Ghosh, D. Roy, Numeric-analytic form of the adomian decomposition method for two-point boundary value problems in nonlinear mechanics, J. Eng. Mech. 133 (2007) 1124–1133. [9] G. Lewis, F. Monasa, Large deflections of cantilever beams of nonlinear materials, Comput. Struct. 14 (1981) 357–360. [10] K. Lee, Large deflections of cantilever beams of non-linear elastic material under a combined loading, Int. J. Non-Linear Mech. 37 (2002) 439–443. [11] T. Yu, W. Johnson, The plastica: the large elastic-plastic deflection of a strut, Int. J. Non-linear Mech. 17 (1982) 195–209. [12] M. Vaz, M.H. Patel, Post-buckling behaviour of slender structures with a bilinear bending moment-curvature relationship, Int. J. Non-Linear Mech. 42 (2007) 470–483. [13] V. Picandet, N. Challamel, S. Hin, Buckling and post-buckling of gradient and nonlocal plasticity columns experiencing softening, Int. J. Solids Struct. 51 (2014) 4052–4067. [14] M. Beck, Die knicklast des einseitig eingespannten, tangential gedrückten stabes, Zeit. Angew. Math. Phys. (ZAMP) 3 (1952) 225–228. [15] B.N. Rao, G.V. Rao, Applicability of the static or dynamic criterion for the stability of a cantilever column under a tip-concentrated subtangential follower force, J. Sound Vib. 120 (1988) 197–200. [16] V.V. Bolotin, Nonconservative Problems of the Theory of Elastic Stability (Moscow, 1961), English translation published by Pergamon Press, New York, 1963. [17] M. Langthjem, Y. Sugiyama, Dynamic stability of columns subjected to follower loads: a survey, J. Sound Vib. 238 (2000) 809–851. [18] I. Elishakoff, Controversy associated with the so-called follower forces: critical overview, Appl. Mech. Rev. 58 (2005) 117–142. [19] B. Shvartsman, Large deflections of a cantilever beam subjected to a follower force, J. Sound Vib. 304 (2007) 969–973. [20] J. Argyris, S. Symeonidis, Nonlinear finite element analysis of elastic systems under nonconservative loading-natural formulation. Part I. Quasistatic problems, Comput. Methods Appl. Mech. Eng. 26 (1981) 75–123. [21] S. Srpčič, M. Saje, Large deformations of thin curved plane beam of constant initial curvature, Int. J. Mech. Sci. 28 (1986) 275–287. [22] A.K. Nallathambi, C.L. Rao, S.M. Srinivasan, Large deflection of constant curvature cantilever beam under follower load, Int. J. Mech. Sci. 52 (2010) 440–445. [23] B. Shvartsman, Analysis of large deflections of a curved cantilever subjected to a tip-concentrated follower force, Int. J. Non-Linear Mech. 50 (2013) 75–80. [24] S. Ghuku, K.N. Saha, A theoretical and experimental study on geometric nonlinearity of initially curved cantilever beams, Eng. Sci. Technol. Int. J. 19 (1) (2016) 135–146. [25] K.N. Karlson, M.J. Leamy, Three-dimensional equilibria of nonlinear precurved beams using an intrinsic formulation and shooting, Int. J. Solids Struct. 50 (2013) 3491–3504.