

Example: Bending Waves in a RodSystemLet us consider bending waves in an isotropic rod (with constant crosssection).
We use the following notation: The xaxis is chosen in direction of the rod, the y,zaxes perpendicular.
FunctionsWe 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];
EquationsThe equations of motion are:
em = { D[s[x,x][x,y,z,t],x]+D[s[x,y][x,y,z,t],y]+ D[s[x,z][x,y,z,t],z] == r D[u[x][x,y,z,t],{t,2}], D[s[x,y][x,y,z,t],x]+D[s[y,y][x,y,z,t],y]+ D[s[y,z][x,y,z,t],z] == r D[u[y][x,y,z,t],{t,2}], D[s[x,z][x,y,z,t],x]+D[s[y,z][x,y,z,t],y]+ 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]  2n(s[x,x][x,y,z,t]+s[y,y][x,y,z,t]+s[z,z][x,y,z,t]), 2 e D[u[z][x,y,z,t],z] == 2(1+n)s[z,z][x,y,z,t]  2n(s[x,x][x,y,z,t]+s[y,y][x,y,z,t]+s[z,z][x,y,z,t]) }; we get the following system of PDEs:
physicalEquations = Join[em, krAndEl];
IndexForm[physicalEquations, Functions > fcts]
ScalingLet us rename the functions and the variables such that we can later reuse the old names for the dimensionless 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 xaxis (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 ParametersWe 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 ConditionsThe 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] };
IndexForm[%, AdditionalFunctions > {no[y][y,z], no[z][y,z]}]
Asymptotic ExpansionThe 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 OrdersWe 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 = PolyOrderListPDE[newEquations, 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 
[Top] 