0
$\begingroup$

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[{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

Improve this question
$\endgroup$
  • $\begingroup$ Can't copy this [{\!\( \*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:35
  • $\begingroup$ Thanks for the comment. I don't know what's going on. I copied it into my Mathematica and it worked fine. Let me test and I will edit. $\endgroup$ – wpilgri May 21 '20 at 20:54
  • $\begingroup$ Before you copy it, make a copy of the cell to another empty notebook (so you do not lose your original formating) Then convert the new cell to INPUT FORM first. Then copy the result as plain text to here. $\endgroup$ – Nasser May 21 '20 at 20:55
  • $\begingroup$ I edited so I hope it works better now. Thanks for pointing out my error! $\endgroup$ – wpilgri May 21 '20 at 21:00
  • $\begingroup$ It is still not working. Reverse[a] there is no a 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
2
$\begingroup$

After small modification we have coincidence for calcdeltaL and deltaL with all digits:

Y = 199*^9;(*Young's modulus in Pa*)\[Rho] = 7860;(*Steel density in \
kg/m^3*)dia = 1/39.37; nu = .0;(*1" dia converted to meters*)c = 
 Sqrt[Y/\[Rho]]; g = 9.81;
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, t] - c^2 D[z[x, t], {x, 2}] - g == 
    NeumannValue[0, x == len], 
   DirichletCondition[z[x, t] == 0, x == 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}, 
 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}, 
  PlotLabel -> "PDE Val and Correct (dashed line)", PlotRange -> All],
  ListLinePlot[xydata, PlotStyle -> {Dashed, Green}]]

Length 0.193735 and 0.193735. Visualization Figure 1

Improve this answer
$\endgroup$
  • $\begingroup$ Absolute genius @Alex-Trounev Thank you very much. I need to go read why g replaces my distributed load but clearly you solved in minutes what I spent countless hours on. Greatly appreciated....Thank you, thank you!!! $\endgroup$ – wpilgri May 22 '20 at 1:43
  • $\begingroup$ Why do you solve this as a stationary PDE and not a time dependant one? $\endgroup$ – user21 May 22 '20 at 3:46
  • $\begingroup$ @wpilgri You are welcome! $\endgroup$ – Alex Trounev May 22 '20 at 11:19
  • $\begingroup$ @user21 For time dependence there is no reason, since we are looking stationary solution to compare with textbook. $\endgroup$ – Alex Trounev May 22 '20 at 11:31
  • $\begingroup$ @user21 Baby steps. I tried to solve it and look at the stationary solution before getting more complicated. If I couldn't get stationary right, why proceed is my thinking. $\endgroup$ – wpilgri May 22 '20 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.