Tip-Loaded Cantilever Beam

Note

The project file for this example may be viewed/run in PFC.[1] The data files used are shown at the end of this example.

Problem Statement

A tip-loaded cantilever beam is modeled using PFC2D and PFC3D. The computed deformations and the maximum shear and normal stresses are compared with the analytical beam-theory solution. This problem demonstrates the capabilities of the parallel-bond logic in PFC.

The physical model of the beam, with a length of 200 inches and a circular cross-section with radius, \(R\), of 8 inches in 3D, is shown in portion (a) of Fig. 1 below. In 2D, the beam has a rectangular section of unit thickness and height \(2R\). The beam material is 2024-T3 aluminum with a Young’s modulus, \(E\), of 1 × 10⁴ ksi and a shear modulus, \(G\), of 3759.4 ksi. The loading consists of transverse load, \(F\), axial load, \(P\), (in-plane) bending moment, \(M\), and (out-of-plane) twisting moment, \(T\) (in 3D only), applied at the beam tip.

Closed-Form Solution

A closed-form expression for the deformation and maximum shear and normal stresses of the cantilever beam is given by simple beam theory. The solution is expressed in terms of the Cartesian coordinate system shown in Fig. 1 below. The origin of this coordinate system is at the fixed end of the beam, and the \(x\)-axis is directed toward the beam tip and lies along the neutral axis.

Transverse load, \(F\), and bending moment, \(M\), produce a transverse deflection, \(Δ_z\), and a rotation about the \(y\)-axis, \(θ_y\). Axial load, \(P\), produces an axial extension, \(Δ_x\). The twisting moment, \(T\), produces a rotation, \(θ_y\), about the \(x\)-axis. These deformations are given by McGuire and Gallagher (1979):

(1)\[\begin{split}\Delta_y &= { Fx^2 \over 2EI } \left( L - {x\over3} \right) + {Mx^2 \over 2EI} \\ \\ \Delta_x &= { Px \over EA } \\ \\ \theta_z &= { Fx \over EI } \left( L - {x\over2} \right) + {Mx \over EI} \\ \\ \theta_x &= { Tx \over GJ } ~\mbox{(3D only)}\end{split}\]

where:

(2)\[\begin{split}A &= \begin{cases} \pi R^2 ~\mbox{, in 3D} \\ 2 R ~\mbox{, in 2D} \end{cases} \\ I &= \begin{cases} \textstyle{1\over4} \pi R^4 ~\mbox{, in 3D} \\ \textstyle{2\over3} \pi R^3 ~\mbox{, in 2D} \end{cases} \\ \\ \\ J &= 2I ~\mbox{ (3D only)}\end{split}\]

The total contribution from the bending load and applied twisting moment to the maximum shear and normal stresses is given by

(3)\[\begin{split}\sigma_{\rm max} &= { FLR \over I } \left( 1 - {x\over L} \right) + { MR \over I } \\ \\ \tau_{\rm max} &= { F \over A } + \underbrace{ TR \over J }_{\text{3D only}}\end{split}\]

The extension of the beam subjected to an axial load, \(P\), can be found using \(P/K\), where \(K\) is the equivalent stiffness of a single spring. If the beam is modeled by a collection of \(n\) springs, each with a stiffness of \(k^{(i)}\), \(K\) can be determined based on whether the collection of springs is acting in series or in parallel, and is given by

(4)\[\begin{split}K = \begin{cases} \left( \sum\limits_{i=1}^n {1 / k^{(i)}} \right)^{-1},& \mbox{springs acting in series} \\ \\ \sum\limits_{i=1}^n k^{(i)},& \mbox{springs acting in parallel} \end{cases}\end{split}\]

For the special case in which all of the springs have the same stiffness, \(k\), the expression for the equivalent stiffness simplifies to

(5)\[\begin{split}K = \begin{cases} k/n, &\mbox{springs acting in series} \\ \\ nk, &\mbox{springs acting in parallel } \end{cases}\end{split}\]

The analytical solution for the beam modeled using \(n\) parallel bonds and \(n\) contacts is found by noting that the parallel bonds act in parallel with the contact springs to carry the load, and that the equivalent stiffnesses of the parallel bonds and the ball contacts are given by Equation (5), such that the equivalent stiffness of the entire beam is given by

