|
>
|
restart:with(linalg):with(plots):
|
Equilibrium Configurations of Cantilever Beam under Terminal Load
Milan Batista, University of Ljubljana, Slovenia, EU
Introduction
This worksheet implements the calculation of the equilibrium shapes of initially straight inextensible and unshearable elastic cantilever beam subject to terminal force. The parameter of the modela are:
|
•
|
load parameter omega ( w nondimensional)
|
|
•
|
the CW angle beta ( b in degrees) between the x-axis and the direction of force
|
|
•
|
the CW angle alpha (a in degrees) between tangent to the beam base curve at the tip point and direction of the terminal force
|
Three types of load conditions are implemented:
|
•
|
Follower load - the load parameter w and angle a are given. This problem has one solution.
|
|
•
|
Load parameter problem - both a and b are given. This problem has infinitely many solutions.
|
|
•
|
Conservative load - load parameter w and angle b are given. This problem has finite number of solutions.
|
Basic procedures
There are five basic procedures that implement the parametric equations of the beam base curve which are expressed in terms of Jacobi's elliptical functions. First two procedures calculate the coordinates of beam base curve in local coordinate system which is rotated CW with respect to space coordinate system by the angle b
|
>
|
#========================================
|
|
>
|
xi := proc( s, omega) global k, K, K1, E; if omega = 0 then (1 - s)*(1 - 2*k**2) else (2*E*K1-1)*(1 - s) + 2/omega*(JacobiZeta(omega,k) - JacobiZeta( omega*s, k) - k**2*(JacobiSN(omega,k)*JacobiCN(omega,k)/JacobiDN(omega,k)-JacobiSN(omega*s,k)*JacobiCN(omega*s,k)/JacobiDN(omega*s,k))) end if; end proc;
|
|
| (1) |
|
>
|
#========================================
|
|
>
|
eta := proc( s, omega) global k, k1, K; if omega = 0 then (1 - s)*2*k*k1 else 2*k*k1/omega*(JacobiSN(omega,k)/JacobiDN(omega,k) - JacobiSN(omega*s,k)/JacobiDN(omega*s,k)) end if; end proc;
|
|
| (2) |
|
>
|
#========================================
|
Next two procedures calculate the coordinates of beam base curve in space coordinates:
|
>
|
#========================================
|
|
>
|
x:=proc( s, omega) local xx, yy, beta0; global k, K; beta0:=beta(omega); xx:=xi( s, omega); yy:= eta(s, omega);xx*cos(beta0)+yy*sin(beta0);end proc;
|
|
| (3) |
|
>
|
#========================================
|
|
>
|
y:=proc( s, omega) local xx, yy, beta0; global k, K;beta0:=beta(omega); xx:=xi( s, omega); yy:= eta(s, omega);-xx*sin(beta0)+yy*cos(beta0);end proc;
|
|
| (4) |
|
>
|
#==============================================
|
Finaly, calculation of angle b is implemented by the following procedure
|
>
|
#==============================================
|
|
>
|
beta:=proc( omega) global k; 2*arcsin(k*JacobiCN(omega,k)/JacobiDN(omega,k)) end proc;
|
|
| (5) |
|
>
|
#==============================================
|
Global variables k , K and E in above procedures represent modulus of elliptic functions and complete elliptic integral of first kind and second kind respectively. These must be for computational efficiency calculated before use of the above procedures. When a is known then we can for setting k and K use procedure
|
>
|
#===================================================
|
|
>
|
setup := proc( alpha) global k, k1, K, K1, E; k:=evalf(sin(rad(alpha)/2)); k1:= sqrt(1 - k**2); K:=EllipticK(k); E:=EllipticE(k); if k < 1 then K1 := 1/K else K1 := 0 end if; end proc;
|
|
| (6) |
|
>
|
#====================================================
|
For displaying the applied force we need the following procedure
|
>
|
#===================================================
|
|
>
|
force:=(f,beta0)-><-f*cos(beta0),f*sin(beta0)>;
|
|
| (7) |
|
>
|
#====================================================
|
|
>
|
#=======================
|
|
>
|
rad:=alpha->evalf(alpha*Pi/180): # convert deg to rad
|
|
>
|
deg:=alpha->evalf(alpha*180/Pi): # convert rad to deg
|
|
>
|
#====================================
|
The variables that controls plots are:
|
•
|
number of points that is used for approximation of deformed beam
|
|
>
|
NumPts:=100; # number of points for beam shape
|
|
| (8) |
|
•
|
number of points used for approximation of beam tip path
|
|
| (9) |
|
•
|
the scale factor for the displaying the applied force
|
|
>
|
ForceScale:=0.025; # Scale for force
|
|
| (10) |
|
•
|
number of frames used for animations
|
|
>
|
NumFrames:=50; # number of frames for animations
|
|
| (11) |
Note that for larger value of load parameters the number of digits must be increased
|
| (12) |
Follower load
We first give four examples of calculation of beam equilibrium shape in the case of follower load. For this we need two data items. the angle a
|
>
|
#========================================
|
|
| (13) |
|
>
|
setup(alpha0): # Now K is availabe
|
and the load parameter
|
| (14) |
Note that K is globaly available and can thus be used for specifying load parameter.
|
>
|
#ForceScale:=0.025: # scale for force
|
|
>
|
#NumPts:=100: # number of points
|
This figure show the shape of deformed beam in local coordinate system
|
>
|
beta0:=beta(omega0): # we need this to show beam root
|
|
>
|
Beam:=plot([xi(s,omega0),eta(s,omega0),s=0..1],numpoints = NumPts, color=black, thickness=2, scaling=constrained,caption=typeset("omega = ", omega0," alpha = ",alpha0)): # Local coordinates
|
|
>
|
ForceVect:= arrow(<xi(0,omega0),eta(0,omega0)>,<-ForceScale,0>):
|
|
>
|
BeamRoot:=pointplot( [[ForceScale*sin(beta0),-ForceScale*cos(beta0)],[-ForceScale*sin(beta0),ForceScale*cos(beta0)]],connect=true):
|
|
>
|
display( Beam, ForceVect,BeamRoot);
|
|
>
|
#-------------------------------------
|
and this figure the shape of deformed beam in global coordinate system
|
>
|
#-------------------------
|
|
>
|
ForceVect:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,beta0)):
|
|
>
|
Beam:=plot([x(s,omega0),y( s , omega0),s=0..1],numpoints = NumPts, color=black, thickness=2, scaling=constrained,caption=typeset("omega = ", omega0," alpha = ",alpha0," beta = ",evalf(beta0*180/Pi))): # Space coordinates
|
|
>
|
display( Beam, ForceVect);
|
|
>
|
#========================================
|
|
>
|
# Follower load problem
|
|
>
|
#========================================
|
This animation shows continiuous deformation of beam when load factor is constant and a increase from 0 to 180 deg. The dots shows trace of beam tip.
|
>
|
# Constant load parameter -------------------------------------------------
|
|
>
|
xx0:=array(1..NumPtsPath):yy0:=array(1..NumPtsPath):
|
|
>
|
for i from 0 to NumPtsPath-1 do; alpha:=i/NumPtsPath; setup(alpha*180); xx0[i+1]:=x(0,omega0);yy0[i+1]:=y(0,omega0);od:
|
|
>
|
TipPath:=pointplot( {seq( [xx0[i],yy0[i]],i=1..NumPtsPath)}, color=red , symbolsize=10):
|
|
>
|
BaseCurve:=array(1..NumFrames):ForceVect:=array(1..NumFrames):
|
|
>
|
for i from 0 to NumFrames - 1 do alpha:=i/NumPts; setup(alpha*180);beta0:=beta(omega0);BaseCurve[i+1]:= plot([x(s,omega0),y(s,omega0), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=3,caption=typeset("frame #",i,"/",NumFrames," omega = ", omega0," alpha = ",alpha*180.," beta = ",evalf(beta0*180/Pi))); ForceVect[i+1]:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,beta0)):od:
|
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumFrames),insequence = true, scaling=constrained):
|
|
>
|
Force :=display(seq(ForceVect[i],i=1..NumFrames),insequence = true, scaling=constrained):
|
|
>
|
display( Beam, Force,TipPath);
|
Load parameter problem
There are two variants of the problem. In first we as data have the angle b:
|
>
|
#========================================
|
|
>
|
# Load parameter problem
|
|
>
|
#========================================
|
|
>
|
# These functions calculate load parameter as function of alpha and wave number n
|
|
>
|
#-------------------------------------------------------------------------
|
|
>
|
f1:=(alpha,n)->(-InverseJacobiSN(sin(beta0*Pi/360)/sin(alpha*Pi/2),sin(alpha*Pi/2)))+(4*n+1)*EllipticK(sin(alpha*Pi/2)):
|
|
>
|
f2:=(alpha,n)->-(-InverseJacobiSN(sin(beta0*Pi/360)/sin(alpha*Pi/2),sin(alpha*Pi/2)))+(4*n+3)*EllipticK(sin(alpha*Pi/2)):
|
|
>
|
#---------------------------------------------------------
|
|
| (15) |
|
>
|
if beta0 < 0 then beta0 := -beta0 end if:
|
and the signed wave number
|
>
|
m:=2; # Input wave number
|
|
| (16) |
The load parameter that correspond to these data items is represented by branch curve on the bifurcation diagram on (alpha,omega) plane. In the example, as shown in the graph below, the number of load parameters equals the number of frames for animation. These values are shown as circles on the corresponding branch curve.
|
>
|
if beta0 < 0 then beta0 := -beta0 end if;p:=evalf(sin(rad(beta0)/2)): # Check data
|
|
>
|
ma:=abs(m):if modp(abs(m),2) = 1 then n:=floor((ma-1)/2);else n:=floor(ma/2-1);end if:da:=(0.999-beta0/180)/(NumFrames-1):
|
|
>
|
for i from 0 to NumFrames - 1 do if m > 0 then alpha:=evalf(beta0/180+da*i) else alpha := evalf(-beta0/180 - da*i); end if; setup(evalf(alpha*180));omega0:=Re(-InverseJacobiSN(p/k,k));if modp(n,2) = 1 then w[i+1]:=omega0+(4*n+1)*K else w[i+1]:=-omega0 + (4*n + 3)*K end if; aa[i+1]:=alpha;od:
|
|
>
|
ss:=seq([aa[i],w[i]],i=1..NumFrames):
|
|
>
|
ww:=pointplot( {ss},color=black, symbol=circle,symbolsize = 10 ):
|
|
>
|
unassign('alpha');tt:=plot({seq(f1(alpha,n),n=0..2),seq(f2(alpha,n),n=0..2)},alpha=-1..1,y=0..mx,color=red,labels=["alpha/pi","omega"]):
|
Equilibrium shapes belonging to the above load parameter values areshown in the following animation
|
>
|
ForceScale:=0.2: # force scale factor
|
|
>
|
for i from 0 to NumFrames - 1 do alpha:=aa[i+1]; setup(evalf(alpha*180));omega:=w[i+1]; BaseCurve[i+1]:= plot([x(s,omega),y(s,omega), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=3,caption=typeset("frame #",i+1,"/",NumFrames," m = ",m," omega = ", omega," alpha = ",alpha*180.," beta = ",evalf(beta0))); ForceVect[i+1]:= arrow(<x(0,omega),y(0,omega)>,force(ForceScale*omega/w[NumFrames],evalf(beta0*Pi/180))):od:
|
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumFrames),insequence = true, scaling=constrained):
|
|
>
|
Force :=display(seq(ForceVect[i],i=1..NumFrames),insequence = true, scaling=constrained):
|
This example assume that both angles are given
|
| (17) |
|
>
|
alpha0:=145; # Input alpha
|
|
| (18) |
and the wave numbers increase from 1 to
|
>
|
NumWaves:=6; # How many waves
|
|
| (19) |
|
>
|
ForceScale:=0.15: # force scale factor
|
|
>
|
if beta0 < 0 then beta0 := -beta0 end if:p:=evalf(sin(beta0*Pi/180/2)):setup(alpha0): # check beta
|
In this graph the corresponding load parameter values are indicated by circles
|
>
|
for i from 1 to NumWaves do omega0:=Re(K-InverseJacobiSN(p/k,k));if modp(abs(i),2) = 1 then n:=(abs(i)-1)/2; w[i]:=omega0+4*n*K else n:=abs(i)/2-1; w[i]:=-omega0 + 4*(n + 1)*K end if;od:
|
|
>
|
ss:=seq([alpha0/180.,w[i]],i=1..NumWaves):
|
|
>
|
ww:=pointplot( {ss},color=black, symbol=solidcircle,symbolsize = 10,labels=["alpha/pi","omega"] ):
|
|
>
|
unassign('alpha');nn:=floor(NumWaves/2)+1:graph:=plot({seq(f1(alpha,n),n=0..nn),seq(f2(alpha,n),n=0..nn)},alpha=-1..1,y=0..w[NumWaves],color=red,labels=["alpha/pi","omega"]):
|
|
>
|
alp:=pointplot( [[alpha0/180,0],[alpha0/180,w[NumWaves]]],connect=true,color=black,thickness=1,labels=["alpha/pi","omega"]):
|
|
>
|
display(graph,ww, alp);
|
and this animation shows coresponding equilibrium configurations (for better view decrease the speed of the animation)
|
>
|
for i from 1 to NumWaves do omega:=w[i]; BaseCurve[i]:= plot([x(s,omega),y(s,omega), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=2,caption=typeset("wave = ",i," omega = ", omega," alpha = ",alpha0," beta = ",evalf(beta0))); ForceVect[i]:= arrow(<x(0,omega),y(0,omega)>,force(ForceScale*omega/w[NumWaves],evalf(beta0*Pi/180))):od:
|
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumWaves),insequence = true, scaling=constrained):
|
|
>
|
Force :=display(seq(ForceVect[i],i=1..NumWaves),insequence = true, scaling=constrained):
|
This picture display all equilibrium configurations
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumWaves), scaling=constrained):
|
|
>
|
Force :=display(seq(ForceVect[i],i=1..NumWaves), scaling=constrained):
|
Conservative Load
In this kind of problem the load parameter and the space direction of force are constant. The load parameter for this example is
|
>
|
#========================================
|
|
>
|
# (3) Conservative Load Problem
|
|
>
|
#========================================
|
|
| (20) |
|
| (21) |
|
>
|
#---------------------------------------------
|
|
>
|
# This is variant of procedure beta where also k and K are calculates
|
|
>
|
f:=proc(alpha,omega) global k, K; k:=evalf(sin(alpha/2));K:=EllipticK(k); beta(omega); end proc:
|
|
>
|
#------------------------------------------
|
|
>
|
# This is first derivative of function f which is needed to calculate extrem values
|
|
>
|
df := (k,omega)->2*(JacobiSN(omega,k)*omega*EllipticK(k)*k^2-JacobiSN(omega,k)*omega*EllipticK(k)+JacobiCN(omega,k)*JacobiDN(omega,k)*EllipticK(k)+JacobiSN(omega,k)*JacobiZeta(omega,k)*EllipticK(k)+JacobiSN(omega,k)*omega*EllipticE(k))/JacobiDN(omega,k)^2/EllipticK(k)/((-1+k^2)/(-1+k^2*JacobiSN(omega,k)^2))^(1/2):
|
|
>
|
#---------------------------------------------
|
|
>
|
# These two functions enclose n-th root of function f
|
|
>
|
f1:=n->2*arcsin(sqrt(1-exp(Pi-omega0*2/(2*n-1)))):f2:=n->2*arcsin(sqrt(1-16*exp(-omega0*2/(2*n-1)))):
|
|
>
|
#----------------------------------
|
|
>
|
# Calculate roots of f in the interval -Pi..Pi
|
|
>
|
#--------------------------------------------
|
For calculation of zeros the number of digits is increased becouse the intervals near the end of domain of force inclination angle becomes very small
|
>
|
dg:=Digits: Digits:=22; # Increase number of digits
|
|
>
|
N:=floor(omega0/Pi+1/2): # half number of possible nontrivial eqilibrium configurations
|
|
| (22) |
The maximum number of possible equilibrium configurations for this case is
|
>
|
M:=2*N+1; # max. number of equilibrium configurations
|
|
| (23) |
|
>
|
z0:=array(1..M): # initialize array
|
|
>
|
unassign('alpha'):for n from 1 to N do alpha1:=evalf(f1(n)); alpha2:=evalf(f2(n));z0[n]:=evalf(fsolve(f(alpha,omega0),alpha=alpha1..alpha2)/Pi);od:
|
|
>
|
# obtain rest of roots by symetry
|
|
>
|
z0[N+1]:=0:for n from 1 to N do z0[N+1+n]:=-z0[N+1-n] od:
|
|
>
|
#------------------------
|
|
>
|
zero0:=(seq( pointplot( [z0[i],0], symbol = circle, symbolsize = 10),i=1..M)):
|
|
>
|
graph:=plot({evalf(f(alpha*Pi,omega0)/Pi)},alpha=-1..1, color=black, numpoints=400): # note that number of points should be relatively large
|
|
>
|
#-------------------------------------
|
|
>
|
# Calculate extems of function f
|
|
>
|
#------------------------------------
|
|
>
|
s:=array(1..M+1):r:=array(1..M+1):
|
|
>
|
s[1]:=1:r[1]:=evalf(Pi):
|
|
>
|
for n from 1 to M-1 do;if modp(n,2) =1 then xx:=fsolve(df(evalf(sin(alpha*Pi/2)),omega0),alpha=z0[n+1]..z0[n]); else xx:=fsolve(df(evalf(sin(alpha*Pi/2)),omega0),alpha=z0[n+1]..z0[n]); end if; s[n+1]:=xx;r[n+1]:=f(xx*Pi,omega0) od:s[M+1]:=-1:r[M+1]:=evalf(-Pi):
|
|
>
|
#for i from 1 to M do z0[i],r[i]; od;
|
|
>
|
ext0:=(seq( pointplot( [s[i],evalf(r[i]/Pi)], symbol = circle, symbolsize = 10, color=blue),i=1..M+1)):
|
|
>
|
#display(graph,zero0,ext0);
|
The inclination angle of the fors iss
|
>
|
beta0:=0; # Input <---------
|
|
| (24) |
|
>
|
#unassign('alpha'):fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=z0[1]..r[1]);z0[1];r[1];
|
|
>
|
if beta0 < 0 then beta0 := -beta0 end if; # Check
|
|
>
|
#-------------------------------------- Here the calulation phase is done !
|
|
>
|
if (beta0 > 0) then n:=0;for i from 1 to M do zz:=fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=z0[i]..s[i]): if type(zz,numeric) then n:=n+1;zzz[n]:=zz end if; zz:=fsolve(f(evalf(alpha*Pi),omega0)-rad(beta0),alpha=s[i+1]..z0[i]): if type(zz,numeric) then n:=n+1;zzz[n]:=zz end if; od: else n:=M; for i from 1 to M do zzz[i]:=z0[i] od; end if:
|
The actual number of equilibrium configurations is
|
| (25) |
and the angle a (in deg) for each of these configurations is
|
>
|
for i from 1 to NumShape do i,180*zzz[i] od;
|
|
| (26) |
These values are displayed as red circles on the following graph (the blue cicles are extrems of function)
|
>
|
zero:=(seq( pointplot( [zzz[i],evalf(beta0/180)], symbol = solidcircle, symbolsize = 12, color=red),i=1..n)):
|
|
>
|
display(graph,zero0,ext0,zero,plot(beta0/180,t=-1..1),labels=["alpha/pi","beta/Pi"]);
|
The coresponding equilibrium configurations are shown on the following animation ((for better view decrease the speed of the animation)
|
>
|
NumPts:=50:for i from 1 to NumShape do alpha:=zzz[i]; setup(alpha*180); BaseCurve[i]:= plot([x(s,omega0),y(s,omega0), s=0..1],scaling=constrained,numpoints = NumPts, color=black, thickness=2,caption=typeset("shape = ",i,"/",NumShape," omega = ", omega0," alpha = ",alpha*180.," beta = ",evalf(beta0))); ForceVect[i]:= arrow(<x(0,omega0),y(0,omega0)>,force(ForceScale,evalf(beta0*Pi/180))):od:
|
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumShape), scaling=constrained, insequence = true):
|
|
>
|
Force :=display(seq(ForceVect[i],i=1..NumShape), scaling=constrained, insequence = true):
|
On this picture all equilibrium configurations are shown
|
>
|
Beam := display(seq(BaseCurve[i],i=1..NumShape), scaling=constrained):
|
|
>
|
#Force :=display(seq(ForceVect[i],i=1..NumShape), scaling=constrained):
|
|
>
|
Digits:=dg: # reset number of digits
|
|
