I am trying to setup Mathematica to analyze a vertical round rod under its own weight, fixed on one end free on the other. I have the 1D wave equation and a distributed load to represent the self weight of the round rod.
The problem is when I compare the Mathematica solution to the textbook solution the two do not agree.
Sample problem is given below.
Y = 199*^9; (*Young's modulus in Pa *)
\[Rho] = 7860; (* Steel density in kg/m^3*)
dia = 1/39.37; (* 1" dia converted to meters*)
c = Sqrt[Y/\[Rho]];
len = 1000; (*length in meters*)
tmax = 5; (* Max time for analysis*)
area = \[Pi]*dia^2/4; (*Round rod cross sectional area*)
wtfactor = \[Rho]*9.81*area/len;
frwt[x_] := \[Rho]*
area*9.81*(1 -
x/len); (*Rod Self weight imposed as a distributed load*)
nsol6 = NDSolve[{D[z[x, t], {t, 2}] ==
c^2*D[z[x, t], {x, 2}] + frwt[x] + NeumannValue[0, x == len],
z[0, t] == 0},
z[x, t], {x, 0, len}, {t, 0, tmax},
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 10}}]
fnnsol6[x_, t_] = nsol6[[1, 1, 2]]
Plot3D[fnnsol6[x, t], {x, 0, len}, {t, 0, tmax},
PlotLabels -> Automatic, AxesLabel -> Automatic]
deltaL = ((\[Rho]*9.81*len^2)/(
Y*2)) (*Textbook elongation for a vertical rod under self weight*)
calcdeltaL =
fnnsol6[len,
5] (*Calculated delta Length from PDE solution. Should match \
textbook*)
deltaLfunc[x_, l_] := \[Rho]*9.81*
x*(2*len - x)/(2*Y) (*Verified Correct*)
xydata = Thread[{Range[0, 1000, 100],
deltaLfunc[x, 1] /. {x -> Range[0, 1000, 100]}}];
Show[Plot[fnnsol6[x, 0], {x, 0, len}, PlotLabels -> {"PDE Val"},
PlotRange -> All
],
ListLinePlot[xydata, PlotStyle -> Green, PlotLabels -> {"Correct"}]]
If you've read this far, thank you.
In summary my question is: Is this a Mathematica issue or a PDE setup problem? The PDE is right out of a textbook so I don't think that's the problem but Mathematica gives no errors and I am out of troubleshooting ideas so looking for some help.
Thank You
[{\!\( \*SubscriptBox[\(\[PartialD]\), \({t, 2}\)]\(z[x, t]\)\) == c^2*\!\( \*SubscriptBox[\(\[PartialD]\), \({x, 2}\)]\(z[x, t]\)\) + frwt[x] +
it shows up all mangled on my notebook. try to copy/paste using plain text from the InputForm version of the cell. Make sure to make copy of it first. $\endgroup$ – Nasser May 21 '20 at 20:35Reverse[a]
there is noa
anywhere to reverse. It is always best to try your code before posting it from clean kernel. $\endgroup$ – Nasser May 21 '20 at 21:02