(6)\[K = \left( A \bar k_n \over n \right) + \left( k_n \over 2n \right)\]

where \(\bar k_n\) and \(k_n\) are the normal stiffnesses of each parallel bond and each ball, respectively. Note that the contact stiffness for two balls with identical stiffnesses is equal to half the ball stiffness.

Numerical Model

The cantilever beam is modeled in PFC by joining together a row of 11 balls of equal radius with parallel bonds, as shown in portion (b) of Fig. 1 below. The ball radii are chosen to be larger than the radius of the cantilever beam in order to demonstrate the relation between each parallel bond and an equivalent elastic cylinder. The radius of each ball is 10 inches, but the radius of each parallel bond is set to 8 inches to match the physical beam. The normal and shear stiffnesses of the parallel bonds are chosen such that the collection of ten equivalent elastic cylinders, shown in portion (c) of Fig. 1 below, mimics the physical beam. Thus, the chosen values are

(7)\[\begin{split}\bar k^n &\equiv {\bar E \over \bar L} = { 10000 \over 20 } = 500\hbox{ kips/inch}^3 \\ \\ \bar k^s &\equiv {\bar G \over \bar L} = { 3759.4 \over 20 } = 188\hbox{ kips/inch}^3\end{split}\]

Both the normal and shear strength of the parallel bonds are set to the large value of 1 × 10⁵ ksi, so that inertial effects arising from sudden application of the loads do not cause the parallel bonds to break. The normal and shear stiffnesses of each of the balls (2.0106 × 10⁵ kips/inch) are chosen by setting the second term of Equation (6) equal to the overall beam stiffness of \(AE/L\).

The velocities of all six degrees of freedom are fixed to zero for the ball at the origin of the coordinate system, and three of the generalized forces are applied to the ball at the opposite end of the assembly. The model is then cycled until static equilibrium is reached (i.e., the velocities of all balls become negligible). Since the closed-form solution based on simple beam theory is no longer valid when an axial load is applied to the deflected configuration, because the problem becomes geometrically nonlinear, the axial load is applied as a separate case.

../../../../../../../_images/cantilever.png — Figure 1: Cantilever beam and *PFC* model. (\(M\) is an in-plane bending moment. \(T\) is an out-of-plane twisting moment in the 3D model.)

Results and Discussion

PFC3D Model

The deformations at the tip of the beam arising from the application of transverse load, \(F\), of 1.0 kip; bending moment, \(M\), of 100.0 kip × inch; and twisting moment, \(T\), of 2000.0 kip × inch are compared with the closed-form solution in Table 1. The total rotation of the ball at the beam tip is monitored by activating orientation tracking (using the model orientation-tracking command). FISH variables _xdel11, _ydel11 and _xrot11, _zrot11 store accumulated displacements and rotations, respectively, at the beam tip. The computed values of maximum normal and shear stress within the parallel bonds are found by using the FISH function get_maxpbstr. The maximum normal stress is found to be acting in the parallel bond centered at \(x\) = 10. The value of 7.2117 × 10^-1 ksi agrees exactly (within five significant figures) with the analytical solution of 7.2117 × 10^-1 ksi, at \(x\) = 10. The maximum shear stress arises from the transverse load and the twisting moment and is uniform throughout the length of the beam. The value of 2.4918 ksi agrees exactly (within five significant figures) with the analytical solution. These results can be obtained by printing _maxpbnstr, _maxpbsstr, _maxpbnx, and _maxpbsx.

Table 1: Deformation Comparison (3D model) –
Apply F, M, and T
	Numerical	Analytical	Relative Error
\(\Delta_y(L)\)	1.4512 × 10^-1 in.	1.4506 × 10^-1 in.	< 0.1%
\(\theta_z(L)\)	1.2436 × 10^-3 rad.	1.4506 × 10^-3 rad.	< 0.1%
\(\theta_x(L)\)	1.6535 × 10^-2 rad.	1.6537 × 10^-2 rad.	< 0.1%

The results of three separate cases in which only the axial load, \(P\), was applied to the beam are presented here. The cases differ based on whether the load is carried by the parallel bonds, the contact springs, or both.

In the first case, \(P\) = +1000 kips, and the parallel bonds are left intact. This puts the beam in tension, and all of the load is carried by the parallel-bond springs. The computed axial extension of 9.9472 x 10^-2 agrees exactly (within five significant figures) with the analytical solution.

