1D Wave Equation: Vertical Rod and Displacement vs. Textbook Solution


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.

Vertical Rod Layout

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[{\!\( \*SubscriptBox[\(\[PartialD]\), \({t, 2}\)]\(z[x, t]\)\) == c^2*\!\( \*SubscriptBox[\(\[PartialD]\), \({x, 2}\)]\(z[x, t]\)\) + 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]}}]; xydata2 =   Thread[{Range[0, 1000, 100],     Reverse[a]}]; (*Same answer different calc format for debugging*) Show[Plot[fnnsol6[x, 0], {x, 0, len}, PlotLabels -> {"PDE Val"},    PlotRange -> All   ],  ListLinePlot[xydata2, 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