Longitudinal Bar Vibration: Prescribed BC's at end not working

I am working on a vertical 1D bar vibration problem using the wave equation (ie: oil industry sucker rod application). On one end of the bar I have prescribed displacement and force boundary conditions. The plan is to use this model to calculate what the load and motion would be at the other end of the rod, which in the real world in many thousands of feet underground and you can't measure directly.

For debugging and "let's start simple" purposes I applied the prescribed force and displacement boundary conditions on the x=0 end. Displacement as Dirichlet and Force and NeumannValue. NDSolve solves without error, but when I take the derivative of the solution z[x,t] with respect to x to get strain, and calculate the force (Force=Strain x Modulus x Area) and compare it to the input boundary condition the two do not match so I am clearly doing something wrong. Either my BC is setup incorrectly or I am calculating the strain improperly. I am hoping someone can help point out the error I have.

Y = 199*^9; (*Pa*)
\[Rho] = 7860; (*kg/m^3*)
dia = 1/39.37; (* 1" dia converted to meters*)
c = Sqrt[Y/\[Rho]];
g = 9.81;
area = \[Pi]*dia^2/4;

endPrescribedDisp[t_] := Sin[t]; (*Prescribed Displacement BC*)
endPrescribedForce[t_] := 100*Cos[t]; (*Prescribed Force BC*)

(*prescribed Force BC.  Translated to strain for NeumannValue BC*)
eq1 = D[z[x, t], {t, 2}] - c^2*D[z[x, t], {x, 2}] - g ==
NeumannValue[endPrescribedForce[t]/(Y*area), x == 0]
bc = DirichletCondition[z[x, t] == endPrescribedDisp[t], x == 0]

sol = NDSolve[{eq1, bc, z[x, 0] == 0},
z[x, t], {x, 0, 1}, {t, 0, 2*\[Pi]}]
solf[x_, t_] := sol[[1, 1, 2]]

Plot3D[solf[x, t], {x, 0, 1}, {t, 0, 2*\[Pi]},
PlotLabel -> "Solution Displacement"]

strain[x_, t_] :=
Evaluate[D[solf[x, t],
x]] (*Take derivative of solution to get strain*)

Plot[strain[0, t], {t, 0, 2*\[Pi]},
PlotLabel -> "Strain at Prescribed End"]
Plot[strain[0, t]*Y*area, {t, 0, 2*\[Pi]},
PlotLabel -> "Calculated Force at Prescribed End"]
Print["BC force at sample time=", endPrescribedForce[\[Pi]]]
Print["Calculated BC force at sample time=",
strain[0, \[Pi]]*Y*area, ".  Doesn't match above!"]



Greatly Appreciated!

• This is a wave equation, you need 2 initial conditions. Also, you've imposed 2 b.c. at $x=0$, which is a bit unusual. Currently what you're actually doing is something similar to this: mathematica.stackexchange.com/q/222371/1871 Jun 29, 2020 at 1:31
• Don't believe you can prescribe both force and displacement at same point. Jun 29, 2020 at 2:52
• In an actual pumping oil well, you measure both position and load vs time at the top of the rod string, so yes, you start with force and displacement at the same point, but they are not strictly independent. The load depends on displacement amplitude, speed, damping, the weight of rod string and the amount of fluid being lifted. You use the wave equation to compute the position and load at the bottom of the rod string where the pump resides, starting with the given conditions at the top. If that is what the OP is really after, a free end bc at the bottom will not give the desired result. Jun 30, 2020 at 0:39
• Another related post: mathematica.stackexchange.com/q/156651/1871 Jun 30, 2020 at 3:07

It has been a long time since I've thought about this problem. It is probably hard to solve with NDSolve, but it works with finite differences. Start with conditions on the surface and work your way down the rod string to the bottom.

General wave equation with viscous damping and gravity.

pde = D[u[x, t], t, t] == a^2*D[u[x, t], x, x] - c*D[u[x, t], t] - g


u = displacement of rod x = position on the rod string a = speed of sound in the rod string (wave velocity) c = viscous damping coefficient g = acceleration due to gravity F = Load (Tension at a point on the rod)

Finite differences as a rule.

fd = {dttu[i, n] -> ( u[i, n + 1] - 2 u[i, n] + u[i, n - 1])/Δt^2,
dxxu[i, n] -> ( u[i + 1, n] - 2 u[i, n] + u[i - 1, n])/Δx^2,
dtu[i, n] -> (u[i, n + 1] - u[i, n - 1])/(2 Δt),
dxu[i, n] -> (u[i + 1, n] - u[i - 1, n])/(2 Δx),
F[i, n] -> -EA (u[i + 1, n] - u[i - 1, n])/(2 Δx), a -> Δx/Δt, c -> γ/Δt}


Get the pde in finite difference form.

pdefd = dttu[i, n] == a^2 dxxu[i, n] - c dtu[i, n] - g /. fd
(*(u[i, n - 1] - 2*u[i, n] + u[i, n + 1])/Δt^2 == -((γ*(u[i, n + 1] - u[i, n-1]))/(2*Δt^2)) +
(u[i - 1, n] - 2*u[i, n] + u[i + 1, n])/Δt^2 - g*)

Solve[pdefd, u[i + 1, n]] // Flatten // Simplify // Collect[#, {u[i, n + 1], u[i, n - 1]}] &
(*{u[i + 1, n] -> (1/2)*(2 - γ)*u[i, n - 1] + (γ/2 + 1)*u[i, n + 1] - u[i - 1, n] + Δt^2*g}*)

{u[i + 1, n] -> (1/2 (2 - γ) // Expand) u[i, n - 1] + (γ/2 + 1) u[i, n + 1] - u[i - 1, n] + Δt^2 g}
(*{u[i + 1, n] -> (1 - γ/2) u[i, n - 1] + (γ/2 + 1) u[i, n + 1] -
u[i - 1, n] + Δt^2 g}*)


Get u[2,n]

u2nRule = % /. i -> 1
(*{u[2, n] -> (1 - γ/2)*u[1, n - 1] + (γ/2 + 1)*u[1, n + 1] - u[0, n] + Δt^2*g}*)


From the F finite difference

(u2nRule /. u[0, n] -> u[2, n] - (2 Δx)/EA F[n])[[1]] /. Rule -> Equal
(*u[2, n] == (1 - γ/2) u[1, n - 1] + (γ/2 + 1) u[1, n + 1] - u[2, n] + (2 Δx F[n])/ EA + Δt^2 g*)

Solve[%, F[n]] // Flatten // Expand // Collect[#, {u[1, n + 1], u[1, n - 1]}] &
(*{F[n] -> u[1, n - 1]*((γ*EA)/(4*Δx) - EA/(2*Δx)) + u[1, n + 1]*(-((γ*EA)/(4*Δx)) - EA/(2*Δx)) + (EA*u[2, n])/Δx - (Δt^2*EA*g)/(2*Δx)}*)


Simplify a little

F[n] == EA/(2 Δx) ((F[n] /. % ) (2 Δx)/EA // Expand) //
Collect[#, {EA/(2 Δx), u[1, n + 1], u[1, n - 1]}] &
(*F[n] == (EA*((γ/2 - 1)*u[1, n - 1] + (-(γ/2) - 1)*u[1, n + 1] + 2*u[2, n] + Δt^2*(-g)))/(2*Δx))


Manually change back to i's and n's. 2->i+1, 1->i

F[n] == (EA (-((1 - γ/2) u[i, n - 1]) - (γ/2 + 1) u[i, n + 1] + 2 u[i + 1, n] + Δt^2 (-g)))/(2 Δx);


Now that we have established the finite difference equations we can implement the finite difference procedure. We will use actual well data I have saved from a well measured in the early 80's.

data = {{0., 0., 12.689}, {0.335, 0.161, 13.933}, {0.446, 0.31, 14.928}, {0.538, 0.459, 15.177}, {0.611, 0.621, 15.177},
{0.865, 1.242, 15.799}, {1.049, 1.863, 18.66}, {1.217, 2.483, 20.65}, {1.368, 3.104, 21.77}, {1.503, 3.725, 22.765},
{1.641, 4.346, 24.009}, {1.78, 4.967, 24.507}, {1.899, 5.588, 22.89}, {2.016, 6.208, 20.65}, {2.132, 6.829, 19.531},
{2.266, 7.45, 18.411}, {2.382, 8.071, 18.162}, {2.498, 8.692, 17.789}, {2.615, 9.313, 17.789}, {2.755, 9.933, 19.655},
{2.916, 10.554, 20.277}, {3.077, 11.175, 20.402}, {3.275, 11.796, 20.028}, {3.334, 11.957, 19.531}, {3.413, 12.106, 19.033},
{3.509, 12.268, 17.54}, {3.717, 12.417, 16.172}, {3.921, 12.268, 15.301}, {3.995, 12.106, 15.052}, {4.068, 11.957, 15.052},
{4.122, 11.796, 14.928}, {4.286, 11.175, 12.938}, {4.415, 10.554, 10.574}, {4.528, 9.933, 8.832}, {4.625, 9.313, 7.34},
{4.723, 8.692, 6.22}, {4.822, 8.071, 5.598}, {4.921, 7.45, 6.469}, {5.017, 6.829, 9.206}, {5.11, 6.208, 11.694},
{5.202, 5.588, 13.062}, {5.293, 4.967, 14.182}, {5.404, 4.346, 15.052}, {5.496, 3.725, 15.177}, {5.61, 3.104, 14.182},
{5.749, 2.483, 11.569}, {5.909, 1.863, 10.325}, {6.091, 1.242, 10.45}, {6.327, 0.621, 13.186}, {6.406, 0.459, 13.808},
{6.504, 0.31, 14.928}, {6.621, 0.161, 13.933}, {6.889, 0., 12.689}};


Actual well data of the polish rod position and load with time. The polish rod is the top rod of the rod string and is above the surface.

Column 1 is the time in seconds

Column 2 is the position in feet.

Column 3 is the load in units of 1000 lbs.

These are practical oilfield units. Not metric

Assign the values to variables

timez = Table[data[[n, 1]], {n, Length[data]}];
posz = Table[data[[n, 2]], {n, Length[data]}];
loadz = Table[data[[n, 3]], {n, Length[data]}];
tmax = timez[[Length[data]]];


Interpolation functions for position and load

posT = Interpolation[Table[{timez[[n]], posz[[n]]}, {n, Length[data]}]];


Plot polish rod position vs polish rod load

ParametricPlot[{posT[t], loadT[t]}, {t, 0, 6.889},
AspectRatio -> 1/GoldenRatio, PlotRange -> {{0, 13}, {0, 25}},


The polish rod ( the top rod that is above ground) is moving clockwise in the above plot. Higher loads occur when the rod string is rising. It is jerky because it is real data.

More well data.

sg = 0.993;(* produced fluid specific gravity*)
a = 1.95538  10^4;(*speed of sound in steel, ft/sec*)
Ey = 3 10^4;(*Young's modulus KSI*)
g = 32.2 ;(*gravitation constant, ft/sec^2*)
n1 = 5; (*initial rod string segments*)
n2 = 5;(*alternate rod string setments if necessary*)
roddia = {1.0, 0.875,
0.75};(*the rod string consists of 3 segements with these \
diameters, inches*)
rodlen = {1950, 2025, 1893};(*rod string segment lengths, feet*)
c = 0.2;(*damping coefficient, 1/second*)
area = π roddia^2/4;(*square inches)


Now move down the rod string.

ll = 0;
nrod = Length[rodlen];
timez = Table[data[[n, 1]], {n, Length[data]}];
posz = Table[data[[n, 2]], {n, Length[data]}];
loadz = Table[data[[n, 3]], {n, Length[data]}]; posT =
Interpolation[Table[{timez[[n]], posz[[n]]}, {n, Length[data]}]];
Do[
ll = ll + rodlen[[m]];
If[m < nrod, buoy = 0.433 sg ll (area[[m]] - area[[m + 1]])/1000];
EA = Ey area[[m]];
Δx = -rodlen[[m]]/n1;
n3 = n2;
Δt = -Δx/a;
γ = c Δt;
j = 40;
k = Floor[tmax/Δt + j];
u = Table[0, {ii, 7}, {jj, k}];
t = Table[0, {ii, k}];
F = Table[0, {ii, k}];
Do[
t[[n]] = (n - j/2 - 1) Δt;
time =
If[t[[n]] >= 0, If[t[[n]] <= tmax, t[[n]], t[[n]] - tmax],
t[[n]] + tmax];
u[[1, n]] = posT[time];
, {n, 1, k}
];
Do[
u[[2, n]] =
0.5 ((1 - γ/2) u[[1, n - 1]] + (1 + γ/2) u[[1,
n + 1]] + 2 Δx/EA F[[n]] +
g Δt^2);
, {n, 2, k - 1}
];
Do[
Do[
u[[i + 1, n]] = (1 + γ/2) u[[i, n + 1]] -
u[[i - 1, n]] + (1 - γ/2) u[[i, n - 1]] +
g Δt^2;
, {n, i + 1, k - i}
];
, {i, 2, n2 + 1}
];
If[m < nrod,
timez = Table[0, {ii, k}];
posz = Table[0, {ii, k}];
];
Do[
F[[n]] =
EA/Δx/
2 (2 u[[n2 + 2, n]] - (1 - γ/2) u[[n2 + 1,
n - 1]] - (1 + γ/2) u[[n2 + 1, n + 1]] -
g Δt^2);
t[[n]] = (n - j/2 - 1) Δt;
If[m < nrod && t[[n]] >= 0 && t[[n]] <= tmax + 0.1,
loadz[[n - n2 - 1]] = F[[n]] + buoy;
timez[[n - n2 - 1]] = t[[n]];
posz[[n - n2 - 1]] = u[[n2 + 1, n]];
];
If[t[[n]] < 0, t[[n]] = t[[n]] + tmax];
If[t[[n]] > tmax, t[[n]] = t[[n]] - tmax];
, {n, n2 + 2, k - n2 - 1}
];
timez = DeleteCases[timez, 0];
posz = DeleteCases[posz, 0];
posT = Interpolation[
Table[{timez[[n]], posz[[n]]}, {n, Length[timez]}]];
nn = k - 2 (n2 + 1);
n2 = n1;
, {m, 1, 3}
];


Now we have positions and load conditions at the bottom of the rod string where the pump is. Make new interpolation functions so we can plot the bottom conditions.

tf = Table[t[[n]], {n, j/2 + 1, k - j/2, 2}];
loadf = Table[(F[[n - 3]] + F[[n - 2]] + F[[n - 1]] + F[[n]] +
F[[n + 1]] + F[[n + 2]] + F[[n + 3]])/7, {n, j/2 + 1, k - j/2,
2}];
posf = Table[u[[n3 + 1, n]], {n, j/2 + 1, k - j/2, 2}];
posT = Interpolation[Table[{tf[[n]], posf[[n]]}, {n, Length[tf]}]];

AspectRatio -> 1/GoldenRatio, AxesLabel -> {"Position", "Load"}]


The ideal bottom hole plot is a perfect rectangle, and this well is in pretty good shape. Again, the jerky plot is typical of real data. I smoothed the loads by taking a 7 point average. The damping coefficient in general is not known with great accuracy. Fortunately, the shape of the bottom hole curve is relatively insensitive to the damping coefficient, and it is the shape, that determines well problems. It is much more difficult to diagnose well problems from the surface data than it is with bottom hole data.

Again, in time the pump moves clockwise, the higher loads occur as the pump rises. This routine also takes into account the change in buoyancy force moving to a smaller diameter rod string.

This procedure is adapted from a FORTRAN program I wrote in about 1982. I have adapted it into Mathematica code, but some purists may think it is not very good Mathematica code and they will be right. This program is very fast on modern computers and I am not about to spend a bunch of hours messing with it. Be glad you don't have to run it on a 2 Mhz 286 computer without a math coprocessor.

• Er… what's the i.c. and b.c. imposed here? I can spot Derivative[1, 0][u][0, t] == F[t]/EA, but can't distinguish others… Jul 4, 2020 at 8:47
• @xzczd The b.c are the polish rod position and load vs time at the surface which are specified in the data array. My first plot shows the b.c. The force equation you mention is just the F finite difference equation with i = 1. This is a steady-state situation, so we do not have a true i.c. We are not turning on the pumping unit at a time 0 and then computing its motion. We have take the first point in the data array to be at time 0, but with steady-state periodic motion, we could have taken time 0 to be any point on the position-load curve. Jul 4, 2020 at 17:55
• First, thank you to everyone for their feedback and comments, they are all very helpful. @BillWatts Thank you very much as this is exactly what I am looking for. A lot to process for me but this is fantastic. Like a kid at Christmas right now reading the comments! Jul 6, 2020 at 16:24
• I had a chance to review the code and play with my own test data It works very well, although I'm not 100% sure what J=40 represents. Thank you @BillWatts! Greatly appreciated and great answer! Jul 12, 2020 at 13:05
• @wpilgri I am glad its working out for you. j does not have to be 40, but because the finite difference equations have n+1 and n-1 expressions,j allows extra terms for times beyond 0 and tmax in the finite difference equations. We throw them out when we are done, but they need to be there for the calculation. Also, I have found a wide range of values for the speed of sound in various types of steel, so you can fool with that. The value I used was on the higher end. I guess you need to be in the oil business to understand what your question was asking. Jul 12, 2020 at 17:37

As recommended by PaulCommentary and xzczd we put bc = DirichletCondition[z[x, t] == 0, x == 1]; to the end x=1 and apply force to the end x=0. Also we use physical normalization and change phase of force in accordance with zero initial condition, then we have

Y = 199*^9;(*Pa*)ρ = 7860;(*kg/m^3*)dia =
1/39.37;(*1" dia converted to meters*)c = Sqrt[Y/ρ];
g = 9.81;
area = π*dia^2/4;
endPrescribedForce[t_] :=
100*Sin[t];(*Prescribed Force BC*)(*prescribed Force BC.Translated \
to strain for NeumannValue BC*)eq1 =
D[z[x, t], {t, 2}] - c^2*D[z[x, t], {x, 2}] - g ==
NeumannValue[endPrescribedForce[t]/(ρ*area), x == 0];
bc = DirichletCondition[z[x, t] == 0, x == 1];

sol = NDSolve[{eq1, bc, z[x, 0] == 0},
z[x, t], {x, 0, 1}, {t, 0, 2*π}];
solf[x_, t_] := sol[[1, 1, 2]]
strain[x_, t_] :=
Evaluate[D[solf[x, t],
x]] (*Take derivative of solution to get strain*)

{Plot3D[solf[x, t], {x, 0, 1}, {t, 0, 2*π},
PlotLabel -> "Displacement", Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic, Boxed -> False],
Plot[strain[0, t], {t, 0, 2*π},
PlotLabel -> "Strain at Prescribed End"],
Plot[-strain[0, t]*c^2 area ρ, {t, 0, 2*π},
PlotLabel -> "Calculated Force at Prescribed End"]}


Update 1. We can improve code by using "MethodOfLines". Then we can calculate displacement for arbitrary time (not only up to $$2\pi$$):

Y = 199*^9;(*Pa*)ρ = 7860;(*kg/m^3*)dia =
1/39.37;(*1" dia converted to meters*)c = Sqrt[Y/ρ];
g = 9.81;
area = π*dia^2/4;
endPrescribedForce[t_] :=
100 Sin[t];(*Prescribed Force BC*)(*prescribed Force BC.Translated \
to strain for NeumannValue BC*)eq1 =
D[z[x, t], {t, 2}]/c^2 - D[z[x, t], {x, 2}] - g/c^2 ==
NeumannValue[endPrescribedForce[t]/(ρ*area), x == 0]/c^2
bc = DirichletCondition[z[x, t] == 0, x == 1];
reg = Line[{{0}, {1}}];
tm = 7.1;
sol = NDSolve[{eq1, z[x, 0] == 0, Derivative[0, 1][z][x, 0] == 0,
z[1, t] == 0}, z[x, t], {t, 0, tm}, {x} ∈ reg,
Method -> {"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 5},
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {z -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0002}}}},
MaxStepSize -> 0.1];
solf[x_, t_] := sol[[1, 1, 2]]
strain[x_, t_] :=
Evaluate[D[solf[x, t],
x]] (*Take derivative of solution to get strain*)

{Plot3D[solf[x, t], {x, 0, 1}, {t, 0, tm},
PlotLabel -> "Displacement", Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic, Boxed -> False],
Plot[strain[0, t], {t, 0, tm},
PlotLabel -> "Strain at Prescribed End", PlotRange -> All],
Plot[-strain[0, t]*c^2 area ρ, {t, 0, tm},
PlotLabel -> "Calculated Force at Prescribed End", PlotRange -> All]}


• Well, this does give a solution, but without adding another initial condition, NDSolve will automatically add zero Neumann value at the end time (in this case $t=2 \pi$) and solve the problem as a pure BVP, while pure BVP of wave equation is a well-known ill-posed problem. Jun 30, 2020 at 2:51
• @xzczd It is not look so bad as with unphysical condition bc = DirichletCondition[z[x, t] == endPrescribedDisp[t], x == 0]. Also this solution gives right force scale. We can improve it by adding initial condition. Jun 30, 2020 at 11:15
• Partly because c is large in this case. The following is an example in which the instablity is more obvious: test[tend_] := NDSolveValue[{D[z[x, t], {t, 2}] - D[z[x, t], {x, 2}] - 10 == 0, z[x, 0] == 0, z[0, t] == 0, z[1, t] == 0, z[x, tend] == 0}, z, {x, 0, 1}, {t, 0, tend}, Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.001}]; Plot[{test[2 Pi][1/2, t], test[1.01 2 Pi][1/2, t]} // Evaluate, {t, 0, tend}] As we can see, we've only varied the value of tend for $1\%$, but the solution changes drastically. Jun 30, 2020 at 12:00
• @xzczd See Update 1 to my answer. Jun 30, 2020 at 18:11
• Well, now you've imposed 3 b.c.s in $x$ direction. A quick test shows that NDSolve seems to simply ignore the b.c. NeumannValue[-Derivative[0, 1][z][x, t], x == 1]. Jul 1, 2020 at 2:11