In the second case, \(P\) = -1000 kips, and the parallel bonds are removed. This puts the beam in compression, and the load is carried by the contact springs between the balls. The computed axial extension of -9.9473 × 10^-2 is within 0.1 percent of the analytical solution of -9.9472 × 10^-2.

In the third case, \(P\) = -1000 kips, and the parallel bonds are left intact. This puts the beam in compression, and the load is shared by both the parallel bond springs and the contact springs between the balls. This case does not correspond to the physical problem, since both the parallel bond and contact stiffnesses were chosen assuming that only one was acting at a time. The analytical solution for this case is found using Equation (6). The computed axial extension of -4.9736 × 10^-2 agrees exactly (within five significant figures) with the analytical solution.

PFC2D

The deformations at the tip of the beam arising from the application of transverse load, \(F\), of 1.0 kip, and bending moment, \(M\), of 100.0 kip × inch, are compared with the closed-form solution in Table 2. The total rotation of the ball at the beam tip is monitored by activating orientation tracking (using the model orientation-tracking command). FISH variables _xdel11, _ydel11, and _rot11 store accumulated displacements and rotations, respectively, at the beam tip. The computed values of maximum normal and shear stress within the parallel bonds is found by using the FISH function get_maxpbstr. The maximum normal stress is found to be acting in the parallel bond centered at \(x\) = 10. The value of 7.2117 × 10^-1 ksi agrees exactly (within five significant figures) with the analytical solution of 7.2117 × 10^-1 ksi, at \(x\) = 10. The maximum shear stress arises from the transverse load and the twisting moment and is uniform throughout the length of the beam. The value of 2.4918 ksi agrees exactly (within five significant figures) with the analytical solution. These results can be obtained by printing _maxpbnstr, _maxpbsstr, _maxpbnx, and _maxpbsx.

Table 2: Deformation Comparison (2D Model) -
Apply F and M
	Numerical	Analytical	Relative Error
\(\Delta_y(L)\)	1.3685 in.	1.3672 in.	< 0.1%
\(\theta_z(L)\)	1.1718 × 10^-2 rad.	1.1719 × 10^-2 rad.	< 0.1%

In the first case, \(P\) = +1000 kips, and the parallel bonds are left intact. This puts the beam in tension, and all of the load is carried by the parallel-bond springs. The computed axial extension of 1.2499 agrees exactly, within four significant figures, with the analytical solution of 1.25.

In the second case, \(P\) = -1000 kips, and the parallel bonds are removed. This puts the beam in compression, and the load is carried by the contact springs between the balls. The computed axial extension of -1.2499 agrees exactly, within four significant figures, with the analytical solution of -1.25.

In the third case, \(P\) = -1000 kips, and the parallel bonds are left intact. This puts the beam in compression, and the load is shared by both the parallel bond springs and the contact springs between the balls. This case does not correspond to the physical problem, since both the parallel bond and contact stiffnesses were chosen assuming that only one was acting at a time. The analytical solution for this case is found using Equation (6). The computed axial extension of -6.2497 × 10^-1 agrees exactly (within four significant figures) with the analytical solution of -6.25e-1.

References

McGuire, W., and R. H. Gallagher. Matrix Structural Analysis, p. 87. New York: John Wiley & Sons, 1979.

Data Files

cantilever.dat (3D)

; fname: cantilever.dat (3D)
;
; Modeling of cantilever beam using a single row of balls
; joined by parallel bonds.
;
;=========================================================================
model new
model large-strain on
model random
model domain extent -50 250 -50 50
contact cmat default model linearpbond

fish define make_cantilever
  global _x = 0.0
  global _ballcnt = 1
  loop _ballcnt (1,11)
    command
      ball create id [_ballcnt] radius 10 position-x [_x] 
    end_command
    _x = _x + 20.0
  end_loop
  global bp_11 = ball.find(11)
  _xrot11 = 0
  _zrot11 = 0
end
; ----- Apply load(s) to ball.11 at tip
fish define apply_loads(_f,_m,_t,_p)
  ball.force.app.y(bp_11)  = _f
  ball.moment.app.z(bp_11) = -1.0*_m
  ball.moment.app.x(bp_11) = _t
  ball.force.app.x(bp_11)  = _p
