Center of MechanicsStephan Kaufmann


Example: Bending Waves in a Rod


Let us consider bending waves in an isotropic rod (with constant cross-section).

We use the following notation:
s: components of stress tensor,
u: components of displacements,
e: modulus of elasticity,
n: Poisson's ratio.

The x-axis is chosen in direction of the rod, the y,z-axes perpendicular.


We first construct lists for the components of s and u.

  sFcts =
    Flatten[Table[s[j, i][x, y, z, t], {i, 3}, {j, i}]] /.
    Thread[Range[3] -> {x, y, z}]

  uFcts = Table[u[i][x, y, z, t], {i, 3}] /.
    Thread[Range[3] -> {x, y, z}]

  fcts = Join[sFcts, uFcts];


The equations of motion are:

  em =
      D[s[x,z][x,y,z,t],z] == r D[u[x][x,y,z,t],{t,2}],
      D[s[y,z][x,y,z,t],z] == r D[u[y][x,y,z,t],{t,2}],
      D[s[z,z][x,y,z,t],z] == r D[u[z][x,y,z,t],{t,2}]

Together with the kinematic relations and the equation of linear elasticity

  krAndEl =
    e D[u[x][x,y,z,t],x] == s[x,x][x,y,z,t]-
      n (s[y,y][x,y,z,t]+s[z,z][x,y,z,t]),
    e (D[u[x][x,y,z,t],y]+D[u[y][x,y,z,t],x]) ==
      2 (1+n) s[x,y][x,y,z,t],
    e (D[u[x][x,y,z,t],z]+D[u[z][x,y,z,t],x]) ==
      2 (1+n) s[x,z][x,y,z,t],
    e (D[u[y][x,y,z,t],z]+D[u[z][x,y,z,t],y]) ==
      2 (1+n) s[y,z][x,y,z,t],
    2 e D[u[y][x,y,z,t],y] == 2(1+n)s[y,y][x,y,z,t] -
    2 e D[u[z][x,y,z,t],z] == 2(1+n)s[z,z][x,y,z,t] -

we get the following system of PDEs:

  physicalEquations = Join[em, krAndEl];

  IndexForm[physicalEquations, Functions -> fcts]


Let us rename the functions and the variables such that we can later reuse the old names for the dimension-less functions and variables.

  sFcts1 = sFcts /.
    f_[x, y, z, t] -> f[x1, y1, z1, t1] /.
    s -> s1

  uFcts1 = uFcts /. 
     f_[x, y, z, t] -> f[x1, y1, z1, t1] /. 
     u -> u1

  fcts1 = Join[sFcts1, uFcts1];

We set up equations for scaling the stresses with respect to a reference stress s0.

  scaleS = Thread[sFcts1 == sFcts / s0]

The reference displacement shall be u0.

  scaleU = Thread[uFcts1 == uFcts / u0]

  scaleFcts = Join[scaleS, scaleU];

The x-axis (in direction of the rod) is scaled by the length of the rod l0. For the perpendicular axes y and z, we choose the height h0 of the rod. The time is scaled by a reference time t0 (to be chosen later).

  scaleVars = {x1 == x/l0, y1 == y/h0, z1 == z/h0,
    t1 == t/t0}

We can now transform the equations to the new functions and variables. To keep the names short, we immediately reuse the old names for the scaled functions.

  scaledEq1 = Transform[physicalEquations,
    fcts, fcts1, scaleFcts, scaleVars] /.
    {s1 -> s, u1 -> u, x1 -> x, y1 -> y, z1 -> z,
     t1 -> t};

  IndexForm[scaledEq1, Functions -> fcts]

Replacing Parameters

We now choose the parameter according to physical considerations.

The reference stress is set to e u0/h0.

  (scaledEq2 = SimplifyEquation[scaledEq1, fcts,
    s0 == e u0/h0, s0]) // IndexForm

A new parameter b is set to r/e (h0/t0)2.

  (scaledEq3 = SimplifyEquation[scaledEq2, fcts,
    b == r/e (h0/t0)^2, {h0, t0, r, e}]) // IndexForm

And finally, the small parameter eps is chosen to be h0/l0 (length/height of the rod).

  (scaledEq = SimplifyEquation[scaledEq3, fcts,
    eps == h0/l0, {h0, l0}]) // IndexForm

Boundary and Initial Conditions

The initial displacements uyt0, uzt0 and their time derivatives uty0p, utz0p, together with the y and z components no of the normal vector to the surface of the rod, are assumed to be known. To induce bending waves, we apply a displacement at the end (x=0) of the rod, with given components uy0 and uz0. Thus, we have the initial and boundary conditions

  boundsAndInits =
    s[x,y][x,yr,zr,t] no[y][yr,zr]+
      s[x,z][x,yr,zr,t] no[z][yr,zr] == 0,
    s[y,y][x,yr,zr,t] no[y][yr,zr]+
      s[y,z][x,yr,zr,t] no[z][yr,zr] == 0,
    s[y,z][x,yr,zr,t] no[y][yr,zr]+
    s[z,z][x,yr,zr,t] no[z][yr,zr] == 0,
    u[y][0,y,z,t] == uy0[z,t],
    u[z][0,y,z,t] == uz0[y,t],
    s[x,x][0,y,z,t] == 0,
    u[x][x,y,z,0] == 0,
    (D[u[x][x,y,z,t],t] /. t -> 0) == 0,
      u[y][x,y,z,0] == uyt0[x,z],
    u[z][x,y,z,0] == uzt0[x,y],
      (D[u[y][x,y,z,t],t] /. t -> 0) == uyt0p[x,z],
    (D[u[z][x,y,z,t],t] /. t -> 0) == uzt0p[x,y]

    AdditionalFunctions -> {no[y][y,z], no[z][y,z]}]

Asymptotic Expansion

The expansion of all functions as power series in eps produces for the scaled equations, initials, and boundary conditions produces a very long list of PDEs.

  expansion =
    PolyOrderListPDE[Join[scaledEq, boundsAndInits],
      fcts, eps, 4];

  IndexForm[expansion, TabForm -> False,
    AdditionalFunctions -> {no[y][y,z], no[z][y,z]}]

Finding the Lowest Orders

We can determine the lowest orders of the functions involved.

  (ords = FindOrders[Join[scaledEq, boundsAndInits],
    fcts, {uy0[z,t],uz0[y,t]},
    {x, y, z}, eps]) // IndexForm

The order of the parameter b is also fixed.

  FindParameterOrder[scaledEq, ords, b, eps]

Let us substitute the result into the equations.

  newEquations = Join[scaledEq, boundsAndInits] /. %;

Based on this, we have to generate a new expansion.

  newExpansion =
      fcts, eps, 4];

We now drop the vanishing functions from the list of equations.

  IndexForm[newExpansion /. OrderRules[ords, fcts],
    TabForm -> False,
    AdditionalFunctions -> {no[y][y,z], no[z][y,z]}]

The PDEs now have to be solved. This is out of reach for Mathematica (at the moment). For such general boundary conditions, the situation will probably not change very soon. This leaves some work for the skilled applied mathematician.

  On[General::spell, General::spell1];

Up to PerturbationPDE


10-Feb-2016 | Stephan Kaufmann | ZfM | ETH