Torsional Elasticity (Twisting)
Consider a beam extended in z. If we fix the z = 0 surface and apply a moment in the x- or y- direction on the z = Lz surface this causes the beam to bend in the y- or -x-direction, respectively. If instead we apply a moment in the z-direction on this face the beam will twist about the z-axis:
Torsion is deformation resulting from torque (picture from https://en.wikipedia.org/wiki/Torsion_(mechanics))
Copyright By PowCoder代写 加微信 powcoder
While bending moments are applied by axial stress gradients, twisting moments are applied by shear stress gradients: to produce a moment in z on the z=Lz face we need a changing shear stress profile on this surface:
Torsion forms the last Hooke’s law example, completing our study of elasticity:
Hooke’s Law
Hooke’s Law
(macroscopic)
Normal (axial; stretching/compressing)
Shear (sliding)
Flexion (bending)
Torsion (twisting)
Torsional stress and strain
A twisting moment applied on the z = Lz surface will cause the beam to experience a torsional strain: increasing twist angle per step in z: , where G is the shear modulus and is the torsion constant of the surface. Collectively, is called the torsional rigidity.
Note that if the twist angle is uniform at the cross section then the shear stress and strain increase linearly with distance from the rotation axis, reaching max values at the outer surface: , meaning the shear stress magnitude is .
If you assume no warping (i.e., a plane cross-section before twisting remains plane after twisting and a straight line across the face remains a straight line) then torsion constant is the second moment of area , where r is the radial distance from the centroid of the surface. However, this turns out to only be true for a circle; other cross sections warp in torsion and we need to use numerical methods (i.e., FEM like FlexPDE) or experiments to find torsion constants.
Cross Section
Torsion Constant
Where r2 is the outer radius and r1 is the inner radius
where is the major radius, and is the minor radius
where is half the side length
where & b are half the long & short side lengths, and
Alternatively, you can use this equation which is supposed to be accurate to within 4% error
(modified from https://en.wikipedia.org/wiki/Torsion_constant )
e.g., find the second moment of area of a square and circular cross section and compare with torsion constant
For the circle we can work in polar coordinates and say:
> restart:
Jcircle:=2*Pi*int(r^2*r, r=0..R);
which is exactly the torsion constant.
For the square, the tricky part is that the bounds in x & y are easy to set up but the integral itself is nasty. Fortunately the centroid is easy. Setting the origin at the centroid the cross section goes from x = -a to a and y = -a to a:
> restart:
r:=sqrt(x^2+y^2);
Jsquare:=int(int(r^2,x=-a..a),y=-a..a); evalf(%);
which is larger than the actual torsion constant of a square. This means that square cross sections don’t twist in a way that keeps straight lines straight before and after twisting, but instead find some easier way to twist by warping a bit.
Circular cross sections in FlexPDE
FlexPDE can create circles using the arc command rather than the line to commands. e.g.,
START(-diam/2,0)
arc(center=0,0) angle=360 !Note that the angle spec is in degrees by default, and note the American spelling of ‘center’
LINE TO CLOSE !Note that ‘line’ is unnecessary; “TO CLOSE” does the same thing
Polar coordinate calculation: relating twist angle to coordinate displacement and strain
Consider a beam extending in z and twisting about the z-axis. A cross section of the beam in the xy-plane has points with original Cartesian coordinates (x,y) related to original polar coordinates (r,) via
If the beam at this particular cross section has a twist angle of , then a point originally at in polar coordinates moves to ; that is, it stays at the same distance from the rotation axis (which we’re taking to be the origin) but rotates to a new angle :
This set of new coordinates corresponds to Cartesian displacements of:
In pure twisting from a uniform z moment along the beam extending in z with the z = 0 surface fixed, there’s no z displacement (if we ignore warping), so the shear strains are expressed via only x & y displacements:
Now we can sub in the u & v we found earlier to get how the stress is related to twist angle:
The torque about the centre produced by this shear is
Looking at this the other way around, it’s possible to determine the angle that occurred once you know the displacement:
So in FlexPDE we can find the twist angle at any cross section by:
One final note on small displacements – returning to , if is very small, these simplify as follows:
So, for small twist angles, the y displacement is proportional to the x position while the x displacement is proportional to the negative y-position:
(After rotation by a small angle, the point shifts up proportional to its original x value, and shifts left proportional to its original y value)
Example: Beam torsion
A beam centred on the z-axis with E = 5 GPa, = 0.3, and Lz = 4 m is fixed on its z = 0 surface and subjected to a torsional shear stress on its z = Lz surface given by . Use Maple and FlexPDE to determine the twist angle that results if the beam cross section is:
a) a circle with a 1 m diameter.
b) a square with a 1 m side length.
Objects before rotation:
Because the beam is centred on the z-axis (i.e., extends the same distance in x&-x and in y&-y) this shear stress distribution creates no net shear force and the torsional displacement should be all we need to account for. At any point (x,y), this shear force distribution creates a torque about(0,0) on that face given by:
> restart:
E:=5e9: nu:=.3: G:=E/(2*(1+nu)): Lz:=4:
tauzx:=-y*1e6: tauzy:=x*1e6:
R:=1/2: #radius
Mz:=int(int(y*(-tauzx)+x*tauzy, y=-sqrt(R^2-x^2)..sqrt(R^2-x^2)), x=-R..R);
JT:=3.14159/2*R^4; #Torsion constant
thetaTip:=Mz*Lz/(G*JT);
Mz:=int(int(y*(-tauzx)+x*tauzy, y=-R..R), x=-R..R);
JT:=2.25*R^4; #Torsion constant
thetaTip:=Mz*Lz/(G*JT);
In FlexPDE we can calculate the angle using
(magnification of 1000x)
(magnification of 1000x)
Interestingly, Flex predicts an increase in area of the top face because of this twisting but nevertheless correctly predicts the angle that the top face would rotate to.
Flex Code:
‘Eg9.1 circular torsion’
errlim=1e-8
spectral_colors
COORDINATES
cartesian3
u !Displacement in x
v !Displacement in y
w !Displacement in z
DEFINITIONS
diam=1 !diameter
G=E/(2*(1+nu))
r = sqrt(x^2+y^2)
phi=atan2(y,x)
!Calculate the twist angle in flexPDE directly:
thetatest = atan2(y+v,x+u)-phi !twist angle, but possibly outside the range -Pi to Pi
theta = if(thetatest<-pi) then thetatest+2*pi else if (thetatest>pi) then thetatest-2*pi else thetatest !twist angle
mag = 1000 !magnification
C11 =E*(1-nu)/((1+nu)*(1-2*nu))
C12 = E*nu/((1+nu)*(1-2*nu))
!Axial Strain
!Engineering Shear Strain
gxy=(dx(v)+dy(u))
gyz=(dy(w)+dz(v))
gxz=(dz(u)+dx(w))
!!Stress via Hooke’s law
!Axial Stress
sx = C11*ex+C12*ey+C13*ez
sy = C21*ex+C22*ey+C23*ez
sz = C31*ex+C32*ey+C33*ez
!Shear stress
u: dx(sx)+dy(sxy)+dz(sxz)=0
v: dx(sxy)+dy(sy)+dz(syz)=0
w: dx(sxz)+dy(syz)+dz(sz)=0
surface ‘bottom’ z=0
surface ‘top’ z=Lz
BOUNDARIES
surface ‘bottom’
value(u)=0
value(v)=0
value(w)=0
surface ‘top’
load(u)=-(y)*1e6
load(v)=(x)*1e6
START(-diam/2,0)
arc(center=0,0) angle=360
!Square cross section
{ start(-diam/2,-diam/2)
line to(diam/2,-diam/2) to (diam/2,diam/2) to (-diam/2,diam/2) to close}
grid(x+u, y+v, z+w)
grid(x+u*mag, y+v*mag, z+w*mag)
contour(sxz) painted on surface z=Lz
contour(syz) painted on surface z=Lz
contour(sz) painted on surface z=Lz
contour(u) painted on surface z=Lz
contour(theta) painted on surface z=Lz
elevation(u,v,w) from (0,0,0) to (0,0,Lz)
elevation(theta) from (diam/4,0,0) to (diam/4,0,Lz)
report val(u,diam/2,0,Lz)
report val(v,diam/2,0,Lz)
report val(w,diam/2,0,Lz)
report val(theta,diam/2,0,Lz)
Repeat the previous example with E = 5 MPa rather than 5 GPa
Now Maple gives 1000x the angles, but FlexPDE does not; here’s the non-magnified plots:
(no magnification!)
(no magnification!)
-→ rad rather than the 2.4 rad predicted by Maple.
This seems like it has to do with the coordinate system and the degree of rotation; the stresses applied at the top don’t rotate as the beam rotates:
surface ‘top’
load(u)=-(y)*1e6
load(v)=(x)*1e6
value(w)=0
These are always in the same directions even after the beam rotates meaning if the rotation is significant it’s not always in the direction around the beam. This means in this chapter FlexPDE will only be correct for small twist angles.
That is, when we take some material and apply forces to create torsion:
When the object twists, we probably meant for the stresses to twist with it if we want to keep making just torsional stresses:
But FlexPDE isn’t a mind reader; it keeps the stresses in the direction we pointed them:
This means if the rotation is significant the stresses will start just pulling the beam apart.
This also means if you reduce stiffness even more (say, to 5e3 Pa), the beam will reach a limit of rotation at an angle where they’d no longer produce a torque; 45°, i.e., 1.57 rad:
Confirmed!
Moment from Distributions and Vice- a +x-surface, a shear stress distribution of creates a torque about (y,z)=(0,0) of:
Note: To see why the directions and signs are this way consider the following image:
A positive x-torque about the origin would be created by + if it’s at -z (i.e., the purple ) or by if it’s at +y. The red and orange shears are at locations that create torques in the -x direction about the origin.
If each shear stress distribution is linear and centred on the centroid : , then we obtain:
This torque is independent of position and therefore creates a pure moment in the x direction.
/*******Proof:
Centroids are at:
Now consider the torque about an arbitrary position:
The second integral vanishes because
and therefore for any z0.
*********/
Therefore,
and therefore a shear distribution that produces the moment Mx is , where is the distance from the centroid on the x surface.
For a z-surface, we could then say that the stress distribution that produces a given is , where .
e.g.: Finding a shear stress distribution for a moment
Use the above formula to find the shear stress produced on a rectangular cross section Lx by Ly centred at x=y=0 by moment M, and check that it does produce the stated moment.
> restart:
#Lx:=.3: Ly:=.1:
r:=sqrt(x^2+y^2):
tauzx:=-y*M/int(int(r^2,x=-Lx/2..Lx/2),y=-Ly/2..Ly/2);
tauzy:=x*M/int(int(r^2,x=-Lx/2..Lx/2),y=-Ly/2..Ly/2);
Mz:=simplify(int(int(-y*tauzx+x*tauzy,y=-Ly/2..Ly/2),x=-Lx/2..Lx/2));
Continuity
Complicated pure torsion problems often involve asserting piecewise continuity in displacement, force, and torque and working out how they can change when point-moments are applied. e.g., consider the following problem:
(Take G for this steel to be 80 GPa).
In this problem, twist angle and internal moment (same as internal torque about the axis) will depend on position. Define the positive direction for each at any point as down and to the right; here’s the original isometric drawing, a view of the gears from the side, and a top view:
Internal moment Mz is defined at any point as:
Assertions:
-Torsion follows Hooke’s law for torsion:
-Twist angle at any radius within a gear is the same, and displacement of the gear teeth is the same for each gear, meaning
-Torque produced by the gear about the shaft is , and force applied at the gear teeth by C on D is equal and opposite to that applied by D on C. Assigning the upward direction of force (out of the page for the top view) at that intersection as positive, we have:
(Note that while the forces between the teeth, FCD & FDC, are Newton’s third law action-reaction pairs, the torques of each gear on the other TCD& TDC aren’t; they end up pointing in the same direction and having different magnitude from the gear ratio.)
-The angular displacement at point F is 0 (but it can apply any torque), while the torque applied by point A is zero (but it can have any rotational angle).
Internal moment is 0 from A to A’, then between C and A’ it’s 12 N-m, which means that , and therefore . After the point moment applied at G, this means the moment between F and G must be .
Or, directly,
We can work our way to the same answer using Maple in perhaps a more robust approach. Assigning z = 0 at point C for ABC, we have:
> restart:
rD:=.12: rC:=.04: rshaft:=.03/2:
JT:=3.14159/2*rshaft^4: G:=80e9:
MzDF:=TF-piecewise(z>1.5, 4)-piecewise(z>2, TCD):
MzAC:=-TDC-piecewise(z>1.25, 12):
evalf(subs(z=1.5, MzAC)=0); #No torque at point A;
TDC:=solve(%);
TCD:=rD/rC*TDC; #steps up the torque by the radius
evalf(subs(z=2.1, MzDF));
TF:=solve(%);
plot([MzDF, MzAC], z=0..2);
Confirms our earlier calculation of internal moments; it shows the moment for DF is much larger, and largest near D, while that for AC is largest between C and A’.
The peak shear stress is at the surface: , so
> tauDFMax:=evalf(subs(z=1.6, MzDF)*rshaft/JT);
tauACMax:=evalf(subs(z=0, MzAC)*rshaft/JT);
To work out the angles, use , and start with for DF:
> thetaDF:=int(MzDF/(JT*G),z=0..z):
thetaD:=evalf(subs(z=2, thetaDF));
Then use the gear ratio to determine the angle of twist this causes at C: . This is the starting point for the twist angle from C to A:
> thetaC:=-rD/rC*thetaD;
thetaAC:=thetaC+int(MzAC/(JT*G),z=0..z):
thetaA:=evalf(subs(z=1.5, thetaAC));
plot([thetaDF, thetaAC], z=0..2);
Note that the maximum shear stress that occurs here may not be in the rods at all, but in the teeth of the gear. Suppose that for this problem the entire load is carried by only one gear tooth, and the gear tooth base area is 1 cm thick by 2 cm width. Then the gear tooth base shear is
> taugt:=TCD/rD/(.01*.02);
which is lower than the peak shear on the shaft in this case.
Combination Elasticity
For arbitrary loading a beam will experience axial, shear, flexure, and torsion strain simultaneously. If strains stay small enough for Hooke’s law to be valid then the linearity of Hooke’s law means superposition is valid: To find the total deformation a beam experiences from arbitrary loading, add the displacements calculated from axial, shear, flexure, and torsion strains together (“superimpose” them).
Example of Combination Elasticity
e.g., a cube with E = 70 GPa and = 0.3 of side-length L = 1 m extending from (0,0,0) to (L,L,L) is fixed on its x = 0 surface and subjected to a force of and a moment of , both evenly distributed over its surface.
a) Find the reaction force and moment produced at the support
b) Find the displacement of the (L,L,L) corner of the cube.
Solution Theory
(directions shown in this image are arbitrary). To find the reaction force and moment at the x = 0 surface, we can use Newton’s second law in statics, which now includes 3 force and 3 torque equations.
solving these for RF and RM gives:
> restart:
with(LinearAlgebra):
Fend:=<200,300,400>:
Mend:=<300,-200,500>:
RF:=-Fend;
RM:=-CrossProduct(<1,0,0>,Fend)-Mend;
At any point x from the fixed support to the end we can split the cube and determine the internal force and moment:
> Fin:=-RF;
Min:=-RM-CrossProduct(
This internal force and moment then let us determine the displacements using our elasticity information:
Forces directly relate to average axial and shear stresses ,
which then relate to strains through Hooke’s law: (if we ignore the axial coupling of other axis-stress into x-stress via Poisson’s ratio, otherwise we’d get additional strains through ; notably if , we can directly calculate . Note that these don’t increase along the length of the beam if is constant.)
Moments cause either bending or twisting depending on the direction. For an x-surface,
How much displacement the twist angle corresponds to in each direction depends on the specific point on the surface Ax we’re talking about:
, where here and are the coordinates of the centre of rotation
(which in this problem is by symmetry), and we’re using a small angle approximation.
Maple Implementation
For this case of the cube then,
> restart: with(LinearAlgebra):
Fend:=<200,300,400>:
Mend:=<300,-200,500>:
RF:=-Fend;
RM:=-CrossProduct(<1,0,0>,Fend)-Mend;
Min:=-RM-CrossProduct(
L:=1: E:=70e9: nu:=.3: G:=E/(2*(1+nu)):
sx:=Fin[1]/Ax; tauxy:=Fin[2]/Ax; tauxz:=Fin[3]/Ax;
u:=int(sx/E,x=0..x); utip:=subs(x=L,u);
vaxialtip:=L/2*(-nu)*sx/E;
waxialtip:=L/2*(-nu)*sx/E;
vshear:=int(tauxy/G, x=0..x);
vsheartip:=subs(x=L, vshear);
wshear:=int(tauxz/G, x=0..x);
wsheartip:=subs(x=L, wshear);
JTx:=2.25*(L/2)^4:
theta:=int(Min[1]/(JTx*G),x=0..x);
thetatip:=subs(x=L, theta);
vtorsiontip:=-L/2*thetatip;
wtorsiontip:=L/2*thetatip;
EIz:=E*L^4/12: EIy:=EIz:
wflexure:=int(int(-Min[2]/EIz,x=0..x),x=0..x);
wflexuretip:=subs(x=L, wflexure);
vflexure:=int(int(Min[3]/EIy,x=0..x),x=0..x);
vflexuretip:=subs(x=L, vflexure);
vtip:=vaxialtip+vsheartip+vflexuretip+vtorsiontip;
wtip:=waxialtip+wsheartip+wflexuretip+wtorsiontip;
FlexPDE Implementation
To check this in FlexPDE, we need to determine how to apply the moments via stress distributions on the tip. From the work in beam bending, a z-moment is produced by an additional x-stress of:
Similarly, a y-moment is produced by an extra x-stress of
(the sign is reversed because -y is to z as +z is to y).
And from our work above, an x-moment is caused by shear stress distributions of
These give:
> sxp:=-Mend[3]*(y-L/2)/int((y-L/2)^2,y=0..L)+Mend[2]*(z-L/2)/int((z-L/2)^2,z=0..L);
Mzcheck=int(int(-sxp*y,y=0..L),z=0..L);
Mycheck=int(int(sxp*z,y=0..L),z=0..L);
zbar:=L/2: ybar:=L/2:
rc:=sqrt((z-zbar)^2+(y-ybar)^2):
tauxy:=-Mend[1]/int(int(rc^2,y=0..L),z=0..L)*(z-zbar);
tauxz:=Mend[1]/int(int(rc^2,y=0..L),z=0..L)*(y-ybar);
Mxcheck=int(int(tauxy*(-z)+tauxz*y,y=0..L),z=0..L);
Coding in FlexPDE gives good agreement with Maple for the v and w at (L,L,L) but not so good for u:
Epilogue: Exploring why Flex and Maple don’t agree
Some reasons u isn’t in as much agreement in our Maple calculation are:
1. The warping that occurs for twisting non-circular cross-sections can cause local increases or decreases in longitudinal displacement,
2. The boundary conditions fixing v&w to 0 at the wall translate to extra u displacements showing up, and
Maple ignored local variation in in calculating displacement
-> In fact, the u we calculated is very close to the u in the centre of the x = Lx surface:
this means that reason 2 is not a large effect in this problem. Reason 3 is expected to be most of the issue because there are local variations that are quite large due to the extra stress we needed to apply to create the bending moments:
To remove the effect of bending moments and compare, repeat the problem with Mend[2]=Mend[3]=0:
This has removed about half of the difference between the u calculated using flex and maple:
Note that the u-displacement ultimately results from local stress variation within the beam along the length. While the x-stress at x = Lx is uniform, the x-stress at an intermediate position is not:
This is becaus
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com