I am a student of mathematics and for my modelling seminar, I need to dip a little bit into engineering in order to model the movements of a suspension bridge under influence of wind. This will be modeled as a one-dimensional vibrating beam and so I am trying to get a grasp on the following partial differential equation:

$$ \rho \frac{\partial^2}{\partial t^2} v(x,t) = -EI \frac{\partial^4}{\partial x^4} v(x,t) -\delta \frac{\partial}{\partial t}v(x,t) +q(x,t)$$

with boundary condition

$$ v(0,t) = v(L,t) = \frac{\partial^2}{\partial x^2} v(0,t) = \frac{\partial^2}{\partial x^2} v(L,t) = 0 $$

Here, $v(x,t)$ is supposed to give the downward deflection at time $t$ of the point at length $x$ in the beam, subject to the load $q(x,t)$. $E$ is Young's modulus, $I$ the second moment of inertia of the beam cross-section, $\rho$ is the mass per unit length and $\delta$ some small dampening coefficient.

I have tried to approximate a solution to this in MATLAB by discreticizing the spatial variable (representing the fourth derivative with finite differences) and introducing a dummy variable for the first derivative in time, thereby transforming the whole thing into an ordinary differential equation in time.

However, my problem is that solutions $v$ seem to exhibit extremely high-frequency oscillations that make the ode45 function take forever to finish and let the resultant graph look almost like a solid block of color even on a time scale of [0 10]. Is a Euler beam supposed to oscillate quite as much? I'm really having trouble even obtaining any sort of useable result.

This is especially the case when trying to insert realistic values for $E$, i.e. at least $10^9$. (I notice though that many numerical simulations across the internet use very small values (<=100) -what might be the justification for that?)

P.S. Apologies if any of this is badly or insufficiently phrased - again, as a mathematician I am a little out of my depth here.


Note that for calculation I have actually divided the equation by $\rho$ on both sides, so we get

$$ \frac{\partial^2}{\partial t^2} v(x,t) = -\frac{EI}{\rho} \frac{\partial^4}{\partial x^4}v(x,t) - \delta' \frac{\partial}{\partial t} v(x,t) + g(x,t) $$

where $g$ is now the downward acceleration (only due to gravity in this case). For actual values, I want to take:

  • $E = 2 \cdot 10^{11} N/m^2$
  • $I = 20 m^4$. This comes from the formula $ I = \frac{bh^3}{12} $ where $b = 30m$ is the width and $h = 2m$ the height of the beam cross-section.
  • $\rho = 2 \cdot 10^5 kg/m$
  • $\delta' = 0.01 kg/ms$
  • $g(x,t) = 10 m/s^2$ constant, obviously a rough approximation of acceleration due to gravity
  • $L = 2000 m$

However, the combined parameter $a := \frac{EI}{\rho}$ which would thus end up as $a = 2 \cdot 10^7 m^4/s^2$ and I am not getting the ODE solver to terminate unless I take something as low as $a = 1 m^4/s^2$ (still takes several minutes).

  • $\begingroup$ Could you add exact length, cross-section dimensions, $E$ and loads? (with units, not just numbers) $\endgroup$ May 29, 2023 at 18:07
  • $\begingroup$ @TomášLétal Of course, I have edited my question accordingly. Please tell me if you need more info, e.g. on the discretization scheme. $\endgroup$
    – Eriol
    May 30, 2023 at 14:02
  • $\begingroup$ It seems like you are using full section, which would be very inefficient and expensive for actual bridges. They have some structure to save material resulting in lower mass per length. Check out "box-girder" bridge. $\endgroup$ May 30, 2023 at 19:01


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.