- Email: [email protected]
Cartesian and piecewise parametric large deflection solutions of tip point loaded Euler–Bernoulli cantilever beams
Author's Accepted Manuscript
Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded Euler-Bernoulli Cantilever Beams Hafez Tari, G.L. Kinzel, D.A. Mendelsohn
www.elsevier.com/locate/ijmecsci
PII: DOI: Reference:
S0020-7403(15)00244-1 http://dx.doi.org/10.1016/j.ijmecsci.2015.06.024 MS3034
To appear in:
International Journal of Mechanical Sciences
Received date: 2 May 2015 Revised date: 7 June 2015 Accepted date: 28 June 2015 Cite this article as: Hafez Tari, G.L. Kinzel, D.A. Mendelsohn, Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded EulerBernoulli Cantilever Beams, International Journal of Mechanical Sciences, http://dx. doi.org/10.1016/j.ijmecsci.2015.06.024 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded Euler-Bernoulli Cantilever Beams Hafez Tari1 , G.L. Kinzel2 , D.A. Mendelsohn3 Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, USA
Abstract This paper focuses on the large deflection solutions of Euler-Bernoulli cantilever beams loaded at the tip with a force and moment. Cartesian large deflection solutions are given for general tip point loading, as well as a new closed form solution for the case of a large end-force only. Additionally, we develop piecewise parametric large deflection solutions, which can efficiently and effectively quantify the inherently nonlinear and complex load, deformation, and energy relationships exhibited by compliant mechanisms, thus, paving the way for the emergence of novel compliant mechanism syntheses. The utilization of the proposed solution methodology can help the precision engineering community in the design of new devices with fixed size and no dependence on the load magnitudes. We coded the solutions as BeamSol, a stand-alone user-friendly black-box solver, which is envisioned to help researchers in engineering analyses and syntheses of beam applications. Keywords: Large Deflection, Euler-Bernoulli Beam, Automatic Taylor Expansion Technique, Piecewise Solution, BeamSol. 1. Introduction Beams have found many applications in engineering, for instance, in atomic force microscopy [1], and design of actuators and energy harvesters in smart materials [2]. In particular, slender beams are a building-block of compliant mechanisms [3]. Such mechanisms transmit or transform motion, load, or energy via the mobility and flexibility of their members. In contrast to traditional or rigid mechanisms, which are made strong through increasing their stiffness, compliant mechanisms are designed to be flexible and strong. In contrast to their rigid counterparts, compliant mechanisms offer part-count reduction; reduced assembly time; ease of manufacturing 1
Email: [email protected], address all correspondence to this author, Tel./Fax: (614) 292-2289/3361 Email: [email protected] 3 Email: [email protected] 2
Preprint submitted to International Journal of Mechanical Sciences
July 6, 2015
processes; increased reliability and precision; wear, weight and maintenance reduction. Despite their advantages, compliant mechanisms exhibit a complex nonlinear behavior, which poses modeling and computational challenges. “Kinetostatics” alleviates this shortcoming by replacing beams in a compliant mechanism with equivalent pseudo-rigid-body models, which consist of rigid segments connected with rotational springs that compensate for the compliance. This approach is advantageous, as the traditional kinematic modeling techniques can be applied to find the geometric dimensions (e.g., segments length) and compliance (e.g., springs stiffness) for prescribed deformation/load/energy constraints. Utilizing kinetostatics, Burns [4], and Burns and Crossley [5] synthesized a four-bar linkage with a flexible coupler link. Winter and Shoups [6] studied the kinetostatic analysis of path-generating flexible-link mechanisms. Extending the classical dimensional synthesis problems, Huang and Roth [7] formulated the kinetostatic synthesis equations for three types of planar linkages with linear springs, and obtained the solutions by resorting to polynomial solvers. Krovi et al. [8] studied the kinetostatic synthesis of planar coupled serial chain mechanisms using a hybrid approach, which combines precision point synthesis and optimization. Dado [9] developed a variable parametric pseudo-rigid-body model for limit position synthesis of compliant four-bar linkage with energy requirements. These attempted to find a unique or a partial set of pseudo-rigid-body model solutions. However, Tari and Su [10] recently proposed a systematic framework that transforms the modeling equations into a system of polynomials, which can be solved for the entire solution set. A parameter continuation was devised to avoid computation of junk solutions for subsequent syntheses. Later, Tari et al. [11] showed that this framework is effective even for the forward analysis of a 3D compliant Stewart-Gough platform. Extensive work has been done on linear beam theory to obtain explicit solutions for small deformation of cantilever beams subjected to a tip point loading. Such solutions are still the premise of certain devices in precision engineering. However, the overall dimensions of such devices must be adjusted with the magnitudes of the applied loads, as linear theory solutions are not valid once the beam’s axis slope is significant. Despite the success of pseudo-rigid-body models, they are valid for a limited range of motion. When the compliant members undergo large deformations, large deflection solutions must be used for accurate modeling. Finite-element solution procedures can overcome these shortcomings. For instance, building on independent approximations of three-dimensional fields of displacements and stresses, Argyris and Kaˇcianauskas [12] presented semi-analytical finite-elements leading to the development of a higher-order theory for beams. Mazars et al. [13] employed higher order interpolation functions to investigate solutions for an enhanced multifiber beam element accounting for shear and torsion. Their constitutive 2
model avoids shear locking phenomena, and also captures nonlinear cross section warping. Jun et al. [14] introduced a dynamic finite-element method for free vibration analysis of generally laminated composite beams. Their model is based on first-order shear deformation theory and accounts for Poisson effect, couplings between extensional, bending and torsional deformations, shear deformation and rotary inertia. Vo and Thai [15] analyzed vibration and buckling of composite beams with arbitrary lay-ups using refined shear deformation theory that accounts for the parabolical variation of shear strains through the depth of the beam. Utilizing classical Timoshenko beam theory and refined zigzag kinematics, O˜ nate et al. [16] presented a simple linear two-noded beam element sufficient for analysis of composite laminated and sandwich beams. Sadowski and Rotter [17] studied the use of solid continuum finite-elements and shell finiteelements in the modeling of nonlinear plastic buckling behavior of cylindrical metal tubes and shells under global bending. Miao et al. [18] enhanced the reverberation-ray matrix approach to analyze the transient response of a laminated composite beam subjected to an impulse force load. Zhang et al. [19] developed an efficient analytic framework to study transverse vibrations of two parallel Timoshenko beams that are connected by discrete springs and coupled with various discontinuities. Very recently, Lim et al. [20] employed the nonlocal elasticity stress theory and investigated the torsional static and dynamic effects of circular nanostructures subjected to concentrated and distributed torques. The foregoing finite-element solution methodologies are robust for arbitrary beam applications. However, efficient solution methods have been desired for special loading and beam cross sectional geometries and properties. For example, Bisshop and Drucker [21] reported exact implicit elliptic integral solutions of an Euler-Bernoulli cantilever beam of uniform cross section and loaded with a tip point vertical force. Howell [22] generalized these solutions for the case of general combined tip point loading. Nonetheless, the implicitness of the foregoing elliptic integral solutions poses numerical evaluation challenges. To alleviate this drawback, Banerjee et al. [23] employed nonlinear shooting and Adomian’s decomposition methods, and reported numerical and approximate large deflection solutions for concentrated intermediate and endpoint loadings. Chen [24] developed a numerical integral approach for complex loading and variable beam’s cross sectional properties. Employing an Automatic Taylor Expansion Technique (ATET), Tari [25] recently developed parametric large deflection solutions for Euler-Bernoulli cantilever beams loaded at the tip with force and moment. These solutions are valid for the entire beam length, and capture both angular and axial inflection points. We refer to these solutions as single-piece parametric solutions, as they are expanded about a single start point for the whole computational domain. Even though this treatment was shown to be effective, rather 3
high approximation orders are needed for large loading inputs. Unfortunately, the higher the approximation orders, the more the solution evaluation time, and more importantly the more the number of extraneous solutions. In this paper, we develop cartesian large deflection solutions, which are independent of beam’s arc length for applied tip force and moment. Also, we develop explicit closed form solutions for the case of a large pure end-force. In addition, we develop piecewise parametric large deflection solutions, which outperform the efficiency and robustness of the single-piece parametric solutions. The proposed large deflection solutions can help in the design of new devices with fixed size and no dependence on the load magnitudes. Our solutions are coded as a stand-alone user-friendly black-box solver, BeamSol, which is envisioned to facilitate analysis and synthesis of of novel compliant mechanisms. The rest of the paper is organized as follows. Next, we review Euler-Bernoulli beams and their mathematical modeling. Section 3 gives cartesian large deflection solutions, and Section 4 gives the closed form solutions for a large pure end-force. Section 5 presents the piecewise parametric solutions together with BeamSol, and is followed by results and discussions. Finally, summary and conclusions are given. 2. A Brief Review on Euler-Bernoulli Beams A schematic view of a cantilever beam subjected to a tip moment M1 and force f with the inclination φ is depicted in Figure 1, in which E and I are the Young’s modulus and the moment of inertia of the beam’s cross section; s and θ are the arc length and the angular deflection of the beam; and x and y are the horizontal and vertical coordinates of the deflected beam, respectively.
x1
M1
Inflection(κ=0)
x
f
θ
s
φ
θ1
y1
y l, E, I
Undeflected
Figure 1: A schematic view of a tip point loaded cantilever beam.
A common beam type is an Euler-Bernoulli beam, whose large deflection is governed by
κ(x) =
M1 + f sin(φ)(x1 − x) − f cos(φ)(y1 − y) M d2 y/dx2 = , = 2 3/2 EI EI [1 + (dy/dx) ]
4
(1)
where κ is the curvature, and M is the bending moment. The boundary conditions for the cantilever beam are dy (0) = 0. dx
y(0) = 0,
(2)
The behavior of the Euler-Bernoulli beam depends on the choice of the material (E), the geometry of the cross section (I), and even the location of the applied load (M ). It is a common practice to nondimensionalize the parameters as s/l → s, x/l → x, and y/l → y; and finally to parameterize the problem with respect to the beam’s arc length. For the case of a beam with constant cross sectional properties, i.e., constant E and I, the resulting equations are (see, for example, Tari [25])
d2 θ(s) = α2 sin(θ(s) − φ), ds2 dθ(1) = β, θ(0) = 0, ds dx(s) = cos(θ(s)), x(0) = 0, ds dy(s) = sin(θ(s)), y(0) = 0, ds where α2 =
f l2 EI
M1 l EI .
and β =
(3) (4) (5) (6)
The exact solutions to this weakly decoupled boundary value
problem can be written in terms of the solutions of the following implicit elliptic integral relations s=± x=±
θ(s) 0 θ(s) 0
y=±
0
θ(s)
dθ β 2 + 2α2 (cos(θ1 − φ) − cos(θ − φ)) cos(θ)dθ β2
+ 2α2 (cos(θ
1
− φ) − cos(θ − φ))
,
(7)
,
(8)
,
(9)
sin(θ)dθ β 2 + 2α2 (cos(θ1 − φ) − cos(θ − φ))
where the plus and minus signs in Eqs. (7)-(9) are inherent to the beam’s multiple behavior, and often in the literature are denoted as elbow up and elbow down solutions. The elliptic integral relations pose some challenges. First, there is an implicit dependence on the unknown tip rotation θ1 , which complicates the integration procedure, even when solved numerically. Second, depending on input loading parameters, the solutions may be very close to singularities, which drastically complicates the numerical integration. Recall that the signs in the elliptic integral relationships change as the input loading creates inflection points, where the curvature changes sign, such as the one schematically shown in Figure 1. As the number of inflection
5
points increases, so does the number of solution branches. Therefore, obtaining solutions that do not share such difficulties is desirable. 3. Cartesian Large Deflection Solutions The parametric large deflection solutions (8) and (9) are in terms of the intermediate parameter s, which is not explicitly at hand. The goal of this section is to eliminate s from the solutions, and obtain a cartesian large deflection solution, i.e., the vertical large deflection component directly in terms of the horizontal large deflection component or vice versa. To this end, we employ the coordinates transformation ⎡ ⎣
⎤
R(φ)
⎡
cos(φ)
X
⎤ ⎡ sin(φ)
⎤ x
⎦ ⎣ ⎦, ⎦= 1 ⎣ l − sin(φ) cos(φ) y Y
(10)
where R(φ) is the orthogonal rotation matrix, which rotates the components x and y clockwise by the angle φ about the origin. Keeping in mind that a rotation does not change the curvature, Eq. (1) simplifies to d2 Y /dX 2 = α2 Y + α1 , [1 + (dY /dX)2 ]3/2
(11)
where α1 = β−α2 (cos(φ)y1 −sin(φ)x1 )/l, α and β are as before, and the new boundary conditions become
Y (0) = 0,
dY (0) + tan(φ) = 0. dX
(12)
We can integrate Eq. (11) once, and employ the second boundary condition of Eq. (12) to write
1 1 + (dY /dX)2
= | cos(φ)| − α1 Y −
α2 2 Y G(Y ). 2
(13)
Then, integrating the foregoing equation and applying the first boundary condition of Eq. (12), we write X=
Y (X)
0
G(Y ) dY. 1 − G(Y )2
(14)
Equation (14) is not amenable to an explicit closed form solution in terms of elementary functions. However, letting G1 = G(Y (1)) and simplifying the results, a standard form of elliptic integral solution may be obtained as 6
X=±
G(Y )
| cos(φ)|
1 G √ dG. 2 2 1−G β − 2α2 (G − G1 )
(15)
Had we replaced R(φ) with its transpose (RT (φ)) in the coordinates transformation (10), a similar result would have been obtained. However, both solutions are qualitatively similar, and have no major advantage over each other. It is worth mentioning that Ref. [26] also reported a solution procedure for obtaining a cartesian large deflection solution. However, their coordinates transformation is not based on rotation and is somewhat more complex than Eq. (15). The inextensibility of the beam imposes a side condition, which would further constrain the trajectory of the beam’s endpoints as 1=±
G1 | cos(φ)|
√
1 1 dG. 2 2 1−G β − 2α2 (G − G1 )
(16)
It is clear that the new cartesian elliptic integral large deflection solution (15) retains the main shortcomings of the elliptic large deflection solutions (8) and (9). Even though the new solution involves only x and y, it is implicit and highly nonlinear, thus, still requires numerical evaluation techniques. Likewise, the new solution still has two branches (plus and minus solutions), and the identification of the proper branch for any input loading parameters requires careful assessment. The square root in Eq. (13) is bounded below by one, thus, we have 0 < G(Y ) ≤ 1. This paves the way for fairly accurate approximations through neglecting higher powers of the quadratic function G(Y ) in the aforementioned formulae. Such assumptions facilitate obtaining efficient approximate closed form solutions in terms of elementary functions. Nonetheless, we aim to obtain solutions, which are accurate for arbitrary applied loadings. 4. Closed Form Approximate Large Deflection Solutions for Large Pure End-forces Explicit exact large deflection solutions exist for only a pure end-moment loading (zero force, α = 0); see, for example, Ref. [25]. However, no explicit solution is reported in the literature for a general pure end-force loading (zero moment, β = 0). Nonetheless, for a sufficiently large applied end-force with no moment, the results of Ref. [25] suggest that beam’s end angle approaches the force angle, i.e., θ1 ≈ φ. Utilizing this result, Eq. (7) can be simplified and solved analytically. For β = 0 and θ1 = φ, we give the exact solution of Eqs. (3) and (4) as φ θ(s) = φ − 4 arctan(tan( ) exp(−αs)), 4
(17)
which is valid for both elbow up and down solution branches. To obtain the corresponding axial deflection components, we integrate Eqs. (8) and (9), and used the trigonometric identities 7
⎧ ⎨ cos(4 arctan(z)) = 1 − 8 z 22 2 , (1+z ) ⎩
sin(4 arctan(z)) =
2
1−z 4z (1+z 2 )2 ,
to write
RT (φ)
⎤ ⎡ cos(φ) − sin(φ) ⎢ s − x(s) ⎦⎣ ⎣ ⎦ =⎣ sin(φ) cos(φ) y(s) ⎡
⎤
⎡
4 α 4 α
1
1+tan2 (φ/4) exp(−2αs) tan(φ/4) exp(−αs) 1+tan2 (φ/4) exp(−2αs)
−
1 1+tan2 (φ/4)
−
tan(φ/4) 1+tan2 (φ/4)
⎤ ⎥ ⎦,
(18)
in which the manifestation of the rotation matrix R(φ) is in analogy with the cartesian solutions given in the previous section. (a)
4
(b)
α = 24
1.0
0.8
3
0.6
y(s)
θ(s), rad
Numerical Sol. Explicit Sol.
−π/10
−2
Numerical Sol. Explicit Sol.
φ=π
2
0.4 1
0.2 +2
0 0.0
α=6 0.2
0.4
0.6
φ = π/10 0.8
φ=π
0.0 1.0
-1.0
-0.5
s
α = 24
α=6 0.0
0.5
+π/10
φ = π/10 1.0
x(s)
Figure 2: Comparison of the numerical solution and the explicit solutions (17) and (18) for the beam’s (a) angular and (b) axial deflections.
Figure 2 compares the explicit solutions, given by Eqs. (17) and (18), with numerical solutions. For the latter, the numerical differential solver of Mathematica is utilized to solve Eqs. (3) and (4). Then, the obtained numerical angular deflections are substituted into Eqs. (5) and (6), which are integrated numerically to determine the axial deflection components. Therefore, the numerical solution is convergent for only the two smallest load magnitudes (α = 6 and α = 8), and diverges at the vicinity of the clamped point as the load magnitude increases. In contrast, the explicit solutions are convergent for the entire beam length. Because the tip angle converges to the force angle for every case studied (which was the derivation assumption for the explicit solutions), the obtained convergent solutions are, in fact, accurate. The divergent numerical solutions are not shown when plotting the axial solutions in Subfigure 2(b). Note that there is an outrageous discrepancy in the form of the exact circular large deflection solution, given in Ref. [25], and the newly developed closed form approximate large deflection
8
solution. This underlines the complexity of the beam behavior and the need for a robust solution procedure. 5. Piecewise Parametric Large Deflection Solutions Letting ψ(s) θ(s) − φ, we can rewrite Eqs. (3)-(6) as
d2 ψ(s) = α2 sin(ψ(s)), ds2 dψ(1) = β, ψ(0) = −φ, ds dx(s) = cos(ψ(s) + φ), x(0) = 0, ds dy(s) = sin(ψ(s) + φ), y(0) = 0, ds
(19) (20) (21) (22)
which are numerically well-conditioned. For i = 1, . . . , n, let [si , si+1 ] be a subdivision of the computational domain s ∈ [0, 1] such that ∪ni=1 [si , si+1 ] = [0, 1], where s1 = 0 and sn+1 = 1. In addition, let θi (s), xi (s), and yi (s) denote the deflection components corresponding to s ∈ [si , si+1 ]. Employing the automatic Taylor expansion technique [25], the parametric large deflection solutions corresponding to each subdivision are written as expansions about the point s0i as 2
θi (s) ≈ φ + c0i + c1i (s − s0i ) + α
M 1i j=2
xi (s) ≈ x0i +
M 2i j=1
and yi (s) ≈ y0i +
M 3i j=1
1 dj−2 (sin(ψ(s))) j! dsj−2
1 dj−1 (cos(ψ(s) + φ)) j! dsj−1 1 dj−1 (sin(ψ(s) + φ)) j! dsj−1
where s0i ∈ [si , si+1 ], c0i = ψi (s0i ), c1i =
dψi (s0i ) , ds
(s − s0i )j ,
(23)
s=s0i
(s − s0i )j ,
(24)
(s − s0i )j ,
(25)
s=s0i
s=s0i
x0i = xi (s0i ), y0i = yi (s0i ), and M1i , M2i ,
and M3i are the approximation orders. Note that all of the summation terms are at hand, as the higher order derivatives of ψ(s) can be found in terms of c0i and c1i upon the successive differentiation of Eq. (19). Accordingly, there are 4n unknowns, i.e., C4n = {c0i , c1i , x0i , y0i }, which should be identified to fully define the piecewise parametric large deflection solutions. To this end, requiring the continuity in the zeroth and first order derivatives of the piecewise angular solutions at the interfaces, and continuity in piecewise axial deflections, together with the original boundary conditions (20)-(22) constitute the 4n algebraic-trigonometric constraints
9
B(C4n ) =
⎧ ⎪ ⎪ θ1 (0) = 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n ⎨ dθ ds (1) = β
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬
⎪ ⎪ ⎪ x1 (0) = 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ y1 (0) = 0
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
∪
⎧ ⎪ ⎪ θj (sj+1 ) = θj+1 (sj+1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dθj+1 j ⎪ ⎨ dθ ds (sj+1 ) = ds (sj+1 ) ⎪ ⎪ ⎪ xj (sj+1 ) = xj+1 (sj+1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ yj (sj+1 ) = yj+1 (sj+1 )
,
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ j = 1, . . . , n − 1 , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
(26)
which may be solved for the entire unknowns. For instance, taking n = 2 (i.e., i = 1, 2) and M1i = M2i = M3i = 5, the fifth order piecewise ATET solutions for θi (s), xi (s), and yi (s) may be written as ⎧ ⎧ ⎨ φ + c01 + c11 (s − s01 ) + αs1 (s − s01 )2 + c11 αc1 (s − s01 )3 ⎪ ⎪ 1! 2! 3! ⎪ ⎪ , 0 ≤ s ≤ s2 , ⎪ c11 (αc1 η21 −3α4 ) ⎪ αs1 η11 ⎩ 4 5 ⎪ + 4! (s − s01 ) + (s − s ) ⎨ 01 5! θ(s) = ⎧ ⎪ ⎪ ⎨ φ + c02 + c12 (s − s02 ) + αs2 (s − s02 )2 + c12 αc2 (s − s02 )3 ⎪ ⎪ 1! 2! 3! ⎪ ⎪ , s2 ≤ s ≤ 1, ⎪ c12 (αc2 η22 −3α4 ) ⎩ ⎩ αs2 η12 4 5 + 4! (s − s02 ) + (s − s ) 02 5! x(s) = ⎧⎧ c01 ) c01 ) ⎪ ⎪ (s − s01 ) − c11 sin(˜ (s − s01 )2 x01 + cos(˜ ⎪⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ 2 ⎪ ⎪ − c11 cos(˜c01 )+αs1 sin(˜c01 ) (s − s01 )3 − c11 (3αs1 cos(˜c01 )+η11 sin(˜c01 )) (s − s01 )4 ⎪ ⎪ 3! 4! ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ α η sin(˜ c )+η cos(˜ c ) ⎩ s1 31 01 41 01 ⎨ − (s − s01 )5 , 5! ⎧ c02 ) c02 ) ⎪ ⎪ ⎪ ⎪ (s − s02 ) − c12 sin(˜ (s − s02 )2 x02 + cos(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ c2 cos(˜ c02 )+αs2 sin(˜ c02 ) ⎪ ⎪ − 12 (s − s02 )3 − c12 (3αs2 cos(˜c024!)+η12 sin(˜c02 )) (s − s02 )4 ⎪ 3! ⎪ ⎪ ⎪ ⎪ ⎪ ⎩⎪ ⎩ − αs2 η32 sin(˜c025!)+η42 cos(˜c02 ) (s − s02 )5 ,
(27)
0 ≤ s ≤ s2 ,
s2 ≤ s ≤ 1, (28)
and
y(s) = ⎧⎧ c01 ) c01 ) ⎪ ⎪ ⎪ (s − s01 ) + c11 cos(˜ (s − s01 )2 y01 + sin(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ ⎪ c2 sin(˜ c01 )−αs1 cos(˜ c01 ) ⎪ (s − s01 )3 − c11 (3αs1 sin(˜c014!)−η11 cos(˜c01 )) (s − s01 )4 ⎪ − 11 ⎪ 3! ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨⎩ + αs1 η31 cos(˜c015!)−η41 sin(˜c01 ) (s − s01 )5 , ⎧ c02 ) c02 ) ⎪ ⎪ ⎪ (s − s02 ) + c12 cos(˜ (s − s02 )2 ⎪ y02 + sin(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ c2 sin(˜ c02 )−αs2 cos(˜ c02 ) ⎪ ⎪ − 12 (s − s02 )3 − c12 (3αs2 sin(˜c024!)−η12 cos(˜c02 )) (s − s02 )4 ⎪ 3! ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎩ ⎩ + αs2 η32 cos(˜c025!)−η42 sin(˜c02 ) (s − s02 )5 ,
0 ≤ s ≤ s2 ,
s2 ≤ s ≤ 1, (29)
where c˜0i = c0i + φ, αsi = α2 sin(c0i ), αci = α2 cos(c0i ), η1i = αci − c21i , η2i = 4αci − c21i , η3i = αci − 7c21i , and finally η4i = c21i η2i + 3α2si for i = 1, 2. Note that s01 , s02 , and s2 are to be chosen 10
as input parameters, thus, leaving only eight unknowns C8 = {c01 , c11 , c02 , c12 , x01 , x02 , y01 , y02 }, which would be identified upon imposing the boundary conditions B(C8 ), an adaptation of B(C4n ) for this special case. 5.1. Choice of the Number of Subdivisions n The more the number of subdivisions, the more accurate the piecewise solutions. However, more subdivisions introduce a significant number of intermediate unknowns, which would burden the computation. Nevertheless, an increase in the number of subdivisions reduces the need for large approximation orders, which are generally required for attaining a significant accuracy. To get a good balance between solution accuracy and computational cost, we set n = 2, for which the resulting large deflection solutions are highly accurate and compact for low approximation orders. 5.2. Choice of the Segments’ Endpoints si and si+1 Recall that the premise of the piecewise solutions is the Taylor series expansion, which approximates the unknown exact large deflection solutions with polynomials in the arc length of the beam. It is clear that the smaller the computational domain, i.e., the length of the arc length for each piecewise solution, the better the accuracy of the approximations. However, for a fixed number of subdivisions and fixed approximations orders, improving the accuracy of each piecewise solution requires shortening the length of that segment, which would force compromising the accuracy of the remaining segment solutions. Therefore, to avoid this, we subdivide the beam length equally, i.e., [si , si+1 ] = [(i − 1), i]/n for i = 1, . . . , n. 5.3. Choice of the Start Points s0i The piecewise solutions also depend on start points s0i , which must be chosen. Like the parametric solution presented in Ref. [25], start points can simplify the solution procedure. For example, the unknowns {c01 , c1n , x01 , y01 } would become known parameters, if s01 and s0n are taken as s1 = 0 and sn+1 = 1, respectively. In contrast, setting s0i to the center of each segment provides the most accuracy and convergence. 5.4. Independence of the Approximation Orders Any of the approximation orders M1i , M2i , and M3i can be chosen independently. This is significant in two ways. The accuracy of any of the angular and axial deflection components can be exclusively adjusted without disturbing the accuracy of the remaining solution components. Likewise, the accuracy of each segment solution can be adapted without affecting the accuracy
11
of the remaining segment solutions. As an example for the piecewise solutions given by Eq. (27), the accuracy of the angular solution on the domain of s ∈ [0, s2 ] can be controlled independently of that of s ∈ [s2 , 1], and vice versa. 5.5. Solution of the Boundary Constraints B(C4n ) The system of boundary constraints B(C4n ) (26) is a decoupled system of algebraic-trigonometric equations. To solve it efficiently, the angular constraints can be solved for c0i and c1i independently of the axial constraints. Then, the resulting solutions can be inserted into the horizontal and vertical deflection constraints so as to solve for x0i and y0i , respectively. One major advantage of the piecewise solution procedure over the single-piece solution approach of Ref. [25] is that the boundary constraints, i.e., B(C4n ), of the former are far better conditioned than those of the latter. For the single-piece solution, the boundary constraints can be converted into a system of polynomial equations, and solved for the entire solution set, among which the true solutions can be sought. However, this approach is inefficient, and impractical for daily engineering applications. In contrast, we employ the Newton-Raphson technique, which initiates at a given start solution and iteratively updates it to converge to a feasible solution. Even though multiple mathematical solutions may exist, we have found that the Newton-Raphson technique, starting at a zero solution set, has always converged to the true physical solution for B(C4n ). This feature has led us to coding the algorithm as a black-box simulation software, which is discussed in the next section. 6. Results and Discussions For the simulations given in the following, we take advantage of the symmetry of the loading, and assume that the force angle is bounded in the first and second quadrants, i.e., φ ∈ [0, π]. Also, with the piecewise solutions, we mean the two-piece solutions. Additionally, for the sake of simplicity, we plot the nondimensional deflection solutions. Finally, for validation purposes, numerical solutions of the exact governing equations were obtained from the numerical differential solver of Mathematica. 6.1. Effect of the Loading Parameters In general, the larger the load magnitudes, the more the deflection revealed by the beam. Figure 3 illustrates that the compact fifth order piecewise solutions in Eqs. (27)-(29) are accurate for a wide range of the input parameters, namely, α, β ∈ [0, 3] for arbitrary φ. In particular, the
12
axial deflections corresponding to the case of α = β = 3 and φ = π, shown in Subfigure 3(b), are rather large which are captured well with the fifth order of approximation. (a)
4
(b)
0.8
Numerical Sol. Piecewise Sol.
φ=π β=3
0.6
0 /1 +π
/10 +3
y(s)
α=3
/10 +3
θ(s), rad
3
2
+3 /10
0.4
φ=π
β=3
+3/10
α=3
0 0.0
α = 3/10 β = 3/10 φ = π/10 0.2
0.4
0.6
0.8
+π/10
Numerical Sol. Piecewise Sol.
0.2
1
α = 3/10
0.0
1.0
-0.4
-0.2
0.0
0.2
s
0.4
β = 3/10
0.6
φ = π/10 0.8
1.0
x(s)
Figure 3: Comparison of the numerical solution and the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4 for the beam’s (a) angular and (b) cartesian deflections.
The magnitude of the force angle also plays a crucial role in the beam’s deflection. For the case of pure force loading, it is intuitive that the more obtuse the force angle, the more the obtained deflection. As shown schematically in Subfigure 4(a), when the force angle is constant, the tip point traces an elliptic-like curve as the force magnitude is increased. In fact, as shown in Subfigure 4(b), the beam’s tip point locus is almost circular for sufficiently small force magnitudes, but diverges to an elliptic curve for large force magnitudes. It can be seen that with the absence of the external moment (β = 0), the fifth order piecewise solutions are valid even for a wide range, i.e., α ∈ [0, 4].
(a)
(b)
1.0
1.0
Inc
se
0.8
in
0.6
α
0.6
0.4
0.4
0.2
0.2
0.0
φ = 5π/6
φ = 2π/3
α=0 -0.5
0.0
x1
0.5
φ = π/2
φ = π/3 φ = π/6
φ=π
y1
y1
0.8
rea
Numerical Sol. Piecewise Sol. β=0
0.0 1.0
-0.5
0.0
x1
0.5
1.0
Figure 4: Beam’s endpoint trajectory: (a) for a typical pure end-force, and (b) for α ∈ [0, 4] with several load angles obtained using the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4.
6.2. The Post-Buckling Behavior For a pure compressive force loading (β = 0 and φ = π), it is quite well known that the beam starts to buckle at the critical force magnitude α = π/2. However, linear beam theory fails 13
to predict beam post-buckling behavior. Shown in Figure 5, the fifth order piecewise solutions are intelligent in recognizing the critical force magnitude, and also give accurate large deflection solutions past this point. (a)
3.0 2.5
(b)
Numerical Sol. Piecewise Sol.
β=0 φ=π
0.8
α = π/2+1.5 0.6
1.5
y(s)
θ(s), rad
2.0
α = π/2+1.5
0.4
1.0
α ε [0,π/2] 0.2
0.4
0.6
0.8
Numerical Sol. Piecewise Sol.
β=0 φ=π
+0.1
.1 +0
α = π/2+0.1
0.5 0.0 0.0
0.2
α = π/2+0.1
0.0 1.0
-0.4
-0.2
0.0
s
0.2
0.4
0.6
α ε [0,π/2] 0.8
1.0
x(s)
Figure 5: Comparison of the numerical solution and the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4 for the beam’s post-buckling (a) angular and (b) cartesian deflections.
6.3. Efficiency of the Piecewise Solutions Figure 6 compares the proposed piecewise solution procedure against the single-piece solution procedure of Ref. [25]. As demonstrated, to get accurate angular large deflection solutions, the piecewise solutions require only sixth order of approximation, while the single-piece solutions require 16th order of approximation for the chosen loading parameters. Even though cartesian large deflection solutions require, in general, higher approximation orders than the angular solutions, this increase is minimal for the piecewise solutions. The axial single-piece solutions require 24th order of approximation, while the cartesian piecewise solutions require only eighth order of approximation. Obviously, the approximation order of the piecewise solutions can be further reduced, if the beam is subdivided into more segments. 6.4. Robustness of the Piecewise Solutions The beam’s angular derivative may be obtained by integrating Eq. (3) which after utilizing Eq. (4) can be written as
dθ(s) ds
2
= β 2 + 2α2 [cos(θ1 − φ) − cos(θ(s) − φ)] Γ(α, β, φ, θ1 , θ),
(30)
where Γ(α, β, φ, θ1 , θ) is nonnegative and periodic. Recall that θ1 is the unknown tip angle. The model governed by Eq. (1) corresponds to a quasi-static loading. Suppose the beam is loaded at the tip with a fixed force, say f , and the external moment is incremented gradually 14
(a) 3.0
+1/9
y(s)
/9 +5
θ(s), rad
0.6
/9 +1
β=1
2.0 α = 5
0.8
/9 +π
2.5
(b)
Numerical Sol. Single-Piece Sol. Piecewise Sol. φ=π
1.5
0.4
1.0
/9 +5
β=1
α=5
0.2 Numerical Sol. Single Piece Sol. Piecewise Sol.
0.5 0.0 0.0
α = 10/9 β = 2/9 φ = 2π/9 0.2
0.4
s
0.6
0.8
9 +π/
φ=π
0.0 1.0
φ = 2π/9 α = 10/9
-0.5
0.0
β = 2/9
0.5
x(s)
1.0
Figure 6: Comparison of the numerical solution, the single-piece solution (M1 = 16 and M2 = M3 = 24 with s0 = 1/2) of Ref. [25], and the piecewise solution (M1i = 6 and M2i = M3i = 8 with s01 = 1/4, s2 = 1/2, and s02 = 3/4) for the beam’s (a) angular and (b) axial deflections.
(a)
(b) 3π φ+
2α2
dM < 0 ds
Γmax
0
φ
φ+π
φ
M1
Γavg =β2+2α2 cos(θ1-φ)
Γmin
dM > 0 ds
Γmax
f
π φ+
dθ ds = β dM < 0 ds
θ
dM > 0 ds
Γmax φ dM < 0 ds
Γmin
φ+ 2π
Γ
Γmin
Figure 7: Schematic plot of (a) Γ(α, β, φ, θ1 , θ), and (b) a deflected beam undergoing large deformation.
from 0 to its target value. The periodicity shown in Eq. (30) is schematically shown in Figure 7, and originates from the change of the sign of Mf , which is the moment resulting from the force f . Illustrated in Figure 8, the force f like the applied moment M1 − dM1 creates a counter clockwise moment; thus, the two are supportive. However, as the applied moment is incremented to M1 and the beam undergoes further deflection, there may be a point (x, y) on the beam’s axis through which the force’s line of action passes. Therefore, the moment resulting from the force about this point vanishes. Setting Mf in Eq. (1) to zero gives y − y1 = tan(φ), x − x1 which also follows from intuition. This phenomenon is different from that of inflection points, as the nonzero external moment M1 creates a nonzero curvature. For simplicity, we refer to this phenomenon as force’s moment stagnation. Past this stagnation, if the external moment is incremented further, say to M1 +dM1 , the applied force starts to create a clockwise moment, and thus, opposes the applied moment until another stagnation occurs. For the quasi-static loading 15
considered here, this phenomenon occurs on the clamped point, but can occur on intermediate points of the beam as well, if a sufficiently large moment is applied.
f
Mf < 0
M 1+ dM 1
f
M1
Mf = 0
φ
f Mf > 0
M 1 - dM 1 Figure 8: Quasi-static application of external moment, and the occurrence of the force’s moment stagnation.
When the force and applied tip moment are supportive, the large deflection equations are well behaved. However, when the two are in conflict, the aforementioned sign switch poses numerical difficulties. Figure 9 compares the numerical and sixth order piecewise solutions for a constant and fixed applied force as the applied moment is increased in small increments. As illustrated for α = 5 and φ = 2π/3, the numerical differential solver of Mathematica diverges as β passes 5.2, while the piecewise solutions are convergent. To validate the accuracy of the latter solutions for such cases, we substituted the piecewise solutions in Eqs. (3), (5), and (6), and monitored the forward error, which never exceeded 0.01 along the entire beam length.
(a) 0.8
α=5 φ = 2π/3
3
0.6
β=7
Numerical Sol. Piecewise Sol.
α=5 φ = 2π/3
0.4
2 +0.
θ(s), rad
4
Numerical Sol. Piecewise Sol.
y(s)
5
(b)
2
β=7
+0.2
0 0.0
0.2
β=0
1
β=0
0.0 0.2
0.4
0.6
0.8
1.0
-0.4
s
-0.2
0.0
0.2
0.4
x(s)
Figure 9: Comparison of the numerical solution and the sixth order piecewise solution (s01 = 1/4, s2 = 1/2, and s02 = 3/4) for the beam’s (a) angular and (b) axial deflections.
16
6.5. BeamSol: A Large Deflection Beam Solver To further demonstrate the effectiveness of our solution procedure, we coded our algorithm as a black-box solver in Visual C++. Our large deflection beam solver, BeamSol 1.0, provides various user-friendly functionalities, including plotting and sampling the solutions, as well as exporting them to a data file.
Case 8
Case 9
Case 3 Case 10 Case 1 Case 2 Case 5
Case 4 Case 7
Case 6
Figure 10: A snapshot of BeamSol 1.0 corresponding to the ten cases of input data given in Table 1. The annotations in red are added manually for illustration purposes.
For the special case of pure end-moment, BeamSol 1.0 employs the circular exact solution of Ref. [25], and for a general loading it employs, currently, the two-piece fifth order ATET solutions given by Eqs. (27)-(29), but is envisioned to offer various number of pieces and approximation orders in the very near future. For the solution of boundary constraints, as discussed earlier, the Newton-Raphson technique is implemented and a zero solution set is used. To validate the results of BeamSol 1.0, we ran it for many randomly generated input data, and the results compared well with the numerical differential solver of Mathematica. The executable image of BeamSol 1.0 will be posted online for free download. Table 1: Input data used for Figure 10. For all cases: l = 1 (m), I = 10−4 (m4 ), and E = 10 (MPa).
Parameter \ Case #
1
2
3
4
5
6
7
8
9
10
f (kN )
0
0
1
1
16
16
16
16
16
16
φ (radian)
–
–
2
-2
-1
-1
-1
1
1
1
M1 (kN.m)
5
-5
4
-4
4
0
-4
4
0
-4
Figure 10 depicts a snapshot of BeamSol 1.0, which plots the cartesian large deflection 17
components for ten cases that are listed in Table 1. 7. Summary and Conclusions We have presented cartesian large deflection solutions, which are independent of a beam’s arc length. The solutions are implicit and in terms of elliptic integrals, but are envisioned to guide the development of approximate closed form large deflection solutions. For the case of a large pure end-force, a new closed form solution was also presented which incorporates the rotation matrix employed for developing the aforementioned cartesian large deflection solutions. To add further efficiency and robustness to the single-piece parametric solutions, piecewise large deflection solutions have been developed which require low approximation orders for fast convergence. The proposed methodology can quantify the inherently nonlinear and complex load, deformation, and energy relationships revealed by compliant mechanisms efficiently and effectively, thus, paving the way for the emergence of novel compliant mechanisms syntheses. Approximate linear theories are still the premise of certain devices in precision engineering. However, the overall dimensions of such devices must be adjusted with the magnitudes of the applied loads. The utilization of the proposed solution methodology can help in the design of new devices with fixed size and no dependence on the load magnitudes. We coded the solutions as a stand-alone black-box solver, which offers several functionalities for the solution of large deflection of EulerBernoulli cantilever beams. The executable will be posted for free download for those who are interested in beam applications. Acknowledgments The authors are grateful to the anonymous reviewers for the mindful comments. [1] G. K. Bennig, Atomic force microscope and method for imaging surfaces with atomic resolution, U.S. Patent 4,724,318, 1988. [2] D. J. Leo, Engineering Analysis of Smart Material Systems, Wiley InterScience, Wiley, 2007. [3] L. L. Howell, Compliant Mechanisms, John Wiley & Sons, 2001. [4] R. H. Burns, The Kinetostatic Synthesis of Flexible Link Mechanisms, Ph.D. Dissertation, Yale University, 1964. [5] R. H. Burns, F. R. E. Crossley, Kinetostatic Synthesis of Flexible Link Mechanisms, ASME Mechanisms Conference, ASME Paper No 68–MECH–36, 1968. 18
[6] S. J. Winter, T. E. Shoups, The displacement analysis of path-generating flexible-link mechanisms, Mechanism and Machine Theory 7 (4) (1972) 443–451. [7] C. Huang, B. Roth, Dimensional Synthesis of Closed-Loop Linkages to Match Force and Position Specifications, ASME Journal of Mechanical Design 115 (2) (1993) 194–198. [8] V. Krovi, G. K. Ananthasuresh, V. Kumar, Kinematic and Kinetostatic Synthesis of Planar Coupled Serial Chain Mechanisms, Journal of Mechanical Design 124 (2) (2002) 301–312. [9] M. H. F. Dado, Limit position synthesis and analysis of compliant 4-bar mechanisms with specified energy levels using variable parametric pseudo-rigid-body model, Mechanism and Machine Theory 40 (8) (2005) 977–992. [10] H. Tari, H.-J. Su, A Complex Solution Framework for the Kinetostatic Synthesis of a Compliant Four-Bar Mechanism, Mechanism and Machine Theory 46 (8) (2011) 1137–1152. [11] H. Tari, H.-J. Su, J. D. Hauenstein, Classification and complete solution of the kinetostatics of a compliant Stewart–Gough platform, Mechanism and Machine Theory 49 (0) (2012) 177–186. [12] J. Argyris, R. Kaˇcianauskas, Semi-analytical finite elements in the higher-order theory of beams, Computer Methods in Applied Mechanics and Engineering 138 (1–4) (1996) 19–72. [13] J. Mazars, P. Kotronis, F. Ragueneau, G. Casaux, Using multifiber beams to account for shear and torsion: Applications to concrete structural elements, Computer Methods in Applied Mechanics and Engineering 195 (52) (2006) 7264–7281. [14] L. Jun, H. Hongxing, S. Rongying, Dynamic finite element method for generally laminated composite beams, International Journal of Mechanical Sciences 50 (3) (2008) 466–480. [15] T. P. Vo, H.-T. Thai, Vibration and buckling of composite beams using refined shear deformation theory, International Journal of Mechanical Sciences 62 (1) (2012) 67–76. [16] E. O. nate, A. Eijo, S. Oller, Simple and accurate two-noded beam element for composite laminated beams using a refined zigzag theory, Computer Methods in Applied Mechanics and Engineering 213–216 (0) (2012) 362–382. [17] A. J. Sadowski, J. M. Rotter, Solid or shell finite elements to model thick cylindrical tubes and shells under global bending, International Journal of Mechanical Sciences 74 (0) (2013) 143–153.
19
[18] F. Miao, G. Sun, K. Chen, Transient response analysis of balanced laminated composite beams by the method of reverberation-ray matrix, International Journal of Mechanical Sciences 77 (0) (2013) 121–129. [19] Z. Zhang, X. Huang, Z. Zhang, H. Hua, On the transverse vibration of Timoshenko doublebeam systems coupled with various discontinuities, International Journal of Mechanical Sciences 89 (0) (2014) 222–241. [20] C. W. Lim, M. Z. Islam, G. Zhang, A nonlocal finite element method for torsional statics and dynamics of circular nanostructures, International Journal of Mechanical Sciences 9495 (0) (2015) 232–243. [21] K. E. Bisshop, D. C. Drucker, Large Deflection of Cantilever Beams, Quarterly of Applied Mathematics 3 (1945) 272–275. [22] L. L. Howell, The Design and Analysis of Large-Deflection Members in Compliant Mechanisms, M.S. Thesis, Purdue University, 1991. [23] A. Banerjee, B. Bhattacharya, A. Mallik, Large deflection of cantilever beams with geometric non-linearity: Analytical and numerical approaches, International Journal of Non–Linear Mechanics 43 (5) (2008) 366–376. [24] L. Chen, An integral approach for large deflection cantilever beams, International Journal of Non–Linear Mechanics 45 (3) (2010) 301–305. [25] H. Tari, On the parametric large deflection study of Euler-Bernoulli cantilever beams subjected to combined tip point loading, International Journal of Non-Linear Mechanics 49 (2013) 90–99. [26] F. A. Chouery, Exact and Numerical Solutions for Large Deflection of Elastic Non-Prismatic Beams, FAC Systems INC., WA .
20
Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded Euler-Bernoulli Cantilever Beams Hafez Tari, G.L. Kinzel, D.A. Mendelsohn
www.elsevier.com/locate/ijmecsci
PII: DOI: Reference:
S0020-7403(15)00244-1 http://dx.doi.org/10.1016/j.ijmecsci.2015.06.024 MS3034
To appear in:
International Journal of Mechanical Sciences
Received date: 2 May 2015 Revised date: 7 June 2015 Accepted date: 28 June 2015 Cite this article as: Hafez Tari, G.L. Kinzel, D.A. Mendelsohn, Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded EulerBernoulli Cantilever Beams, International Journal of Mechanical Sciences, http://dx. doi.org/10.1016/j.ijmecsci.2015.06.024 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Cartesian and Piecewise Parametric Large Deflection Solutions of Tip Point Loaded Euler-Bernoulli Cantilever Beams Hafez Tari1 , G.L. Kinzel2 , D.A. Mendelsohn3 Department of Mechanical Engineering, The Ohio State University, Columbus, OH 43210, USA
Abstract This paper focuses on the large deflection solutions of Euler-Bernoulli cantilever beams loaded at the tip with a force and moment. Cartesian large deflection solutions are given for general tip point loading, as well as a new closed form solution for the case of a large end-force only. Additionally, we develop piecewise parametric large deflection solutions, which can efficiently and effectively quantify the inherently nonlinear and complex load, deformation, and energy relationships exhibited by compliant mechanisms, thus, paving the way for the emergence of novel compliant mechanism syntheses. The utilization of the proposed solution methodology can help the precision engineering community in the design of new devices with fixed size and no dependence on the load magnitudes. We coded the solutions as BeamSol, a stand-alone user-friendly black-box solver, which is envisioned to help researchers in engineering analyses and syntheses of beam applications. Keywords: Large Deflection, Euler-Bernoulli Beam, Automatic Taylor Expansion Technique, Piecewise Solution, BeamSol. 1. Introduction Beams have found many applications in engineering, for instance, in atomic force microscopy [1], and design of actuators and energy harvesters in smart materials [2]. In particular, slender beams are a building-block of compliant mechanisms [3]. Such mechanisms transmit or transform motion, load, or energy via the mobility and flexibility of their members. In contrast to traditional or rigid mechanisms, which are made strong through increasing their stiffness, compliant mechanisms are designed to be flexible and strong. In contrast to their rigid counterparts, compliant mechanisms offer part-count reduction; reduced assembly time; ease of manufacturing 1
Email: [email protected], address all correspondence to this author, Tel./Fax: (614) 292-2289/3361 Email: [email protected] 3 Email: [email protected] 2
Preprint submitted to International Journal of Mechanical Sciences
July 6, 2015
processes; increased reliability and precision; wear, weight and maintenance reduction. Despite their advantages, compliant mechanisms exhibit a complex nonlinear behavior, which poses modeling and computational challenges. “Kinetostatics” alleviates this shortcoming by replacing beams in a compliant mechanism with equivalent pseudo-rigid-body models, which consist of rigid segments connected with rotational springs that compensate for the compliance. This approach is advantageous, as the traditional kinematic modeling techniques can be applied to find the geometric dimensions (e.g., segments length) and compliance (e.g., springs stiffness) for prescribed deformation/load/energy constraints. Utilizing kinetostatics, Burns [4], and Burns and Crossley [5] synthesized a four-bar linkage with a flexible coupler link. Winter and Shoups [6] studied the kinetostatic analysis of path-generating flexible-link mechanisms. Extending the classical dimensional synthesis problems, Huang and Roth [7] formulated the kinetostatic synthesis equations for three types of planar linkages with linear springs, and obtained the solutions by resorting to polynomial solvers. Krovi et al. [8] studied the kinetostatic synthesis of planar coupled serial chain mechanisms using a hybrid approach, which combines precision point synthesis and optimization. Dado [9] developed a variable parametric pseudo-rigid-body model for limit position synthesis of compliant four-bar linkage with energy requirements. These attempted to find a unique or a partial set of pseudo-rigid-body model solutions. However, Tari and Su [10] recently proposed a systematic framework that transforms the modeling equations into a system of polynomials, which can be solved for the entire solution set. A parameter continuation was devised to avoid computation of junk solutions for subsequent syntheses. Later, Tari et al. [11] showed that this framework is effective even for the forward analysis of a 3D compliant Stewart-Gough platform. Extensive work has been done on linear beam theory to obtain explicit solutions for small deformation of cantilever beams subjected to a tip point loading. Such solutions are still the premise of certain devices in precision engineering. However, the overall dimensions of such devices must be adjusted with the magnitudes of the applied loads, as linear theory solutions are not valid once the beam’s axis slope is significant. Despite the success of pseudo-rigid-body models, they are valid for a limited range of motion. When the compliant members undergo large deformations, large deflection solutions must be used for accurate modeling. Finite-element solution procedures can overcome these shortcomings. For instance, building on independent approximations of three-dimensional fields of displacements and stresses, Argyris and Kaˇcianauskas [12] presented semi-analytical finite-elements leading to the development of a higher-order theory for beams. Mazars et al. [13] employed higher order interpolation functions to investigate solutions for an enhanced multifiber beam element accounting for shear and torsion. Their constitutive 2
model avoids shear locking phenomena, and also captures nonlinear cross section warping. Jun et al. [14] introduced a dynamic finite-element method for free vibration analysis of generally laminated composite beams. Their model is based on first-order shear deformation theory and accounts for Poisson effect, couplings between extensional, bending and torsional deformations, shear deformation and rotary inertia. Vo and Thai [15] analyzed vibration and buckling of composite beams with arbitrary lay-ups using refined shear deformation theory that accounts for the parabolical variation of shear strains through the depth of the beam. Utilizing classical Timoshenko beam theory and refined zigzag kinematics, O˜ nate et al. [16] presented a simple linear two-noded beam element sufficient for analysis of composite laminated and sandwich beams. Sadowski and Rotter [17] studied the use of solid continuum finite-elements and shell finiteelements in the modeling of nonlinear plastic buckling behavior of cylindrical metal tubes and shells under global bending. Miao et al. [18] enhanced the reverberation-ray matrix approach to analyze the transient response of a laminated composite beam subjected to an impulse force load. Zhang et al. [19] developed an efficient analytic framework to study transverse vibrations of two parallel Timoshenko beams that are connected by discrete springs and coupled with various discontinuities. Very recently, Lim et al. [20] employed the nonlocal elasticity stress theory and investigated the torsional static and dynamic effects of circular nanostructures subjected to concentrated and distributed torques. The foregoing finite-element solution methodologies are robust for arbitrary beam applications. However, efficient solution methods have been desired for special loading and beam cross sectional geometries and properties. For example, Bisshop and Drucker [21] reported exact implicit elliptic integral solutions of an Euler-Bernoulli cantilever beam of uniform cross section and loaded with a tip point vertical force. Howell [22] generalized these solutions for the case of general combined tip point loading. Nonetheless, the implicitness of the foregoing elliptic integral solutions poses numerical evaluation challenges. To alleviate this drawback, Banerjee et al. [23] employed nonlinear shooting and Adomian’s decomposition methods, and reported numerical and approximate large deflection solutions for concentrated intermediate and endpoint loadings. Chen [24] developed a numerical integral approach for complex loading and variable beam’s cross sectional properties. Employing an Automatic Taylor Expansion Technique (ATET), Tari [25] recently developed parametric large deflection solutions for Euler-Bernoulli cantilever beams loaded at the tip with force and moment. These solutions are valid for the entire beam length, and capture both angular and axial inflection points. We refer to these solutions as single-piece parametric solutions, as they are expanded about a single start point for the whole computational domain. Even though this treatment was shown to be effective, rather 3
high approximation orders are needed for large loading inputs. Unfortunately, the higher the approximation orders, the more the solution evaluation time, and more importantly the more the number of extraneous solutions. In this paper, we develop cartesian large deflection solutions, which are independent of beam’s arc length for applied tip force and moment. Also, we develop explicit closed form solutions for the case of a large pure end-force. In addition, we develop piecewise parametric large deflection solutions, which outperform the efficiency and robustness of the single-piece parametric solutions. The proposed large deflection solutions can help in the design of new devices with fixed size and no dependence on the load magnitudes. Our solutions are coded as a stand-alone user-friendly black-box solver, BeamSol, which is envisioned to facilitate analysis and synthesis of of novel compliant mechanisms. The rest of the paper is organized as follows. Next, we review Euler-Bernoulli beams and their mathematical modeling. Section 3 gives cartesian large deflection solutions, and Section 4 gives the closed form solutions for a large pure end-force. Section 5 presents the piecewise parametric solutions together with BeamSol, and is followed by results and discussions. Finally, summary and conclusions are given. 2. A Brief Review on Euler-Bernoulli Beams A schematic view of a cantilever beam subjected to a tip moment M1 and force f with the inclination φ is depicted in Figure 1, in which E and I are the Young’s modulus and the moment of inertia of the beam’s cross section; s and θ are the arc length and the angular deflection of the beam; and x and y are the horizontal and vertical coordinates of the deflected beam, respectively.
x1
M1
Inflection(κ=0)
x
f
θ
s
φ
θ1
y1
y l, E, I
Undeflected
Figure 1: A schematic view of a tip point loaded cantilever beam.
A common beam type is an Euler-Bernoulli beam, whose large deflection is governed by
κ(x) =
M1 + f sin(φ)(x1 − x) − f cos(φ)(y1 − y) M d2 y/dx2 = , = 2 3/2 EI EI [1 + (dy/dx) ]
4
(1)
where κ is the curvature, and M is the bending moment. The boundary conditions for the cantilever beam are dy (0) = 0. dx
y(0) = 0,
(2)
The behavior of the Euler-Bernoulli beam depends on the choice of the material (E), the geometry of the cross section (I), and even the location of the applied load (M ). It is a common practice to nondimensionalize the parameters as s/l → s, x/l → x, and y/l → y; and finally to parameterize the problem with respect to the beam’s arc length. For the case of a beam with constant cross sectional properties, i.e., constant E and I, the resulting equations are (see, for example, Tari [25])
d2 θ(s) = α2 sin(θ(s) − φ), ds2 dθ(1) = β, θ(0) = 0, ds dx(s) = cos(θ(s)), x(0) = 0, ds dy(s) = sin(θ(s)), y(0) = 0, ds where α2 =
f l2 EI
M1 l EI .
and β =
(3) (4) (5) (6)
The exact solutions to this weakly decoupled boundary value
problem can be written in terms of the solutions of the following implicit elliptic integral relations s=± x=±
θ(s) 0 θ(s) 0
y=±
0
θ(s)
dθ β 2 + 2α2 (cos(θ1 − φ) − cos(θ − φ)) cos(θ)dθ β2
+ 2α2 (cos(θ
1
− φ) − cos(θ − φ))
,
(7)
,
(8)
,
(9)
sin(θ)dθ β 2 + 2α2 (cos(θ1 − φ) − cos(θ − φ))
where the plus and minus signs in Eqs. (7)-(9) are inherent to the beam’s multiple behavior, and often in the literature are denoted as elbow up and elbow down solutions. The elliptic integral relations pose some challenges. First, there is an implicit dependence on the unknown tip rotation θ1 , which complicates the integration procedure, even when solved numerically. Second, depending on input loading parameters, the solutions may be very close to singularities, which drastically complicates the numerical integration. Recall that the signs in the elliptic integral relationships change as the input loading creates inflection points, where the curvature changes sign, such as the one schematically shown in Figure 1. As the number of inflection
5
points increases, so does the number of solution branches. Therefore, obtaining solutions that do not share such difficulties is desirable. 3. Cartesian Large Deflection Solutions The parametric large deflection solutions (8) and (9) are in terms of the intermediate parameter s, which is not explicitly at hand. The goal of this section is to eliminate s from the solutions, and obtain a cartesian large deflection solution, i.e., the vertical large deflection component directly in terms of the horizontal large deflection component or vice versa. To this end, we employ the coordinates transformation ⎡ ⎣
⎤
R(φ)
⎡
cos(φ)
X
⎤ ⎡ sin(φ)
⎤ x
⎦ ⎣ ⎦, ⎦= 1 ⎣ l − sin(φ) cos(φ) y Y
(10)
where R(φ) is the orthogonal rotation matrix, which rotates the components x and y clockwise by the angle φ about the origin. Keeping in mind that a rotation does not change the curvature, Eq. (1) simplifies to d2 Y /dX 2 = α2 Y + α1 , [1 + (dY /dX)2 ]3/2
(11)
where α1 = β−α2 (cos(φ)y1 −sin(φ)x1 )/l, α and β are as before, and the new boundary conditions become
Y (0) = 0,
dY (0) + tan(φ) = 0. dX
(12)
We can integrate Eq. (11) once, and employ the second boundary condition of Eq. (12) to write
1 1 + (dY /dX)2
= | cos(φ)| − α1 Y −
α2 2 Y G(Y ). 2
(13)
Then, integrating the foregoing equation and applying the first boundary condition of Eq. (12), we write X=
Y (X)
0
G(Y ) dY. 1 − G(Y )2
(14)
Equation (14) is not amenable to an explicit closed form solution in terms of elementary functions. However, letting G1 = G(Y (1)) and simplifying the results, a standard form of elliptic integral solution may be obtained as 6
X=±
G(Y )
| cos(φ)|
1 G √ dG. 2 2 1−G β − 2α2 (G − G1 )
(15)
Had we replaced R(φ) with its transpose (RT (φ)) in the coordinates transformation (10), a similar result would have been obtained. However, both solutions are qualitatively similar, and have no major advantage over each other. It is worth mentioning that Ref. [26] also reported a solution procedure for obtaining a cartesian large deflection solution. However, their coordinates transformation is not based on rotation and is somewhat more complex than Eq. (15). The inextensibility of the beam imposes a side condition, which would further constrain the trajectory of the beam’s endpoints as 1=±
G1 | cos(φ)|
√
1 1 dG. 2 2 1−G β − 2α2 (G − G1 )
(16)
It is clear that the new cartesian elliptic integral large deflection solution (15) retains the main shortcomings of the elliptic large deflection solutions (8) and (9). Even though the new solution involves only x and y, it is implicit and highly nonlinear, thus, still requires numerical evaluation techniques. Likewise, the new solution still has two branches (plus and minus solutions), and the identification of the proper branch for any input loading parameters requires careful assessment. The square root in Eq. (13) is bounded below by one, thus, we have 0 < G(Y ) ≤ 1. This paves the way for fairly accurate approximations through neglecting higher powers of the quadratic function G(Y ) in the aforementioned formulae. Such assumptions facilitate obtaining efficient approximate closed form solutions in terms of elementary functions. Nonetheless, we aim to obtain solutions, which are accurate for arbitrary applied loadings. 4. Closed Form Approximate Large Deflection Solutions for Large Pure End-forces Explicit exact large deflection solutions exist for only a pure end-moment loading (zero force, α = 0); see, for example, Ref. [25]. However, no explicit solution is reported in the literature for a general pure end-force loading (zero moment, β = 0). Nonetheless, for a sufficiently large applied end-force with no moment, the results of Ref. [25] suggest that beam’s end angle approaches the force angle, i.e., θ1 ≈ φ. Utilizing this result, Eq. (7) can be simplified and solved analytically. For β = 0 and θ1 = φ, we give the exact solution of Eqs. (3) and (4) as φ θ(s) = φ − 4 arctan(tan( ) exp(−αs)), 4
(17)
which is valid for both elbow up and down solution branches. To obtain the corresponding axial deflection components, we integrate Eqs. (8) and (9), and used the trigonometric identities 7
⎧ ⎨ cos(4 arctan(z)) = 1 − 8 z 22 2 , (1+z ) ⎩
sin(4 arctan(z)) =
2
1−z 4z (1+z 2 )2 ,
to write
RT (φ)
⎤ ⎡ cos(φ) − sin(φ) ⎢ s − x(s) ⎦⎣ ⎣ ⎦ =⎣ sin(φ) cos(φ) y(s) ⎡
⎤
⎡
4 α 4 α
1
1+tan2 (φ/4) exp(−2αs) tan(φ/4) exp(−αs) 1+tan2 (φ/4) exp(−2αs)
−
1 1+tan2 (φ/4)
−
tan(φ/4) 1+tan2 (φ/4)
⎤ ⎥ ⎦,
(18)
in which the manifestation of the rotation matrix R(φ) is in analogy with the cartesian solutions given in the previous section. (a)
4
(b)
α = 24
1.0
0.8
3
0.6
y(s)
θ(s), rad
Numerical Sol. Explicit Sol.
−π/10
−2
Numerical Sol. Explicit Sol.
φ=π
2
0.4 1
0.2 +2
0 0.0
α=6 0.2
0.4
0.6
φ = π/10 0.8
φ=π
0.0 1.0
-1.0
-0.5
s
α = 24
α=6 0.0
0.5
+π/10
φ = π/10 1.0
x(s)
Figure 2: Comparison of the numerical solution and the explicit solutions (17) and (18) for the beam’s (a) angular and (b) axial deflections.
Figure 2 compares the explicit solutions, given by Eqs. (17) and (18), with numerical solutions. For the latter, the numerical differential solver of Mathematica is utilized to solve Eqs. (3) and (4). Then, the obtained numerical angular deflections are substituted into Eqs. (5) and (6), which are integrated numerically to determine the axial deflection components. Therefore, the numerical solution is convergent for only the two smallest load magnitudes (α = 6 and α = 8), and diverges at the vicinity of the clamped point as the load magnitude increases. In contrast, the explicit solutions are convergent for the entire beam length. Because the tip angle converges to the force angle for every case studied (which was the derivation assumption for the explicit solutions), the obtained convergent solutions are, in fact, accurate. The divergent numerical solutions are not shown when plotting the axial solutions in Subfigure 2(b). Note that there is an outrageous discrepancy in the form of the exact circular large deflection solution, given in Ref. [25], and the newly developed closed form approximate large deflection
8
solution. This underlines the complexity of the beam behavior and the need for a robust solution procedure. 5. Piecewise Parametric Large Deflection Solutions Letting ψ(s) θ(s) − φ, we can rewrite Eqs. (3)-(6) as
d2 ψ(s) = α2 sin(ψ(s)), ds2 dψ(1) = β, ψ(0) = −φ, ds dx(s) = cos(ψ(s) + φ), x(0) = 0, ds dy(s) = sin(ψ(s) + φ), y(0) = 0, ds
(19) (20) (21) (22)
which are numerically well-conditioned. For i = 1, . . . , n, let [si , si+1 ] be a subdivision of the computational domain s ∈ [0, 1] such that ∪ni=1 [si , si+1 ] = [0, 1], where s1 = 0 and sn+1 = 1. In addition, let θi (s), xi (s), and yi (s) denote the deflection components corresponding to s ∈ [si , si+1 ]. Employing the automatic Taylor expansion technique [25], the parametric large deflection solutions corresponding to each subdivision are written as expansions about the point s0i as 2
θi (s) ≈ φ + c0i + c1i (s − s0i ) + α
M 1i j=2
xi (s) ≈ x0i +
M 2i j=1
and yi (s) ≈ y0i +
M 3i j=1
1 dj−2 (sin(ψ(s))) j! dsj−2
1 dj−1 (cos(ψ(s) + φ)) j! dsj−1 1 dj−1 (sin(ψ(s) + φ)) j! dsj−1
where s0i ∈ [si , si+1 ], c0i = ψi (s0i ), c1i =
dψi (s0i ) , ds
(s − s0i )j ,
(23)
s=s0i
(s − s0i )j ,
(24)
(s − s0i )j ,
(25)
s=s0i
s=s0i
x0i = xi (s0i ), y0i = yi (s0i ), and M1i , M2i ,
and M3i are the approximation orders. Note that all of the summation terms are at hand, as the higher order derivatives of ψ(s) can be found in terms of c0i and c1i upon the successive differentiation of Eq. (19). Accordingly, there are 4n unknowns, i.e., C4n = {c0i , c1i , x0i , y0i }, which should be identified to fully define the piecewise parametric large deflection solutions. To this end, requiring the continuity in the zeroth and first order derivatives of the piecewise angular solutions at the interfaces, and continuity in piecewise axial deflections, together with the original boundary conditions (20)-(22) constitute the 4n algebraic-trigonometric constraints
9
B(C4n ) =
⎧ ⎪ ⎪ θ1 (0) = 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n ⎨ dθ ds (1) = β
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬
⎪ ⎪ ⎪ x1 (0) = 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ y1 (0) = 0
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
∪
⎧ ⎪ ⎪ θj (sj+1 ) = θj+1 (sj+1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dθj+1 j ⎪ ⎨ dθ ds (sj+1 ) = ds (sj+1 ) ⎪ ⎪ ⎪ xj (sj+1 ) = xj+1 (sj+1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ yj (sj+1 ) = yj+1 (sj+1 )
,
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ j = 1, . . . , n − 1 , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭
(26)
which may be solved for the entire unknowns. For instance, taking n = 2 (i.e., i = 1, 2) and M1i = M2i = M3i = 5, the fifth order piecewise ATET solutions for θi (s), xi (s), and yi (s) may be written as ⎧ ⎧ ⎨ φ + c01 + c11 (s − s01 ) + αs1 (s − s01 )2 + c11 αc1 (s − s01 )3 ⎪ ⎪ 1! 2! 3! ⎪ ⎪ , 0 ≤ s ≤ s2 , ⎪ c11 (αc1 η21 −3α4 ) ⎪ αs1 η11 ⎩ 4 5 ⎪ + 4! (s − s01 ) + (s − s ) ⎨ 01 5! θ(s) = ⎧ ⎪ ⎪ ⎨ φ + c02 + c12 (s − s02 ) + αs2 (s − s02 )2 + c12 αc2 (s − s02 )3 ⎪ ⎪ 1! 2! 3! ⎪ ⎪ , s2 ≤ s ≤ 1, ⎪ c12 (αc2 η22 −3α4 ) ⎩ ⎩ αs2 η12 4 5 + 4! (s − s02 ) + (s − s ) 02 5! x(s) = ⎧⎧ c01 ) c01 ) ⎪ ⎪ (s − s01 ) − c11 sin(˜ (s − s01 )2 x01 + cos(˜ ⎪⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ 2 ⎪ ⎪ − c11 cos(˜c01 )+αs1 sin(˜c01 ) (s − s01 )3 − c11 (3αs1 cos(˜c01 )+η11 sin(˜c01 )) (s − s01 )4 ⎪ ⎪ 3! 4! ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ α η sin(˜ c )+η cos(˜ c ) ⎩ s1 31 01 41 01 ⎨ − (s − s01 )5 , 5! ⎧ c02 ) c02 ) ⎪ ⎪ ⎪ ⎪ (s − s02 ) − c12 sin(˜ (s − s02 )2 x02 + cos(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ c2 cos(˜ c02 )+αs2 sin(˜ c02 ) ⎪ ⎪ − 12 (s − s02 )3 − c12 (3αs2 cos(˜c024!)+η12 sin(˜c02 )) (s − s02 )4 ⎪ 3! ⎪ ⎪ ⎪ ⎪ ⎪ ⎩⎪ ⎩ − αs2 η32 sin(˜c025!)+η42 cos(˜c02 ) (s − s02 )5 ,
(27)
0 ≤ s ≤ s2 ,
s2 ≤ s ≤ 1, (28)
and
y(s) = ⎧⎧ c01 ) c01 ) ⎪ ⎪ ⎪ (s − s01 ) + c11 cos(˜ (s − s01 )2 y01 + sin(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ ⎪ c2 sin(˜ c01 )−αs1 cos(˜ c01 ) ⎪ (s − s01 )3 − c11 (3αs1 sin(˜c014!)−η11 cos(˜c01 )) (s − s01 )4 ⎪ − 11 ⎪ 3! ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨⎩ + αs1 η31 cos(˜c015!)−η41 sin(˜c01 ) (s − s01 )5 , ⎧ c02 ) c02 ) ⎪ ⎪ ⎪ (s − s02 ) + c12 cos(˜ (s − s02 )2 ⎪ y02 + sin(˜ ⎪ ⎪ 1! 2! ⎪ ⎨ ⎪ ⎪ c2 sin(˜ c02 )−αs2 cos(˜ c02 ) ⎪ ⎪ − 12 (s − s02 )3 − c12 (3αs2 sin(˜c024!)−η12 cos(˜c02 )) (s − s02 )4 ⎪ 3! ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎩ ⎩ + αs2 η32 cos(˜c025!)−η42 sin(˜c02 ) (s − s02 )5 ,
0 ≤ s ≤ s2 ,
s2 ≤ s ≤ 1, (29)
where c˜0i = c0i + φ, αsi = α2 sin(c0i ), αci = α2 cos(c0i ), η1i = αci − c21i , η2i = 4αci − c21i , η3i = αci − 7c21i , and finally η4i = c21i η2i + 3α2si for i = 1, 2. Note that s01 , s02 , and s2 are to be chosen 10
as input parameters, thus, leaving only eight unknowns C8 = {c01 , c11 , c02 , c12 , x01 , x02 , y01 , y02 }, which would be identified upon imposing the boundary conditions B(C8 ), an adaptation of B(C4n ) for this special case. 5.1. Choice of the Number of Subdivisions n The more the number of subdivisions, the more accurate the piecewise solutions. However, more subdivisions introduce a significant number of intermediate unknowns, which would burden the computation. Nevertheless, an increase in the number of subdivisions reduces the need for large approximation orders, which are generally required for attaining a significant accuracy. To get a good balance between solution accuracy and computational cost, we set n = 2, for which the resulting large deflection solutions are highly accurate and compact for low approximation orders. 5.2. Choice of the Segments’ Endpoints si and si+1 Recall that the premise of the piecewise solutions is the Taylor series expansion, which approximates the unknown exact large deflection solutions with polynomials in the arc length of the beam. It is clear that the smaller the computational domain, i.e., the length of the arc length for each piecewise solution, the better the accuracy of the approximations. However, for a fixed number of subdivisions and fixed approximations orders, improving the accuracy of each piecewise solution requires shortening the length of that segment, which would force compromising the accuracy of the remaining segment solutions. Therefore, to avoid this, we subdivide the beam length equally, i.e., [si , si+1 ] = [(i − 1), i]/n for i = 1, . . . , n. 5.3. Choice of the Start Points s0i The piecewise solutions also depend on start points s0i , which must be chosen. Like the parametric solution presented in Ref. [25], start points can simplify the solution procedure. For example, the unknowns {c01 , c1n , x01 , y01 } would become known parameters, if s01 and s0n are taken as s1 = 0 and sn+1 = 1, respectively. In contrast, setting s0i to the center of each segment provides the most accuracy and convergence. 5.4. Independence of the Approximation Orders Any of the approximation orders M1i , M2i , and M3i can be chosen independently. This is significant in two ways. The accuracy of any of the angular and axial deflection components can be exclusively adjusted without disturbing the accuracy of the remaining solution components. Likewise, the accuracy of each segment solution can be adapted without affecting the accuracy
11
of the remaining segment solutions. As an example for the piecewise solutions given by Eq. (27), the accuracy of the angular solution on the domain of s ∈ [0, s2 ] can be controlled independently of that of s ∈ [s2 , 1], and vice versa. 5.5. Solution of the Boundary Constraints B(C4n ) The system of boundary constraints B(C4n ) (26) is a decoupled system of algebraic-trigonometric equations. To solve it efficiently, the angular constraints can be solved for c0i and c1i independently of the axial constraints. Then, the resulting solutions can be inserted into the horizontal and vertical deflection constraints so as to solve for x0i and y0i , respectively. One major advantage of the piecewise solution procedure over the single-piece solution approach of Ref. [25] is that the boundary constraints, i.e., B(C4n ), of the former are far better conditioned than those of the latter. For the single-piece solution, the boundary constraints can be converted into a system of polynomial equations, and solved for the entire solution set, among which the true solutions can be sought. However, this approach is inefficient, and impractical for daily engineering applications. In contrast, we employ the Newton-Raphson technique, which initiates at a given start solution and iteratively updates it to converge to a feasible solution. Even though multiple mathematical solutions may exist, we have found that the Newton-Raphson technique, starting at a zero solution set, has always converged to the true physical solution for B(C4n ). This feature has led us to coding the algorithm as a black-box simulation software, which is discussed in the next section. 6. Results and Discussions For the simulations given in the following, we take advantage of the symmetry of the loading, and assume that the force angle is bounded in the first and second quadrants, i.e., φ ∈ [0, π]. Also, with the piecewise solutions, we mean the two-piece solutions. Additionally, for the sake of simplicity, we plot the nondimensional deflection solutions. Finally, for validation purposes, numerical solutions of the exact governing equations were obtained from the numerical differential solver of Mathematica. 6.1. Effect of the Loading Parameters In general, the larger the load magnitudes, the more the deflection revealed by the beam. Figure 3 illustrates that the compact fifth order piecewise solutions in Eqs. (27)-(29) are accurate for a wide range of the input parameters, namely, α, β ∈ [0, 3] for arbitrary φ. In particular, the
12
axial deflections corresponding to the case of α = β = 3 and φ = π, shown in Subfigure 3(b), are rather large which are captured well with the fifth order of approximation. (a)
4
(b)
0.8
Numerical Sol. Piecewise Sol.
φ=π β=3
0.6
0 /1 +π
/10 +3
y(s)
α=3
/10 +3
θ(s), rad
3
2
+3 /10
0.4
φ=π
β=3
+3/10
α=3
0 0.0
α = 3/10 β = 3/10 φ = π/10 0.2
0.4
0.6
0.8
+π/10
Numerical Sol. Piecewise Sol.
0.2
1
α = 3/10
0.0
1.0
-0.4
-0.2
0.0
0.2
s
0.4
β = 3/10
0.6
φ = π/10 0.8
1.0
x(s)
Figure 3: Comparison of the numerical solution and the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4 for the beam’s (a) angular and (b) cartesian deflections.
The magnitude of the force angle also plays a crucial role in the beam’s deflection. For the case of pure force loading, it is intuitive that the more obtuse the force angle, the more the obtained deflection. As shown schematically in Subfigure 4(a), when the force angle is constant, the tip point traces an elliptic-like curve as the force magnitude is increased. In fact, as shown in Subfigure 4(b), the beam’s tip point locus is almost circular for sufficiently small force magnitudes, but diverges to an elliptic curve for large force magnitudes. It can be seen that with the absence of the external moment (β = 0), the fifth order piecewise solutions are valid even for a wide range, i.e., α ∈ [0, 4].
(a)
(b)
1.0
1.0
Inc
se
0.8
in
0.6
α
0.6
0.4
0.4
0.2
0.2
0.0
φ = 5π/6
φ = 2π/3
α=0 -0.5
0.0
x1
0.5
φ = π/2
φ = π/3 φ = π/6
φ=π
y1
y1
0.8
rea
Numerical Sol. Piecewise Sol. β=0
0.0 1.0
-0.5
0.0
x1
0.5
1.0
Figure 4: Beam’s endpoint trajectory: (a) for a typical pure end-force, and (b) for α ∈ [0, 4] with several load angles obtained using the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4.
6.2. The Post-Buckling Behavior For a pure compressive force loading (β = 0 and φ = π), it is quite well known that the beam starts to buckle at the critical force magnitude α = π/2. However, linear beam theory fails 13
to predict beam post-buckling behavior. Shown in Figure 5, the fifth order piecewise solutions are intelligent in recognizing the critical force magnitude, and also give accurate large deflection solutions past this point. (a)
3.0 2.5
(b)
Numerical Sol. Piecewise Sol.
β=0 φ=π
0.8
α = π/2+1.5 0.6
1.5
y(s)
θ(s), rad
2.0
α = π/2+1.5
0.4
1.0
α ε [0,π/2] 0.2
0.4
0.6
0.8
Numerical Sol. Piecewise Sol.
β=0 φ=π
+0.1
.1 +0
α = π/2+0.1
0.5 0.0 0.0
0.2
α = π/2+0.1
0.0 1.0
-0.4
-0.2
0.0
s
0.2
0.4
0.6
α ε [0,π/2] 0.8
1.0
x(s)
Figure 5: Comparison of the numerical solution and the fifth order piecewise solutions (27)-(29) with s01 = 1/4, s2 = 1/2, and s02 = 3/4 for the beam’s post-buckling (a) angular and (b) cartesian deflections.
6.3. Efficiency of the Piecewise Solutions Figure 6 compares the proposed piecewise solution procedure against the single-piece solution procedure of Ref. [25]. As demonstrated, to get accurate angular large deflection solutions, the piecewise solutions require only sixth order of approximation, while the single-piece solutions require 16th order of approximation for the chosen loading parameters. Even though cartesian large deflection solutions require, in general, higher approximation orders than the angular solutions, this increase is minimal for the piecewise solutions. The axial single-piece solutions require 24th order of approximation, while the cartesian piecewise solutions require only eighth order of approximation. Obviously, the approximation order of the piecewise solutions can be further reduced, if the beam is subdivided into more segments. 6.4. Robustness of the Piecewise Solutions The beam’s angular derivative may be obtained by integrating Eq. (3) which after utilizing Eq. (4) can be written as
dθ(s) ds
2
= β 2 + 2α2 [cos(θ1 − φ) − cos(θ(s) − φ)] Γ(α, β, φ, θ1 , θ),
(30)
where Γ(α, β, φ, θ1 , θ) is nonnegative and periodic. Recall that θ1 is the unknown tip angle. The model governed by Eq. (1) corresponds to a quasi-static loading. Suppose the beam is loaded at the tip with a fixed force, say f , and the external moment is incremented gradually 14
(a) 3.0
+1/9
y(s)
/9 +5
θ(s), rad
0.6
/9 +1
β=1
2.0 α = 5
0.8
/9 +π
2.5
(b)
Numerical Sol. Single-Piece Sol. Piecewise Sol. φ=π
1.5
0.4
1.0
/9 +5
β=1
α=5
0.2 Numerical Sol. Single Piece Sol. Piecewise Sol.
0.5 0.0 0.0
α = 10/9 β = 2/9 φ = 2π/9 0.2
0.4
s
0.6
0.8
9 +π/
φ=π
0.0 1.0
φ = 2π/9 α = 10/9
-0.5
0.0
β = 2/9
0.5
x(s)
1.0
Figure 6: Comparison of the numerical solution, the single-piece solution (M1 = 16 and M2 = M3 = 24 with s0 = 1/2) of Ref. [25], and the piecewise solution (M1i = 6 and M2i = M3i = 8 with s01 = 1/4, s2 = 1/2, and s02 = 3/4) for the beam’s (a) angular and (b) axial deflections.
(a)
(b) 3π φ+
2α2
dM < 0 ds
Γmax
0
φ
φ+π
φ
M1
Γavg =β2+2α2 cos(θ1-φ)
Γmin
dM > 0 ds
Γmax
f
π φ+
dθ ds = β dM < 0 ds
θ
dM > 0 ds
Γmax φ dM < 0 ds
Γmin
φ+ 2π
Γ
Γmin
Figure 7: Schematic plot of (a) Γ(α, β, φ, θ1 , θ), and (b) a deflected beam undergoing large deformation.
from 0 to its target value. The periodicity shown in Eq. (30) is schematically shown in Figure 7, and originates from the change of the sign of Mf , which is the moment resulting from the force f . Illustrated in Figure 8, the force f like the applied moment M1 − dM1 creates a counter clockwise moment; thus, the two are supportive. However, as the applied moment is incremented to M1 and the beam undergoes further deflection, there may be a point (x, y) on the beam’s axis through which the force’s line of action passes. Therefore, the moment resulting from the force about this point vanishes. Setting Mf in Eq. (1) to zero gives y − y1 = tan(φ), x − x1 which also follows from intuition. This phenomenon is different from that of inflection points, as the nonzero external moment M1 creates a nonzero curvature. For simplicity, we refer to this phenomenon as force’s moment stagnation. Past this stagnation, if the external moment is incremented further, say to M1 +dM1 , the applied force starts to create a clockwise moment, and thus, opposes the applied moment until another stagnation occurs. For the quasi-static loading 15
considered here, this phenomenon occurs on the clamped point, but can occur on intermediate points of the beam as well, if a sufficiently large moment is applied.
f
Mf < 0
M 1+ dM 1
f
M1
Mf = 0
φ
f Mf > 0
M 1 - dM 1 Figure 8: Quasi-static application of external moment, and the occurrence of the force’s moment stagnation.
When the force and applied tip moment are supportive, the large deflection equations are well behaved. However, when the two are in conflict, the aforementioned sign switch poses numerical difficulties. Figure 9 compares the numerical and sixth order piecewise solutions for a constant and fixed applied force as the applied moment is increased in small increments. As illustrated for α = 5 and φ = 2π/3, the numerical differential solver of Mathematica diverges as β passes 5.2, while the piecewise solutions are convergent. To validate the accuracy of the latter solutions for such cases, we substituted the piecewise solutions in Eqs. (3), (5), and (6), and monitored the forward error, which never exceeded 0.01 along the entire beam length.
(a) 0.8
α=5 φ = 2π/3
3
0.6
β=7
Numerical Sol. Piecewise Sol.
α=5 φ = 2π/3
0.4
2 +0.
θ(s), rad
4
Numerical Sol. Piecewise Sol.
y(s)
5
(b)
2
β=7
+0.2
0 0.0
0.2
β=0
1
β=0
0.0 0.2
0.4
0.6
0.8
1.0
-0.4
s
-0.2
0.0
0.2
0.4
x(s)
Figure 9: Comparison of the numerical solution and the sixth order piecewise solution (s01 = 1/4, s2 = 1/2, and s02 = 3/4) for the beam’s (a) angular and (b) axial deflections.
16
6.5. BeamSol: A Large Deflection Beam Solver To further demonstrate the effectiveness of our solution procedure, we coded our algorithm as a black-box solver in Visual C++. Our large deflection beam solver, BeamSol 1.0, provides various user-friendly functionalities, including plotting and sampling the solutions, as well as exporting them to a data file.
Case 8
Case 9
Case 3 Case 10 Case 1 Case 2 Case 5
Case 4 Case 7
Case 6
Figure 10: A snapshot of BeamSol 1.0 corresponding to the ten cases of input data given in Table 1. The annotations in red are added manually for illustration purposes.
For the special case of pure end-moment, BeamSol 1.0 employs the circular exact solution of Ref. [25], and for a general loading it employs, currently, the two-piece fifth order ATET solutions given by Eqs. (27)-(29), but is envisioned to offer various number of pieces and approximation orders in the very near future. For the solution of boundary constraints, as discussed earlier, the Newton-Raphson technique is implemented and a zero solution set is used. To validate the results of BeamSol 1.0, we ran it for many randomly generated input data, and the results compared well with the numerical differential solver of Mathematica. The executable image of BeamSol 1.0 will be posted online for free download. Table 1: Input data used for Figure 10. For all cases: l = 1 (m), I = 10−4 (m4 ), and E = 10 (MPa).
Parameter \ Case #
1
2
3
4
5
6
7
8
9
10
f (kN )
0
0
1
1
16
16
16
16
16
16
φ (radian)
–
–
2
-2
-1
-1
-1
1
1
1
M1 (kN.m)
5
-5
4
-4
4
0
-4
4
0
-4
Figure 10 depicts a snapshot of BeamSol 1.0, which plots the cartesian large deflection 17
components for ten cases that are listed in Table 1. 7. Summary and Conclusions We have presented cartesian large deflection solutions, which are independent of a beam’s arc length. The solutions are implicit and in terms of elliptic integrals, but are envisioned to guide the development of approximate closed form large deflection solutions. For the case of a large pure end-force, a new closed form solution was also presented which incorporates the rotation matrix employed for developing the aforementioned cartesian large deflection solutions. To add further efficiency and robustness to the single-piece parametric solutions, piecewise large deflection solutions have been developed which require low approximation orders for fast convergence. The proposed methodology can quantify the inherently nonlinear and complex load, deformation, and energy relationships revealed by compliant mechanisms efficiently and effectively, thus, paving the way for the emergence of novel compliant mechanisms syntheses. Approximate linear theories are still the premise of certain devices in precision engineering. However, the overall dimensions of such devices must be adjusted with the magnitudes of the applied loads. The utilization of the proposed solution methodology can help in the design of new devices with fixed size and no dependence on the load magnitudes. We coded the solutions as a stand-alone black-box solver, which offers several functionalities for the solution of large deflection of EulerBernoulli cantilever beams. The executable will be posted for free download for those who are interested in beam applications. Acknowledgments The authors are grateful to the anonymous reviewers for the mindful comments. [1] G. K. Bennig, Atomic force microscope and method for imaging surfaces with atomic resolution, U.S. Patent 4,724,318, 1988. [2] D. J. Leo, Engineering Analysis of Smart Material Systems, Wiley InterScience, Wiley, 2007. [3] L. L. Howell, Compliant Mechanisms, John Wiley & Sons, 2001. [4] R. H. Burns, The Kinetostatic Synthesis of Flexible Link Mechanisms, Ph.D. Dissertation, Yale University, 1964. [5] R. H. Burns, F. R. E. Crossley, Kinetostatic Synthesis of Flexible Link Mechanisms, ASME Mechanisms Conference, ASME Paper No 68–MECH–36, 1968. 18
[6] S. J. Winter, T. E. Shoups, The displacement analysis of path-generating flexible-link mechanisms, Mechanism and Machine Theory 7 (4) (1972) 443–451. [7] C. Huang, B. Roth, Dimensional Synthesis of Closed-Loop Linkages to Match Force and Position Specifications, ASME Journal of Mechanical Design 115 (2) (1993) 194–198. [8] V. Krovi, G. K. Ananthasuresh, V. Kumar, Kinematic and Kinetostatic Synthesis of Planar Coupled Serial Chain Mechanisms, Journal of Mechanical Design 124 (2) (2002) 301–312. [9] M. H. F. Dado, Limit position synthesis and analysis of compliant 4-bar mechanisms with specified energy levels using variable parametric pseudo-rigid-body model, Mechanism and Machine Theory 40 (8) (2005) 977–992. [10] H. Tari, H.-J. Su, A Complex Solution Framework for the Kinetostatic Synthesis of a Compliant Four-Bar Mechanism, Mechanism and Machine Theory 46 (8) (2011) 1137–1152. [11] H. Tari, H.-J. Su, J. D. Hauenstein, Classification and complete solution of the kinetostatics of a compliant Stewart–Gough platform, Mechanism and Machine Theory 49 (0) (2012) 177–186. [12] J. Argyris, R. Kaˇcianauskas, Semi-analytical finite elements in the higher-order theory of beams, Computer Methods in Applied Mechanics and Engineering 138 (1–4) (1996) 19–72. [13] J. Mazars, P. Kotronis, F. Ragueneau, G. Casaux, Using multifiber beams to account for shear and torsion: Applications to concrete structural elements, Computer Methods in Applied Mechanics and Engineering 195 (52) (2006) 7264–7281. [14] L. Jun, H. Hongxing, S. Rongying, Dynamic finite element method for generally laminated composite beams, International Journal of Mechanical Sciences 50 (3) (2008) 466–480. [15] T. P. Vo, H.-T. Thai, Vibration and buckling of composite beams using refined shear deformation theory, International Journal of Mechanical Sciences 62 (1) (2012) 67–76. [16] E. O. nate, A. Eijo, S. Oller, Simple and accurate two-noded beam element for composite laminated beams using a refined zigzag theory, Computer Methods in Applied Mechanics and Engineering 213–216 (0) (2012) 362–382. [17] A. J. Sadowski, J. M. Rotter, Solid or shell finite elements to model thick cylindrical tubes and shells under global bending, International Journal of Mechanical Sciences 74 (0) (2013) 143–153.
19
[18] F. Miao, G. Sun, K. Chen, Transient response analysis of balanced laminated composite beams by the method of reverberation-ray matrix, International Journal of Mechanical Sciences 77 (0) (2013) 121–129. [19] Z. Zhang, X. Huang, Z. Zhang, H. Hua, On the transverse vibration of Timoshenko doublebeam systems coupled with various discontinuities, International Journal of Mechanical Sciences 89 (0) (2014) 222–241. [20] C. W. Lim, M. Z. Islam, G. Zhang, A nonlocal finite element method for torsional statics and dynamics of circular nanostructures, International Journal of Mechanical Sciences 9495 (0) (2015) 232–243. [21] K. E. Bisshop, D. C. Drucker, Large Deflection of Cantilever Beams, Quarterly of Applied Mathematics 3 (1945) 272–275. [22] L. L. Howell, The Design and Analysis of Large-Deflection Members in Compliant Mechanisms, M.S. Thesis, Purdue University, 1991. [23] A. Banerjee, B. Bhattacharya, A. Mallik, Large deflection of cantilever beams with geometric non-linearity: Analytical and numerical approaches, International Journal of Non–Linear Mechanics 43 (5) (2008) 366–376. [24] L. Chen, An integral approach for large deflection cantilever beams, International Journal of Non–Linear Mechanics 45 (3) (2010) 301–305. [25] H. Tari, On the parametric large deflection study of Euler-Bernoulli cantilever beams subjected to combined tip point loading, International Journal of Non-Linear Mechanics 49 (2013) 90–99. [26] F. A. Chouery, Exact and Numerical Solutions for Large Deflection of Elastic Non-Prismatic Beams, FAC Systems INC., WA .
20