TipLoaded 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 tiploaded cantilever beam is modeled using PFC2D and PFC3D. The computed deformations and the maximum shear and normal stresses are compared with the analytical beamtheory solution. This problem demonstrates the capabilities of the parallelbond logic in PFC.
The physical model of the beam, with a length of 200 inches and a circular crosssection 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 2024T3 aluminum with a Young’s modulus, \(E\), of 1 × 10^{4} ksi and a shear modulus, \(G\), of 3759.4 ksi. The loading consists of transverse load, \(F\), axial load, \(P\), (inplane) bending moment, \(M\), and (outofplane) twisting moment, \(T\) (in 3D only), applied at the beam tip.
ClosedForm Solution
A closedform 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):
where:
The total contribution from the bending load and applied twisting moment to the maximum shear and normal stresses is given by
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
For the special case in which all of the springs have the same stiffness, \(k\), the expression for the equivalent stiffness simplifies to
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
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
Both the normal and shear strength of the parallel bonds are set to the large value of 1 × 10^{5} 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^{5} 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 closedform 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.
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 closedform solution in Table 1. The total rotation of the ball at the beam tip is monitored by activating orientation tracking (using the model orientationtracking
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.
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 parallelbond 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 closedform solution in Table 2. The total rotation of the ball at the beam tip is monitored by activating orientation tracking (using the model orientationtracking
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.
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% 
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 parallelbond 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.25e1.
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 largestrain on
model random
model domain ext 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] rad 10 positionx [_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 ydisplacement of ball11.
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 dens=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 velocityx velocityy velocityz ...
spinx spiny spinz range id=1
;  History variables
fish history _xdel11
fish history _ydel11
fish history _zdel11
fish history _xrot11
fish history _zrot11
model orientationtracking on
model save 'initial'
[apply_loads(1.0,100.0,2000.0,0.0)]
model solve ratiomaximum 1.0e4
[get_maxpbstr]
fish list [_xdel11] [_ydel11] [_xrot11] [_zrot11] ...
[_maxpbnstr] [_maxpbsstr] [_maxpbnx] [_maxpbsx]
model save 'finalcase1'
model restore 'initial'
[apply_loads(0.0,0.0,0.0,1000.0)]
model solve ratiomaximum 1.0e4
fish list [_xdel11]
model save 'finalcase2'
model restore 'initial'
contact method unbond
[apply_loads(0.0,0.0,0.0,1000.0)]
fish list [_xdel11]
model solve ratiomaximum 1.0e4
model save 'finalcase3'
model restore 'initial'
[apply_loads(0.0,0.0,0.0,1000.0)]
fish list [_xdel11]
model solve ratiomaximum 1.0e4
model save 'finalcase4'
;=========================================================================
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 largestrain on
model random
model domain ext 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] rad 10 pos([_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 ydisplacement 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 dens=1000 damp 0.7
ball prop '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 velocityx velocityy spin range id=1
;  History variables
fish history _xdel11
fish history _ydel11
fish history _rot11
model history mechanical age
model orientationtracking on
model save 'initial'
[apply_loads(1.0,100.0,0.0)]
model solve ratiomax 1e5
[get_maxpbstr]
fish list [_xdel11] [_ydel11] [_rot11] [_maxpbnstr] ...
[_maxpbsstr] [_maxpbnx] [_maxpbsx]
model save 'finalcase1'
model restore 'initial'
[apply_loads(0.0,0.0,1000.0)]
model solve ratiomax 1.0e5
model save 'finalcase2'
model restore 'initial'
contact method unbond
[apply_loads(0.0,0.0,1000.0)]
model solve ratiomax 1.0e5
model save 'finalcase3'
model restore 'initial'
[apply_loads(0.0,0.0,1000.0)]
model solve ratiomax 1.0e5
model save 'finalcase4'
; =========================================================================
program return
; eof: cantilever.dat (2D)
Endnotes
⇐ Strength of a FaceCentered Cubic Array of Spheres  Linear Contact Model: Calibrating the Normal Critical Damping Ratio ⇒
Was this helpful? ...  Itasca Software © 2023, Itasca  Updated: Jan 18, 2024 