end
; ----- Find total x- and y-displacement of ball-11.
fish define _xdel11
  _xdel11 = ball.pos(bp_11,1) - 200.0
  _ydel11 = ball.pos(bp_11,2)
  _zdel11 = ball.pos(bp_11,3)
  _xrot11 = ball.euler(bp_11,1)*math.degrad
  _zrot11 = ball.euler(bp_11,3)*math.degrad
end
; ----- Find maximum p.b. stresses and their locations
; out of all p.b. in the model
; OUTPUT VARS:
; _maxpbnstr : max nstress
; _maxpbsstr : max sstress
; _maxpbnx : x coord. of p.b. with max nstress
; _maxpbsx : x coord. of p.b. with max sstress
;
fish define get_maxpbstr
  global _maxpbnstr = -1.0
  global _maxpbsstr = -1.0
  global _maxpbnx = -1.0
  global _maxpbsx = -1.0
  loop foreach local cp contact.list
    if contact.prop(cp,"pb_state") = 3 then
      local rad = ...
         math.min(ball.radius(contact.end1(cp)), ...
                  ball.radius(contact.end2(cp)))
      local br = contact.prop(cp,"pb_rmul") * rad
      local barea = math.pi*br*br           ; treat as spheres  
      local bi = 0.25*barea*br*br           ; treat as spheres
      local bmfac = contact.prop(cp,"pb_mcf")
      local bf    = contact.prop(cp,"pb_force")
      local bfs   = math.sqrt(comp.y(bf)^2+comp.z(bf)^2)
      local bfm   = contact.prop(cp,"pb_moment")
      local tfm   = math.abs(comp.x(bfm)) 
      bfm         = math.sqrt(comp.y(bfm)^2+comp.z(bfm)^2)
      
      local nsmax = - comp.x(bf) / barea  + bmfac * bfm *br/bi
      local ssmax =   bfs / barea + bmfac * tfm * br/(2.0*bi)
      if nsmax > _maxpbnstr then
        _maxpbnstr = nsmax
        _maxpbnx = contact.pos(cp,1)
      end_if
      if ssmax > _maxpbsstr then
        _maxpbsstr = ssmax
        _maxpbsx = contact.pos(cp,1)
      end_if
    end_if
  end_loop
end
; ----- Begin execution here
[make_cantilever]
; ----- Set model properties
; (Note that pb_n and pb_s are set to extremely large
; values so that inertial effects arising from the
; instantaneous application of the final loads do
; not cause the parallel bonds to break.)
ball attribute density 1000 damp 0.7
ball property 'kn'=2.0106e5 'ks'=2.0106e5
model clean
; ----- Add parallel bonds
contact method bond 
contact property pb_ten 1e5 pb_coh 1e5 pb_kn 500 pb_ks 188 pb_rmul 0.8
ball fix velocity-x velocity-y velocity-z ...
         spin-x spin-y spin-z range id=1

; ----- History variables
fish history _xdel11
fish history _ydel11
fish history _zdel11
fish history _xrot11
fish history _zrot11
model orientation-tracking on
model save 'initial'

[apply_loads(1.0,-100.0,2000.0,0.0)]
model solve ratio-maximum 1.0e-4
[get_maxpbstr]
fish list [_xdel11] [_ydel11] [_xrot11] [_zrot11] ...
     [_maxpbnstr] [_maxpbsstr] [_maxpbnx] [_maxpbsx]
model save 'final-case1'


model restore 'initial'
[apply_loads(0.0,0.0,0.0,1000.0)]
model solve ratio-maximum 1.0e-4
fish list [_xdel11]
model save 'final-case2'

model restore 'initial'
contact method unbond 
[apply_loads(0.0,0.0,0.0,-1000.0)]
fish list [_xdel11]
model solve ratio-maximum 1.0e-4
model save 'final-case3'

model restore 'initial'
[apply_loads(0.0,0.0,0.0,-1000.0)]
fish list [_xdel11]
model solve ratio-maximum 1.0e-4
model save 'final-case4'
;=========================================================================
program return
; eof: cantilever.dat (3D)

cantilever.dat (2D)

; fname: cantilever.dat (2D)
;
; Modeling of cantilever beam using a single row of balls
; joined by parallel bonds.
;
; =========================================================================
model new
model large-strain on
model random

