# Compute the curvature of a bent stack of paper

My goal is to link the curvature of a bent stack of paper in a wind tunnel to its bending modulus $$B$$, knowing all the other physical properties. To this end, I would like to reproduce the numerical simulation in P. Buchak, C. Eloy and P. M. Reis, Phys. Rev. Lett. 105, 194301 (2010), also available here.

With $$s$$ the curvilinear abscissa along the stack, and $$\theta(s)$$ the angle between the stack and horizontal at $$s$$, the curvature of the stack should be governed by:

Fx'[s] ==
-rhoa*U^2*W*Sin[theta[s]]*(N[Pi,10]/4*W*Cos[theta[s]]^2*theta'[s] -
0.5*Cd*Sin[theta[s]]^2)

Fy'[s] ==
-rhoa*U^2*W*Cos[theta[s]]*(0.5*Cd*Sin[theta[s]]^2 -
N[Pi,10]/4*W*Cos[theta[s]]^2*theta'[s])+sigmap*9.81*W

B*W*theta''[s] == Sin[theta[s]]*Fx[s] - Cos[theta[s]]*Fy[s]


with the following boundary conditions:

Fx[L] == 0,Fy[L] == 0,theta[0] == 0,theta'[L] == 0


where $$L$$ is the length of the stack.

The other parameters are: $$U$$ (wind speed), $$\rho_a$$ (density of air), $$W$$ (stack width), $$C_d$$ (drag coefficient), $$\sigma_p$$ (area density of stack) — all of which I know experimentally — and $$B$$ (bending modulus) which I am looking for.

I tried to compute these equations using Mathematica's NDSolve.

sol =
NDSolve[{(*equations*),(*boundary conditions*)}, {Fx, Fy, theta}{s, 0, L},
MaxSteps -> Infinity, InterpolationOrder-> All];


and plotted the expected shape of the stack.

The result is that... it doesn't match at all! (the following figures correspond to $$U=5m\cdot s^{-1}$$)

Furthermore, the code can lead to very surprising results for high wind speed! (the following figure corresponds to $$U = 10m\cdot s^{-1}$$ and $$L = 0.3m$$, all the other parameters being unchanged)

Are there other computational methods which would work to solve this ODE system, or more efficient and accurate options to use with NDSolve?

-------- EDIT --------

As requested, I provide here a complete Mathematica code.

The governing equations of my system can be rewritten as

Fx'[s] == Sin[theta[s]](0.5*Cd*Cy*Sin[theta[s]]^2-
0.25*N[Pi,10]*lambda*Cy*Cos[theta[s]]^2*theta'[s])
Fy'[s] == -Cos[theta[s]](0.5*Cd*Cy*Sin[theta[s]]^2-
0.25*N[Pi,10]*lambda*Cy*Cos[theta[s]]^2*theta'[s])+Delta
theta''[s] == Sin[theta[s]]*Fx[s]-Cos[theta[s]]*Fy[s]


and have the same boundary conditions as previously.

The algorithm now has 4 adimensionned parameters: $$C_d = 1.76$$ (drag coefficient), $$C_y = \frac{\rho_aU^2L^3}{B}$$ (Cauchy number, $$C_y > 18.5$$), $$\lambda = \frac{W}{L}$$ (aspect ratio, $$0.1 < \lambda < 0.5$$) and $$\Delta = \frac{\sigma_pgL^3}{B}$$ (elastrogravity number, $$0.1 < \Delta < 20$$).

(* PARAMETERS *)
Cd = 1.76; (* drag coefficient *)
Cy = 90; (* Cauchy number *)
lambda = 0.4; (* aspect ratio *)
Delta = 1; (* elastogravity number *)

(* SIMULATION *)
sol = NDSolve[{Fx'[s]==Sin[theta[s]](0.5*Cd*Cy*Sin[theta[s]]^2-
0.25*N[Pi,10]*lambda*Cy*Cos[theta[s]]^2*theta'[s]),
Fy'[s]==-Cos[theta[s]](0.5*Cd*Cy*Sin[theta[s]]^2-
0.25*N[Pi,10]*lambda*Cy*Cos[theta[s]]^2*theta'[s])+Delta,
theta''[s]==Sin[theta[s]]*Fx[s]-Cos[theta[s]]*Fy[s],
Fx[1]==0,Fy[1]==0,theta[0]==0,theta'[1]==0},{Fx,Fy,theta},
{s,0,1},MaxSteps->Infinity,InterpolationOrder->All];
(* résolution du système *)

(* PLOT *)
ParametricPlot[{Integrate[Cos[Evaluate[theta[t] /. sol][[1]]], {t,0,s}],
Integrate[Sin[Evaluate[theta[t] /. sol][[1]]],{t,0,s}]},
{s, 0, 1},AspectRatio->1,AxesLabel->{"abscissa","ordinate"},
PlotLabel->"Computed shape"]


The (* PLOT *) section returns a graph of the computed shape of the stack (ie, $$y = \int_0^1 ds \sin(\theta)$$ vs. $$x = \int_0^1 ds \cos(\theta)$$).

Executing this script gave me the following graph

which does not match the experimental observations (see 1st figure of the question).

• Please give all the numerical values for your parameters in your code, and show the exact code snipped that you used to generate those plots. Feb 24, 2017 at 23:55
• I'm going to upvote this question and at the same time vote to put this on hold as long as you don't provide the additional numerical values for the parameters. It is a very interesting question, although I'm not sure someone can help you without studying the paper (which is too much for the common user). Nevertheless, except for the missing parameters, your questions is very well written and extremely interesting. Feb 25, 2017 at 4:09
• The numerical values I used to test the code where the ones provided on figure 1 (from the paper): rhoa=1.225[kg.m-3], W=0.029[m], Cd=1.76, sigmap=0.0391[kg.m-2], B=0.0082759[N.m], L=0.15 to 0.3[m], U=5 to 10[m.s-1]. Thanks for your interest for my question! Feb 25, 2017 at 5:46
• @yketa you should edit your question and include are the relevant information.
– zhk
Feb 25, 2017 at 6:39
• This is exactly what I have done. My equations are not linear, therefore there might be different solutions for the same boundary conditions. NDSolve may always return the same solution, but is there a way to look for others? Feb 27, 2017 at 8:07

yketa's guess is right. This nonlinear boundary value problem (BVP) has multiple solutions, and we need to manually set proper initial guess for "Shooting" method to get the desired result. Let me reproduce FIG.2 in the paper:

(* Parameters are taken from the paper. *)
Ulst = {111/10, 65/10, 88/10};
g = 981/100;
rhoa = 129/100;
L = 15 10^-2;
W = 29/10 10^-2;
BW = 24 10^-5;
rhophW = 114/1000/10^3 10^2;
rhoph = rhophW/W;
B = BW/W;
Delta = rhoph g L^3/B;
lambda = W/L;
Cd = 176/100;
Cy = rhoa U^2 L^3/B;

k = 1/2 Cd Cy Sin[theta[s]]^2 - Pi/4 lambda Cy Cos[theta[s]]^2 theta'[s];

(* The integration can be computed inside NDSolve, too. *)
add = {intSin'[s] == Sin@theta[s], intSin[0] == 0, intCos'[s] == Cos@theta[s],
intCos[0] == 0};

Clear[sol, U]
shoot[ic_] := {"Shooting", "StartingInitialConditions" -> ic};
sol[var_] := Block[{U = var}, First@NDSolve[{
Fx'[s] == Sin[theta[s]] k,
Fy'[s] == -Cos[theta[s]] k + Delta,
theta''[s] == Sin[theta[s]] Fx[s] - Cos[theta[s]] Fy[s],
Fx[1] == 0, Fy[1] == 0, theta[0] == 0, theta'[1] == 0, add},
{Fx, Fy, theta, intSin, intCos}, {s, 0, 1},
Method -> shoot@{Fx[1] == 0, Fy[1] == 0, theta[1] == 3,
theta'[1] == 0, intSin[1] == 0, intCos[1] == 0}]];

ParametricPlot[{intCos@s, intSin@s} /. sol /@ Ulst // Evaluate, {s, 0, 1},
AxesLabel -> {"abscissa", "ordinate"}, PlotLabel -> "Computed shape"]


In the code piece above, I simply dope out the proper initial guess intuitively, but this can actually be done in a systematic way for your problem. Let me show you the way with the parameters given in your question as an example.

We first set up a initial value problem with ParametricNDSolveValue:

Cd = 1.76;(*drag coefficient*)
Cy = 90;(*Cauchy number*)
lambda = 0.4;(*aspect ratio*)
Delta = 1;(*elastogravity number*)
k = 1/2 Cd Cy Sin[theta[s]]^2 - Pi/4 lambda Cy Cos[theta[s]]^2 theta'[s];
add = {intSin'[s] == Sin@theta[s], intSin[0] == 0, intCos'[s] == Cos@theta[s],
intCos[0] == 0};
test = ParametricNDSolveValue[{
Fx'[s] == Sin[theta[s]] k,
Fy'[s] == -Cos[theta[s]] k + Delta,
theta''[s] == Sin[theta[s]] Fx[s] - Cos[theta[s]] Fy[s],
Fx[1] == 0, Fy[1] == 0,(*,theta[0]\[Equal]0*)theta[1] == par,
theta'[1] == 0}, {Fx, Fy, theta}, {s, 0, 1}, par];


Then we need to find out all the possible values of par that result in theta[0] == 0. (It is not hard for your specific problem because 0 < par <2 Pi according to the model. ) This job can be done by e.g. RootSearch, but here I'll just do it in a more basic way i.e. by observing the graph:

Plot[test[par][[3]][0], {par, 0, 2 Pi}, PlotRange -> 1]


From the graph we can see 3 (4?) possible solutions. We try them from right to left:

rst = NDSolve[{
Fx'[s] == Sin[theta[s]] k,
Fy'[s] == -Cos[theta[s]] k + Delta,
theta''[s] == Sin[theta[s]] Fx[s] - Cos[theta[s]] Fy[s],
Fx[1] == 0, Fy[1] == 0, theta[0] == 0, theta'[1] == 0, add}, {Fx, Fy,
theta, intSin, intCos}, {s, 0, 1},
Method -> shoot@{Fx[1] == 0, Fy[1] == 0, theta[1] == 2.6,
theta'[1] == 0, intSin[1] == 1, intCos[1] == 1}];

ParametricPlot[{intCos@s, intSin@s} /. rst // Evaluate, {s, 0, 1},
AxesLabel -> {"abscissa", "ordinate"}, PlotLabel -> "Computed shape"]


OK, the right-most root is the one we need.

• Thank you so very much for your answer and your time! All seems a lot clearer now! Feb 27, 2017 at 15:12
• @yketa You're welcome. I wouldn't have found out the answer so quickly if you hadn't pointed out the possibility of existence of multiple solutions :) Feb 27, 2017 at 16:25