I am solving a simple structural beam problem using slope-deflection approach - 3 support A(left- fixed), B(middle- pin), C(right - pin). Beam AB - uniform load, Beam BC - concentrated load in middle.
I was successful in getting the complete full length beam moment equations using a piecewise and unitstep function. I was able to get the moment at specific points along the beam length and get plots comparing the 2 approaches.
Now when I need to take the 1 st derive to get a shear plot or evaluate the shear at a specific point, the derived function does not work - no matter if I take the derivative of either moment function using the Piecewise or UnitStep functions. ?? For some reason, the derived function is not considering the x distance point as a variable ? My ultimate goal is to get shear, deflection, slope plots and evaluate the values as specific points.
(*3 supports A(left end- fixed), B(middle-pinned), C(right end-pinned), uniform and concentrated load *)
qc = 1;
qu = 1;
lab = 10;
lbc = 10;
deltaab = 0;
deltabc = 0;
anglea = 0;
ei = 1;
momableft = (2.*ei)/lab*(2*anglea + angleb - (3*deltaab)/lab) + (-qu*
lab^2)/12;
mombaright = (2.*ei)/lab*(anglea + 2*angleb - (3*deltaab)/lab) + (+qu*
lab^2)/12;
mombcleft = (2.*ei)/lbc*(2*angleb + anglec - (3*deltabc)/lbc) + (-qc*
lbc)/8;
momcbright = (2.*ei)/lbc*(angleb + 2*anglec - (3*deltabc)/lbc) + (+qc*
lbc)/8;
ssol1 = Solve[
(mombaright + mombcleft) == 0 &&
(momcbright) == 0 ,
{angleb, angleb, anglec} ];
anglebValue = angleb /. ssol1[[1]];
anglecValue = anglec /. ssol1[[1]];
angleaValue = anglea;
Print[" The angle for support A: (ei=1) ", angleaValue]
Print[" The angle for support B: (ei=1) ", anglebValue]
Print[" The angle for support C: (ei=1) ", anglecValue]
momableft = (2.*ei)/
lab*(2*angleaValue + anglebValue - (3*deltaab)/lab) + (-qu*lab^2)/
12;
mombaright = (2.*ei)/
lab*(angleaValue + 2*anglebValue - (3*deltaab)/lab) + (+qu*lab^2)/
12;
mombcleft = (2.*ei)/
lab*(2*anglebValue + anglecValue - (3*deltabc)/lbc) + (-qc*lab)/
8;
momcbright = (2.*ei)/
lab*(anglebValue + 2*anglecValue - (3*deltabc)/lbc) + (+qc*lab)/8;
Print[" ********* Moments beam A-B & B-C ********* "];
Print[" Moment AB: ", momableft, " ft-kips"];
Print[" Moment BA: ", mombaright, " ft-kips"];
Print[" Moment BC: ", mombcleft, " ft-kips"];
Print[" Moment CB: ", momcbright, " ft-kips"];
********* Moments beam A-B & B-C *********
Moment AB: -10.1786 ft-kips
Moment BA: 4.64286 ft-kips
Moment BC: -4.64286 ft-kips
Moment CB: 0. ft-kips
(* now for each beam in the model and using the end moments just \
calculated and loading on the beam , determine the reaction loads *)
solutionA =
Solve[(momableft + mombaright - (qu*lab^2)/2 +
rAbeamabqusolve*lab) == 0,
{rAbeamabqusolve}]; (*moment about point B *)
rAbeamabqusolveValue = rAbeamabqusolve /. solutionA[[1]];
solutionB =
Solve[(momableft + mombaright + (qu*lab^2)/2 -
rBbeamabqusolve*lab) == 0,
{rBbeamabqusolve}]; (*moment about point A *)
rBbeamabqusolveValue = rBbeamabqusolve /. solutionB[[1]];
Print[" "];
Print[" ********* Reaction beam A-B ********* "];
Print[" Beam A-B, REACTION A: ", rAbeamabqusolveValue, " kips"];
Print[" Beam A-B, REACTION B: ", rBbeamabqusolveValue, " kips"];
Print[" ZERO FORCE BALANCE CHECK ",
rAbeamabqusolveValue + rBbeamabqusolveValue - qu*lab, " kips"];
solutionBbc =
Solve[(mombcleft + momcbright - (qc*lbc)/2 + rBbeambcqcsolve*lbc) ==
0,
{rBbeambcqcsolve}]; (*moment about point C *)
rBbeambcqcsolveValue = rBbeambcqcsolve /. solutionBbc[[1]];
solutionCbc =
Solve[(mombcleft + momcbright + (qc*lbc)/2 - rCbeambcqcsolve*lbc) ==
0,
{rCbeambcqcsolve}]; (*moment about point B *)
rCbeambcqcsolveValue = rCbeambcqcsolve /. solutionCbc[[1]];
Print[" "];
Print[" ********* Reaction beam B-C ********* "];
Print[" Beam B-C, REACTION B: ", rBbeambcqcsolveValue, " kips"];
Print[" Beam B-C, REACTION C: ", rCbeambcqcsolveValue, " kips"];
Print[" ZERO FORCE BALANCE CHECK ",
rBbeambcqcsolveValue + rCbeambcqcsolveValue - qc , " kips"];
********* Reaction beam A-B *********
Beam A-B, REACTION A: 5.55357 kips
Beam A-B, REACTION B: 4.44643 kips
ZERO FORCE BALANCE CHECK 0. kips
********* Reaction beam B-C *********
Beam B-C, REACTION B: 0.964286 kips
Beam B-C, REACTION C: 0.0357143 kips
ZERO FORCE BALANCE CHECK 0. kips
(* TOTAL MODEL REACTIONS *)
Ra = rAbeamabqusolveValue;
Rb = rBbeamabqusolveValue + rBbeambcqcsolveValue;
Rc = rCbeambcqcsolveValue;
Print[" TOTAL MODEL REACTIONS A: ", Ra, " kips"]
Print[" TOTAL MODEL REACTIONS B: ", Rb, " kips"]
Print[" TOTAL MODEL REACTIONS C: ", Rc, " kips"]
Print[" TOTAL MODEL REACTION check: ",
Ra + Rb + Rc - qc - qu*lab, " kips"]
TOTAL MODEL REACTIONS A: 5.55357 kips
TOTAL MODEL REACTIONS B: 5.41071 kips
TOTAL MODEL REACTIONS C: 0.0357143 kips
TOTAL MODEL REACTION check: 1.77636*10^-15 kips
(*Define the total beam model functions *)
momABCtotal[x_] := Piecewise[{
{momableft - (qu*x^2)/2 + (Ra*x), x < lab},
{momableft + mombaright +
mombcleft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab),
x == lab },
{momableft + mombaright +
mombcleft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab),
x > lab && x < lab + (lbc/2)},
{momableft + mombaright + mombcleft +
momcbright - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab) +
Rc*(x), x == lbc + lab}
}]
momABCtotalUSTEP3[x_] :=
(momableft - (qu*x^2)/2 + (Ra*x))*UnitStep[lab - x] +
(momableft - (qu*lab)*(x - lab/2) + (Ra*x) + Rb*(x - lab))*
UnitStep[x - lab]*UnitStep[(lab + lbc/2) - x] +
(momableft - (qu*lab)*(x - lab/2) - qc*(x - (lab + lbc/2)) + Ra*x +
Rb*(x - lab))*UnitStep[x - (lab + (lbc/2))] +
UnitStep[(lab + lbc) - Infinity]
momABCtotal[3] (* see if function can evaluate value at x = 3 \
point *)
momABCtotalUSTEP3[3]
Plot[momABCtotal[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Bending_Moment (k - ft)[x]}]
Plot[momABCtotalUSTEP3[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, USTEP3_Bending _Moment (k - ft)[x]}]
Plot[{momABCtotal[x], momABCtotalUSTEP3[x]}, {x, 0, lab + lbc},
PlotLegends -> {"Piecewise", "UnitStep"}, PlotRange -> All]
shearmomABCtotal[x_][x_] := D[momABCtotal[x], x]
shearABCtotalUSTEP3[x_] := D[momABCtotalUSTEP3[x], x]
shearmomABCtotal[3] (* see if function can evaluate value at x = \
3 point *)
shearABCtotalUSTEP3[3]
simplifiedshearmomABCtotal[x_] := Simplify[shearmomABCtotal[x]]
simplifiedShearABCtotalUSTEP3[x_] := Simplify[shearABCtotalUSTEP3[x]]
simplifiedshearmomABCtotal[3]
simplifiedShearABCtotalUSTEP3[3]
(*Evaluate the shear force function at a specific point*)
shearValueAt1 = simplifiedShearABCtotalUSTEP3[1]
shearValueAt9 = simplifiedShearABCtotalUSTEP3[9]
(*Output the value*)
Print["shearValueAt1 = ", shearValueAt1]
Print["shearValueAt9 = ", shearValueAt9]
(*Display the shear force function*)
simplifiedShearABCtotalUSTEP3[x]
(*Evaluate the shear force function at specific points*)
shearValueAt9 = simplifiedShearABCtotalUSTEP3[9]
shearValueAt4 = simplifiedShearABCtotalUSTEP3[4]
(*Plot the shear force function*)
Plot[simplifiedShearABCtotalUSTEP3[x], {x, 0, 20},
PlotLabel -> "Shear Force Function", AxesLabel -> {"x", "V(x)"}]
Plot[simplifiedShearABCtotalUSTEP3[x_], {x, 0, 20},
PlotLabel -> "Shear Force Function", AxesLabel -> {"x", "V(x)"}]
Plot[shearForce[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Shear (k)[x]}]
Plot[shearABCtotalUSTEP3[x], {x, 0, lab + lbc}, PlotRange -> All,
AxesLabel -> {x, Shear (k)[x]}]
Plot[simplifiedShearABCtotalUSTEP3[x], {x, 0, lab + lbc},
PlotRange -> All, AxesLabel -> {x, Shear (k)[x]}]
moment plots are correct and match hand calcs.
g[x_] := D[f[x], x]; however, when calledxhas a value and things likeD[f[3], 3]make no sense. Instead use the formg[x_] := Derivative[1][f][x]so that the derivative is taken before the value is assigned to the variable. Alternatively, useSetinstead ofSetDelayed, e.g.,g[x_] = D[f[x], x]$\endgroup$