model domain extent -50 250 -50 50
contact cmat default model LinearPBond

fish define make_cantilever
  global _x = 0.0
  global _ballcnt = 1
  loop _ballcnt (1,11)
    command
      ball create id [_ballcnt] radius 10 position ([_x],[0]) 
    end_command
    _x = _x + 20.0
  end_loop
  global bp_11 = ball.find(11)
end

; ----- Apply load(s) to ball 11 (at tip)
fish define apply_loads(_f,_m,_p)
  ball.force.app.y(bp_11) = _f
  ball.moment.app(bp_11)  = _m 
  ball.force.app.x(bp_11) = _p
end
; ----- Find total x- and y-displacement of ball 11.
fish define _xdel11
  _xdel11  = ball.pos(bp_11,1) - 200.0
  _ydel11  = ball.pos(bp_11,2)
  _rot11   = ball.rotation(bp_11)*math.degrad
end
; ----- Find maximum p.b. stresses and their locations
; out of all p.b. in the model
; OUTPUT VARS:
; _maxpbnstr : max nstress
; _maxpbsstr : max sstress
; _maxpbnx : x coord. of p.b. with max nstress
; _maxpbsx : x coord. of p.b. with max sstress
;
fish define get_maxpbstr
  global _maxpbnstr = -1.0
  global _maxpbsstr = -1.0
  global _maxpbnx = -1.0
  global _maxpbsx = -1.0
  loop foreach local cp contact.list
    if contact.prop(cp,"pb_state") = 3 then
      local rad = math.min(ball.radius(contact.end1(cp)), ...
                           ball.radius(contact.end2(cp)))
      local br = contact.prop(cp,"pb_rmul") * rad
      local barea = math.pi*br*br           ; treat as spheres  
      local bi = 0.25*barea*br*br           ; treat as spheres
      local bmfac = contact.prop(cp,"pb_mcf")
      local bf    = contact.prop(cp,"pb_force")
      local bfm   = contact.prop(cp,"pb_moment")
      local nsmax = - comp.x(bf) / barea  + bmfac * math.abs(bfm) *br/bi
      local ssmax =   math.abs(comp.y(bf)) / barea
      if nsmax > _maxpbnstr then
        _maxpbnstr = nsmax
        _maxpbnx = contact.pos(cp,1)
      end_if
      if ssmax > _maxpbsstr then
        _maxpbsstr = ssmax
        _maxpbsx = contact.pos(cp,1)
      end_if
    end_if
  end_loop
end
; ----- Begin execution here
[make_cantilever]
; ----- Set model properties
; (Note that pb_n and pb_s are set to extremely large
; values so that inertial effects arising from the
; instantaneous application of the final loads do
; not cause the parallel bonds to break.)
ball attribute density=1000 damp 0.7
ball property 'kn'=16000 'ks'=16000
model clean

; ----- Add parallel bonds
contact method bond
contact property pb_ten 1e5 pb_coh 1e5 pb_kn 500 pb_ks 188 pb_rmul 0.8
ball fix velocity-x velocity-y spin range id=1

; ----- History variables
fish history _xdel11
fish history _ydel11
fish history _rot11
model history mechanical age
model orientation-tracking on

model save 'initial'

[apply_loads(1.0,100.0,0.0)]
model solve ratio-maximum 1e-5
[get_maxpbstr]
fish list [_xdel11] [_ydel11] [_rot11] [_maxpbnstr] ...
          [_maxpbsstr] [_maxpbnx] [_maxpbsx]
model save 'final-case1'


model restore 'initial'
[apply_loads(0.0,0.0,1000.0)]
model solve ratio-maximum 1.0e-5
model save 'final-case2'

model restore 'initial'
contact method unbond
[apply_loads(0.0,0.0,-1000.0)]
model solve ratio-maximum 1.0e-5
model save 'final-case3'

model restore 'initial'
[apply_loads(0.0,0.0,-1000.0)]
model solve ratio-maximum 1.0e-5
model save 'final-case4'

; =========================================================================
program return
; eof: cantilever.dat (2D)

Endnote

⇄

⇐ Strength of a Face-Centered Cubic Array of Spheres | Linear Contact Model: Calibrating the Normal Critical Damping Ratio ⇒

Was this helpful? ...

Updated: Oct 31, 2024