Catenary in Wind

The catenary problem but with a twist

\[ \newcommand{\unit}[1]{\,\mathrm{#1}} \DeclareMathOperator{\D}{d} \]

I came across this problem at work when we were prototyping a tethered drone system to supply in principle endless power to a multirotor drone. This can be used in natural catastrophes or military areas to cover it with a communication network. The problem described here is not really necessary for the prototyping so I took it as a personal project.

I had two main questions that I was constantly asking myself:

1. Can I derive the shape of the cable?
2. How much extra downward force is added to the drone and how does it change with wind?

Derivation

We look at the cable that has rigid end points \((0, 0)\) an \((x_0, y_0)\). There are two external forces acting on the cable, wind \(F_w\) and gravity \(F_g\). Figure 1 is depicting the situation. On the right there is a closeup of a small segment \(\D s\) of the cable. The forces acting on this differential segment are shown in blue.

Figure 1: Forces Acting on a Hanging Cable

We write the force balance for both \(x\) and \(y\)-axis \[ \begin{align*} &x: &\quad T_x(x+\D x, y+\D y) - T_x(x, y) &= \rho g \D s \\ &y: &\quad -T_y(x+\D x, y+\D y) + T_y(x, y) &= f_w \D x. \end{align*} \tag{1}\] Here we separate the forces into perpendicular components \(T_x = T \cos\alpha\) and \(T_y = T \sin\alpha\). The parameters are mass density per unit length \(\rho\) and wind pressure \(f_w = F_w / l\), where \(l\) is the total length of the rope.

Naive approach

We seek for the solution \(y(x)\), exactly one \(y\) for each \(x\). From the Pythagoras we have \[ \D s = \D x \sqrt{ 1 + \left( \frac{\D y}{\D x} \right)^2 }. \] Notice that when we divide both equations of system Equation 1 by \(\D x\) we have the definition of the derivative on the left side \[ \begin{align*} \frac{\D T_x}{\D x} &= \rho g \sqrt{ 1 + \left( \frac{\D y}{\D x} \right)^2 } \\ \frac{\D T_y}{\D x} &= -f_w. \end{align*} \] The second equation is easily solved and the solution is \[ T_y = -f_w (x - \bar{x}), \] where \(\bar{x}\) is the integration constant. The solution is a linear function that crosses zero at \(\bar{x}\). Since in the rope the forces always have to be tangential to it, the only place where \(y\) component will be zero is the extremal point. From the defintions of the force components we have that \(T_x = T_y/\tan\alpha = T_y / (dy/dx)\). Using the chain rule the first equation becomes \[ \frac{\D T_y}{\D x} \frac{\D y}{\D x} - T_y \frac{\D^2y}{\D x^2} = \rho g \sqrt{ 1 + \left( \frac{\D y}{\D x} \right)^2 } \left( \frac{\D y}{\D x} \right)^2. \] Since we know \(T_y\) we can replace it and have our final equation \[ \begin{gathered} -\frac{\D y}{\D x} + (x - \bar{x}) \frac{\D^2y}{\D x^2} = \kappa \sqrt{ 1 + \left( \frac{\D y}{\D x} \right)^2 } \left( \frac{\D y}{\D x} \right)^2, \\ l = \int_0^{x_0} \D s. \end{gathered} \tag{2}\] where \(\kappa = \rho g / f_w\) is a constant. This is a pretty complicated differential equation! And notice that we still do not know \(\bar{x}\) but it has a very clear physical meaning. So in order to get a solvable system we need the length constraint.

Problems

Solving this equation was a nightmare! I tried every possible method I knew for boundary value problems.

I started by using finite differences. The optimiser optimistix.IndirectLevenbergMarquardt worked reasonably well, but I had to lower the tolerances to \(10^{-3}\). It works for some very specific parameters but the final solution depends heavily on the initial guess. I used a second order polynomial with roots at the boundaries and the solution differs when scaling the initial guess.

To my surprise, shooting method worked better in the sense that it found only one correct solution. But it worked in an even smaller region of parameters. The reason I was surprised was because if you write the Equation 2 in the form of first order differential equation system standard form, then it involves division by \((x-\bar{x})\) which includes division by zero at \(\bar{x}\). But the solver handled it nicely if I added \(\varepsilon=10^{-4}\) to both the nominator and denominator. From the physics we know that \(\D y/\D x=0\) at \(\bar{x}\).

Parametric equation

Since there were a few problems with Equation 2, I needed to approach the problem from a different angle. So instead of solving for \(y(x)\), now we solve for \((x(s), y(s))\), where \(s\) is the cable length. This is beneficial because with high wind speed or very loose cable, the lower end will naturally want to drop down which would mean there are two \(y\)-s for each \(x\).

Instead of expressing the derivatives in \(x\) we now divide the force balance Equation 1 by \(\D s\) \[ \begin{align} \frac{\D T_x}{\D s} &= \rho g \\ \frac{\D T_y}{\D s} &= f_w \frac{\D x}{\D s}. \end{align} \] From Figure 1 one can easily derive the geometric constraints for the Cartesian coordinates \[ \frac{\D x}{\D s} = \cos{\alpha} \qquad \frac{dy}{ds} = \sin{\alpha}. \] Since the tension in the rope and the vector \((\D x, \D y)\) both have to be tangential to the rope, a similar relationship holds for tension \[ \cos{\alpha} = \frac{T_x}{\sqrt{T_x^2 + T_y^2}} \qquad \sin{\alpha} = \frac{T_y}{\sqrt{T_x^2 + T_y^2}}. \] So putting it all together we get a system of four equations \[ \begin{aligned} \frac{\D x}{\D s} &= \frac{T_x}{T} & \frac{\D y}{\D s} &= \frac{T_y}{T} \\ \frac{\D T_x}{\D s} &= \rho g & \frac{\D T_y}{\D s} &= f_w \frac{T_x}{T}. \end{aligned} \tag{3}\] This system looks much easier than the one for finding \(y(x)\) Equation 2.

Wind force

For aerodynamic drag we can use a simple equation \[ F_w = \frac 12 \rho v^2 C_d S, \] where \(\rho\) is the air density, \(v\) wind speed, \(C_d\) drag coefficient and \(S\) the cross-sectional area. Knowing that the cable has diameter \(d=1.2\unit{mm}\) and drag coefficient of a cylinder is \(C_d = 1.17\) the equation for wind pressure is \[ f_w = \frac{C_d}{2} \rho v^2 d. \]

Solution

The problem has 4 natural boundary conditions \[ \begin{aligned} x(0) &= 0 & x(L) &= x_0 \\ y(0) &= 0 & y(L) &= y_0. \end{aligned} \]

I used Jax to solve this with a shooting method.

# parametric_shooting.py
import jax.numpy as jnp
import diffrax
import equinox as eqx
from jaxtyping import Array, Float64, Int
import optimistix as optx
import matplotlib.pyplot as plt


class EquationParameters(eqx.Module):
    n: Int
    s: Float64[Array, "n"]  # noqa: F821
    end_point: Float64[Array, "2"]
    length: Float64
    mass: Float64
    wind_pressure: Float64

    @property
    def g(self) -> Float64:
        return 9.81

    def __init__(self, n, end_point, length, mass, wind_pressure):
        self.n = n
        self.s = jnp.linspace(0.0, length, num=n)
        self.end_point = end_point
        self.length = length
        self.mass = mass
        self.wind_pressure = wind_pressure


class InitialCondition(eqx.Module):
    """Initial conditions for solving the system of equations. Values at s=0"""

    x: Float64
    y: Float64
    Tx: Float64
    Ty: Float64


def system(_s, S, params: EquationParameters):
    _x, _y, Tx, Ty = S

    T = jnp.sqrt(Tx**2 + Ty**2)
    cos = Tx / T
    sin = Ty / T

    return jnp.array([cos, sin, params.mass * params.g, params.wind_pressure * cos])


def solve_ivp(ic: InitialCondition, params: EquationParameters):
    term = diffrax.ODETerm(system)
    solver = diffrax.Tsit5()
    controller = diffrax.PIDController(
        rtol=1e-6,
        atol=1e-6,
    )
    y_init = jnp.array([ic.x, ic.y, ic.Tx, ic.Ty])
    saveat = diffrax.SaveAt(ts=params.s)
    sol = diffrax.diffeqsolve(
        term,
        solver,
        t0=params.s[0],
        t1=params.s[-1],
        dt0=0.1,
        y0=y_init,
        args=params,
        stepsize_controller=controller,
        saveat=saveat,
        adjoint=diffrax.DirectAdjoint(),
    )
    return sol


def shooting_objective(tension: Float64[Array, "2"], params: EquationParameters):
    ic = InitialCondition(0.0, 0.0, tension[0], tension[1])
    sol = solve_ivp(ic, params)
    assert sol.ys is not None, "Solver did not return any values"
    sol_end = sol.ys[-1, [0, 1]]
    return sol_end - params.end_point


def solve(tension_guess: Float64[Array, "2"], params: EquationParameters):
    solver = optx.IndirectLevenbergMarquardt(
        rtol=1e-5,
        atol=1e-5,
        verbose=frozenset({"step", "accepted", "loss", "step_size"}),
    )

    sol_root = optx.root_find(
        shooting_objective,
        solver,
        tension_guess,
        args=params,
        max_steps=1000,
        throw=False,
    )

    best_tension = sol_root.value
    ic = InitialCondition(0.0, 0.0, best_tension[0], best_tension[1])
    final_sol = solve_ivp(ic, params)
    return final_sol

The example cable tensions and shape is on Figure 2.

# example_sol.py
import jax.numpy as jnp
import matplotlib.pyplot as plt
from parametric_shooting import *

fw = (9**2) / 2 * 1.225 * 1.17 * (1.5 / 1000)
params = EquationParameters(200, jnp.array([100.0, 0.0]), 120, 0.004, fw)
sol = solve(jnp.array([1.0, -1.0]), params)
assert sol.ys is not None, "Solver did not return any values"

x = sol.ys[:, 0]

fig, ax = plt.subplots(2, 1, sharex=True)
ax[0].plot(x, sol.ys[:, 1])
ax[0].set_ylabel("y [m]")

ax[1].plot(x, sol.ys[:, 2], label="Tx")
ax[1].plot(x, sol.ys[:, 3], label="Ty")
ax[1].set_ylabel("Tension [N]")
ax[1].set_xlabel("x [m]")
ax[1].legend()
plt.savefig("cable_100-120.png")
plt.show()

Figure 2: Example solution for wind speed of \(9 \unit{m/s}\), cable length \(120 \unit{m}\) and height of \(100 \unit{m}\).

Figure 3 shows that we would like to keep like ~30% slack on the cable to minimise the forces acting on the drone.

# length-tension.py
import jax.numpy as jnp
import jax
from jaxtyping import Float
import matplotlib.pyplot as plt
from parametric_shooting import EquationParameters, solve


@jax.jit
def tension(L: Float, fw: Float) -> Float:
    params = EquationParameters(200, jnp.array([100.0, 0.0]), L, 0.004, fw)
    sol = solve(jnp.array([1.0, -1.0]), params)
    assert sol.ys is not None, "Solver did not return any values"

    T_end = jnp.sqrt(sol.ys[-1, 2] ** 2 + sol.ys[-1, 3] ** 2)
    return T_end


fw = (9**2) / 2 * 1.225 * 1.17 * (1.5 / 1000)
lengths = jnp.linspace(101, 150, 20)

tension_vmap = jax.vmap(tension, (0, None))
T = tension_vmap(lengths, fw)

plt.plot(lengths, T)
plt.xlabel("length [m]")
plt.ylabel("Tension [N]")
plt.savefig("length-tension.png")
plt.show()

Figure 3: Force acting on the drone depending on the slack with wind speed of \(9 \unit{m/s}\) and height of \(100 \unit{m}\).