Chapter 8
Continuous Systems with Quadratic and Cubic Nonlinearities

In this chapter, we use the method of multiple scales to determine second-order uniform asymptotic solutions of continuous or distributed-parameter systems with quadratic and cubic nonlinearities. We consider some of the internal resonances arising from the cubic nonlinearities, namely, one-to-one and three-to-one internal resonances. Two-to-one internal resonances are treated in Chapter 9. We directly attack these continuous systems because treatment of reduced-order models obtained with the Galerkin procedure might lead to erroneous results. We consider buckled beams in Section 8.1, circular cylindrical shells in Section 8.2, and near-square plates in Section 8.3.

Preliminaries

Off[General::spell1, Integrate::gener]

Needs["Utilities`Notation`"]

To determine second-order uniform asymptotic expansions of the solutions of  buckled beams, circular cylindrical shells, and near-square plates by using the method of multiple scales, we introduce the three time scales (Subscript[T, 0]) = t, (Subscript[T, 1]) = ϵ t, and (Subscript[T, 2]) = (ϵ^2) t, where ϵ is a small nondimensional parameter and the order of the amplitude of oscillations. Moreover, we symbolize these scales according to

Symbolize[(Subscript[T, 0])];Symbolize[(Subscript[T, 1])];Symbolize[(Subscript[T, 2])];

timeScales = {(Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])};

In terms of the time scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), the time derivatives can be expressed as

dt[0][expr_] := expr; 
    dt[1][expr_] := 
      Sum[(ϵ^i) D[expr, timeScales[[i + 1]]], {i, 0, 2}]; 
    dt[2][expr_] := 
      (dt[1][dt[1][expr]] // Expand) /. 
        (ϵ^(i_ /; i > 3)) -> 0;

In the course of the analysis, we need the complex conjugates of A and other variables. We define them by using the following rule:

\(\(conjugateRule = {A -> A\&_, A\&_ -> A, 
        Complex[m_, n_] -> Complex[m, \(-n\)]}; \)\)

To speed up manipulations of the integrals involved in the governing equations with Mathematica, we introduce the rules

\(intRule1 = {int[fun_, arg2__] :> int[Expand[fun], arg2]}; 

    intRule2 = {int[a_ + b_, arg2__] :> int[a, arg2] + int[b, arg2], 
        int[ϵ\^n_. \ fun_, arg2__] :> ϵ\^n\ int[fun, arg2], 
        int[a_\ fun_, a1_, b1___] :> 
          a\ int[fun, a1, b1] /; FreeQ[a, First[a1]], 
        int[int[a1_, a2_]\ fun_. , a3__] :> 
          int[a1, a2]\ int[fun, a3] /; FreeQ[Rest[a2], First[a2]]}; \)

To display the outputs in easily read expressions, we introduce the rules

\(\(displayRule = {
        \(\(Derivative[a_, b__]\)[w_\_i_]\)[_, T\_0, __] :> 
          SequenceForm[
            If[arg1 = 
                Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {b}]; 
              arg1 =!= 1, arg1, ""], 
            w\_i\^SequenceForm\ @@\ Table["'", {a}]], 
        \(\(Derivative[a__]\)[A\_i__]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\_i], \(\(Derivative[a__]\)[A\&_\_i__]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\&_\_i], w_\_i_[_, T\_0, __] -> w\_i, A\_i__[__] -> A\_i, 
        A\&_\_i__[__] -> A\&_\_i, 
        Exp[a_.  + b_. \ Complex[0, m_]\ T\_0 + 
              c_. \ Complex[0, n_]\ T\_0] -> 
          Exp[a + \((m*b + n*c)\) I\ T\_0]}; \)\)

8.1 Buckled Beams

In this section, we consider the nonlinear transverse response of a buckled beam possessing either a one-to-one or a three-to-one internal resonance to a principal parametric resonance of the higher mode. The analysis assumes a unimodal  buckled deflection.

The equation governing nonlinear transverse vibrations of a parametrically excited homogeneous beam is (Nayfeh and Mook, 1979)

eq81a = 
      m (Subscript[w, t, t]) + c (Subscript[w, t]) + (Subscript[E, s]) (Subscript[I, s]) (Subscript[w, x, x, x, x]) + 
            (P - (Subscript[f, p]) Cos[Ω t] - 
                  (((Subscript[E, s]) A)/(2 ℓ)) 
                    (∫_0^ℓ)((((Subscript[w, x])))^2) 
                        \[DifferentialD]x) (Subscript[w, x, x]) == 0 /. 
        (Subscript[w, m])__ :> (∂_m) w[x, t];

where P is a static axial load, (Subscript[f, p]) Cos[Ω t] is an axial harmonic excitation, w[x, t] is the transverse deflection of the beam, x is the distance along the undeflected beam, and t is time. Other dimensional parameters are the damping coefficient c and the beam length ℓ, mass per unit length m, cross-sectional area A, modulus of elasticity (Subscript[E, s]), and cross-sectional area moment of inertia (Subscript[I, s]). We note that the subscript s is used with I and E so that Mathematica will not confuse them with (Sqrt[-1])and the exponential function . Equation eq81a is supplemented with the four linear homogeneous boundary conditions

bc81a = {(Subscript[ℬ, 1])[w[0, t]] == 0, 
        (Subscript[ℬ, 2])[w[0, t]] == 0, 
        (Subscript[ℬ, 3])[w[ℓ, t]] == 0, 
        (Subscript[ℬ, 4])[w[ℓ, t]] == 0};

where the (Subscript[ℬ, i]) are linear homogeneous operators.

It is good practice to cast the governing equations and boundary conditions in nondimensional form. Thus, we introduce the following nondimensional variables:

T = (Sqrt[((m (ℓ^4))/( (Subscript[E, s]) (Subscript[I, s])))]);

\(\(nondimRule = {w -> \((ℓ\ w\&~[#1/ℓ, #2/T]\ &)\), 
        x -> ℓ\ x\&~, t -> T\ t\&~, 
        Ω -> Ω\&~/T, 
        P -> P\&~\ E\_s\ I\_s/ℓ\^2, 
        f\_p -> f\&~\ E\_s\ I\_s/ℓ\^2, 
        c -> c\&~\ \@\(E\_s\ I\_s\ m\)/ℓ\^2, 
        A -> 2\ α\ I\_s/ℓ\^2}; \)\)

where the tilde indicates a nondimensional variable. Substituting nondimRule into eq81a and dropping the tildes, we obtain

\(eq81b = 
      \((\(\(\(\(\(ℓ\^3\/\(E\_s\ I\_s\)\) eq81a[\([1]\)] /. 
                      Integrate -> int\) /. nondimRule\) /. s_\&~ -> s\) /. 
                int[a_, {x\ ℓ, __}] -> 
                  ℓ\ int[a, {x, 0, 1}] // Expand\) // PowerExpand)
          \) == 0; 
eq81b /. int -> Integrate\)

8.1.1 Postbuckling Deflection

The first step in analyzing the nonlinear vibrations of a buckled beam is the determination of the critical Euler buckling loads (Subscript[P, k]) and their corresponding mode shapes (Subscript[ϕ, k])[x]. We symbolize the (Subscript[P, k]) as

Symbolize[(Subscript[P, k])];

Dropping the time derivative and neglecting the forcing and nonlinear terms in eq81b, we define the linear buckling problem as

buckEq = 
      eq81b /. int -> (0 &) /. {P -> (Subscript[P, k]), f -> 0, 
          w -> ((Subscript[ϕ, k])[#1] &)}

Solving buckEq for (Subscript[ϕ, k]^((4)))[x] yields

buckRule=Solve[buckEq,(Subscript[ϕ, k]^((4)))[x]][[1]]

Next, we increase the axial load P beyond (Subscript[P, k]) and assume that the postbuckling displacement is b (Subscript[ϕ, k]), where b is a nondimensional measure of the buckling level. Substituting this displacement into eq81b, dropping the forcing term and the time derivative, and using intRule2 and buckRule, we obtain

bEq = 
      Factor[eq81b[[1]] /. {f -> 0, w -> (b (Subscript[ϕ, k])[#1] &)}
                   /. intRule2 /. buckRule] == 0

It follows from bEq that either b = 0 (unbuckled case) or b is given by

(bRule = (b^2) -> ((b^2) /. Solve[bEq, b][[2]])) /. 
      int -> Integrate

And solving bEq for P yields

(PRule = Solve[bEq, P][[1]]) /. int -> Integrate

Because P is greater than (Subscript[P, k]), PRule yields the amplitude of the kth buckling mode.

Having solved the postbuckling problem, we assume that the beam deflection is the sum of the static buckled displacement b (Subscript[ϕ, k])[x], corresponding to the kth buckled mode, and a time-dependent relative deflection u[x, t]; that is,

wRule = {w -> (b (Subscript[ϕ, k])[#1] + u[#1, #2] &)};

Substituting wRule into eq81b and using intRule1, intRule2, buckRule, and PRule, we find that the nonlinear transverse vibrations of the beam around its buckled configuration are given by

EOM = 
      eq81b /. wRule /. intRule1 //. intRule2 /. buckRule /. 
          PRule // ExpandAll; EOM /. int -> Integrate

BC = bc81a /. {w -> u, ℓ -> 1}

In the following analysis, we focus on hinged-hinged, first mode-buckled beams. Then, the  boundary conditions for the static deflection (Subscript[ϕ, k])[x] can be written as

bBC={(Subscript[ϕ, k])[0]==0,(Subscript[ϕ, k]^″)[0]==0,(Subscript[ϕ, k])[1]==0,(Subscript[ϕ, k]^″)[1]==0};

and the boundary conditions for the dynamic deflection u[x, t] can be written as

BC1 = {u[0, t] == 0, (∂_(x, x))u[x, t] == 0 /. x -> 0, 
        u[1, t] == 0, (∂_(x, x))u[x, t] == 0 /. x -> 1}

For the case of first buckled mode, it follows from buckEq and bBC that the first critical buckling load and associated modeshape are

\(bucklingLoad = P\_k -> π\^2; 

    modeshape = {ϕ\_k -> \((Sin[π\ #]\ &)\)}; \)

where the buckling modeshape was normalized so that (Subscript[ϕ, k])[1/2] = 1.

Linear Vibrations

The next step in analyzing the nonlinear vibrations of a buckled beam is the calculation of its linear natural frequencies and corresponding modeshapes. To this end, we linearize EOM in u[x, t], drop the forcing, damping, and nonlinear terms, and obtain

\(linearEq = 
      \(\(\(\(EOM /. {c -> 0, f -> 0}\) /. 
              u -> \((ϵ\ u[#1, #2]\ &)\)\) /. intRule2\) /. 
          ϵ\^\(n_ /; n > 1\) -> 0\) /. ϵ -> 1; 

    linearEq /. int -> Integrate\)

The boundary conditions are the same as BC1, which we list below as

linearBC = BC1

Then, using separation of variables, we seek the solution of linearEq and linearBC in the form

uSol = u -> Function[{x, t}, Φ[x] Exp[I ω t]];

where ω is the nondimensional linear natural frequency and Φ[x] is the corresponding modeshape. Substituting uSol into linearEq and linearBC and using intRule2 and Thread, we obtain the equation governing Φ[x] as

\(ΦEq = 
      \(\(linearEq /. uSol\) /. intRule2 // 
          Thread[E\^\(\(-I\)\ ω\ t\)\ #, Equal]\ &\) // ExpandAll; 

    ΦEq /. int -> Integrate\)

and the boundary conditions as

ΦBC = linearBC /. u -> (Φ[#1] &)

To find the general solution of ΦEq and ΦBC, we follow Nayfeh, Kreider, and Anderson (1995) and treat the definite integral as a constant B; thus, we have a nonhomogeneous linear ordinary-differential equation. Hence, its general solution can be expressed as a linear combination of homogeneous and particular solutions. Because (Subscript[ϕ, k])[x] is a solution of buckEq, the particular solution is proportional to ((Subscript[ϕ, k])^″)[x]. Consequently, we express the general solution of ΦEq and ΦBC as

cList = Table[(Subscript[c, i]), {i, 5}];

shapeFunc = {Cos[(Subscript[λ, 1]) x], Sin[(Subscript[λ, 1]) x], 
        Cosh[(Subscript[λ, 2]) x], Sinh[(Subscript[λ, 2]) x], 
        ((Subscript[ϕ, k])^″)[x]};

ΦSol = 
      Φ -> 
        Function[x, cList . shapeFunc /. modeshape // Evaluate]

where the (Subscript[λ, i]) satisfy the characteristic equation

charEq = 
      ΦEq /. b -> 0 /. bucklingLoad /. 
            Φ -> (Exp[λ #] &) // 
          Thread[(E^(-λ x)) #, Equal] & // ExpandAll

λRule = {-(Subscript[λ, 1]^2), (Subscript[λ, 2]^2)} -> 
          ((λ^2) /. Solve[charEq, λ] // Union) // Thread

Substituting bucklingLoad, modeshape, and ΦSol into the left-hand sides of ΦBC and ΦEq and using intRule1 and intRule2, we obtain

expr1 = 
      Append[#[[1]] & /@ ΦBC /. ΦSol, 
        ΦEq[[1]] /. bucklingLoad /. modeshape
                   /. ΦSol /. intRule1 //. intRule2 /. 
          int -> Integrate]

The coefficient matrix of cList in expr1 can be obtained as

coefMat = Outer[Coefficient, expr1, cList]

Substituting charEq into coefMat[[5]] and eliminating the common factor Sin[π x], we have

coefMat[[5]] = Coefficient[coefMat[[5]], Sin[π x]]

Setting the determinant of coefMat equal to zero yields the following characteristic equation for the natural frequencies:

Factor[Det[coefMat]] == 0

There are two possibilities: (b^2) (π^4) α - (ω^2) = 0 and Sin[(Subscript[λ, 1])]=0. Hence, either

ωSol1 = 
      Solve[(b^2) (π^4) α - (ω^2) == 0, ω][
        [2]]

or

ωSol2 = 
      Solve[-(Subscript[λ, 1]^2) + (((n π))^2) == 0 /. λRule, 
          ω][[2]]

for n > 1. The lowest five natural frequencies are

frequencies = {(Subscript[ω, 1]) -> (ω /. ωSol1), 
        (Subscript[ω, 2]) -> (ω /. ωSol2 /. n -> 2), 
        (Subscript[ω, 3]) -> (ω /. ωSol2 /. n -> 3), 
        (Subscript[ω, 4]) -> (ω /. ωSol2 /. n -> 4), 
        (Subscript[ω, 5]) -> (ω /. ωSol2 /. n -> 5)}

The first two frequencies vs. the buckling level are plotted as follows:

\(\(Plot[
      ω /. {ωSol1 /. b\ \@α -> s, 
            ωSol2 /. n -> 2} // Evaluate, {s, 0, 15}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1]}, Frame -> True, 
      FrameLabel -> {\*""b\!\(\@α\)"", "ω"}, 
      RotateLabel -> False, 
      Epilog -> {Text["C", {2 \@ 3, 2 \(\@ 3\) π\^2 + 5}], 
          Text[\*""\!\(b\^\(1 : 3\)\)"", {2/\@3, 3}], 
          Text[\*""\!\(b\^\(1 : 1\)\)"", {2 \@ 3, 3}], 
          Text[\*""\!\(b\^\(3 : 1\)\)"", {6 \@ 3, 3}]}]; \)\)

where the low-frequency mode (red line) is symmetric and the second mode (blue line) is antisymmetric. As the buckling level increases from zero, the frequency (Subscript[ω, 1]) of the first mode increases from zero and crosses the frequency (Subscript[ω, 2]) of the second mode at point C. Hence, there are two possible three-to-one internal resonances: (Subscript[ω, 2]) = 3 (Subscript[ω, 1]) when b (Sqrt[ α]) = (b^(1 : 3)) and (Subscript[ω, 1]) = 3 (Subscript[ω, 2]) when (b^(3 : 1)). There is also a possible one-to-one internal resonance (Subscript[ω, 2]) = (Subscript[ω, 1]) when b (Sqrt[ α]) = (b^(1 : 1)). Whether or not these candidates for internal resonances will be activated depends on the corresponding modeshapes. In this section, we consider the case of one-to-one internal resonance between the first and second modes. In Section 8.1.4, we consider three-to-one internal resonances between the first two modes. In Sections 8.1.5 and 8.1.6, we consider one-to-one and three-to-one internal resonances between the first and third modes and the first and fouth modes, respectively.

In the next section, we consider a combination of a one-to-one internal resonance between the lowest two modes and a principal parametric resonance of the second mode; that is,

\(omgList = {ω\_1, ω\_2}; 

    ResonanceConds = {
        ω\_2 == ω\_1 + ϵ\^2\ σ\_1, 
        Ω == 2  ω\_2 + ϵ\^2\ σ\_2}; \)

where the (Subscript[σ, i]) are detuning parameters that describe the nearness of the resonances. The buckling level for which the natural frequencies of the first and second modes are equal (i.e., (Subscript[ω, 2]) = (Subscript[ω, 1])) is given by

values1 = 
      Solve[(Subscript[ω, 1]) == (Subscript[ω, 2]) /. frequencies, b][[1]]

The corresponding (Subscript[λ, m]) are

λRuleN = 
      Solve[λRule /. ω -> (Subscript[ω, 1]) /. frequencies
                 /. values1 /. Rule -> Equal, {(Subscript[λ, 1]), 
            (Subscript[λ, 2])}][[4]]

In order that the influence of the damping and nonlinearity balance the influence of the forcing, we scale the damping and forcing terms as

scaling = {c -> 2  (ϵ^2) μ , f -> (ϵ^2) f};

Substituting for the natural frequency (Subscript[ω, 1]) back into coefMat, we determine the modeshape as

NullSpace[coefMat /. ω -> (Subscript[ω, 1]) /. frequencies][
          [1]] . shapeFunc /. modeshape

which we normalize as

shape1 = 
      c Sin[π x] /. 
        Solve[(∫_0^1) (c^2) Sin[π x]^2 
                  \[DifferentialD]x == 1, c][[2]]

Substituting (Subscript[ω, 2]) into coefMat, we determine the modeshape as

NullSpace[
            coefMat /. λRuleN /. ω -> (Subscript[ω, 2]) /. 
              frequencies][[1]] . shapeFunc /. λRuleN

which we also normalize as

shape2 = 
      c Sin[2  π x] /. 
        Solve[(∫_0^1) (c^2) Sin[2  π x]^2 
                  \[DifferentialD]x == 1, c][[2]]

In a similar fashion, we find that the next three normalized modeshapes are

shape3 = (Sqrt[2]) Sin[3 π x];

shape4 = (Sqrt[2]) Sin[4 π x];

shape5 = (Sqrt[2]) Sin[5 π x];

To treat a one-to-one internal resonance between the first and second modes, we define (Subscript[ϕ, k]) and the first two eigenmodes of the buckled beam as

modeshapes = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape1 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape2 // Evaluate]}];

8.1.2 Perturbation Analysis

We use the method of multiple scales to directly attack EOM and BC1. To transform the time derivatives in EOM in terms of the scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), we let

multiScales = {u[x_, t] -> u[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
        Derivative[m_, n_][u_][x_, t] -> 
          dt[n][D[u[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], {x, m}]], t -> (Subscript[T, 0])};

Then, we seek a second-order approximate solution in the form

solRule = 
      u -> (Evaluate[Sum[(ϵ^j) (Subscript[u, j])[#1, #2, #3, #4], {j, 3}]] &)

where #1 stands for x and #2, #3, and #4 stand for (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), respectively.

Substituting multiScales, solRule, and scaling into EOM and BC1, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq81c = 
      (Join[{EOM}, BC1] /. multiScales /. solRule /. scaling
                     /. intRule1 //. intRule2 // TrigToExp // ExpandAll)
           /. (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

\(eqEps = 
      \(\(Thread[CoefficientList[Subtract\ @@\ #, ϵ] == 0]\ &\)\ /@\ 
            eq81c // Transpose\) // Rest; 
eqEps /. displayRule\)

First-Order Solution

Because in the presence of damping all of the modes that are not directly or indirectly excited decay with time, the solution of eqEps[[1]] is taken to consist of the two modes that might be involved in the internal resonance; that is,

\(\(sol1 = {
        u\_1 -> Function[{x, T\_0, T\_1, T\_2}, 
            Sum[A\_i[T\_1, T\_2]\ Φ\_i[x]\ 
                    Exp[I\ ω\_i\ T\_0] + 
                  A\&_\_i[T\_1, T\_2]\ Φ\_i[x]\ 
                    Exp[\(-I\)\ ω\_i\ T\_0], {i, 2}] // Evaluate]}; 
    \)\)

where the (Subscript[Φ, i]) are the eigenmodes, which satisfy the orthonormality condition \(TraditionalForm
      \`∫\_0\%1 Φ\_n[x]\ Φ\_m[x]\ 
            dx = δ\_\(n, m\)\).

Second-Order Solution

Substituting sol1 into eqEps[[2,1]] and using intRule1 and intRule2 yields

order2Eq = 
      (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 2])) == 
                (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 2])) - 
                  eqEps[[2, 1, 1]] /. sol1 /. intRule1 //. 
          intRule2 // ExpandAll; order2Eq /. displayRule

To collect the terms that might lead to secular terms from the right-hand side of order2Eq, we define the rules:

OmgRule = 
      Solve[ResonanceConds, 
              Complement[omgList, {#}]~Join~{Ω}][[1]] &
           /@ omgList

expRule1[i_] := 
      Exp[arg_] :> 
        Exp[Expand[arg /. OmgRule[[i]]] /. (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])]

Collecting the terms that may lead to secular terms, the terms proportional to (E^(I (Subscript[ω, i]) (Subscript[T, 0]))), we have

\(ST1 = 
      \(Coefficient[order2Eq[\([2]\)] /. expRule1[#], 
            E\^\(I\ ω\_#\ T\_0\)]\ &\)\ /@\ {1, 2}; 

    ST1 /. displayRule\)

Because the homogeneous parts of order2Eq and corresponding boundary conditions have a nontrivial solution, the corresponding nonhomogeneous problem has a solution only if solvability conditions are satisfied. These conditions demand that ST1 be orthogonal to every solution of the adjoint homogeneous problem. In this case, the problem is self-adjoint and hence solutions of the adjoint problem are given by the (Subscript[Φ, i])[x]. Demanding that ST1[[i]] be orthogonal to the (Subscript[Φ, i])[x], we obtain the solvability conditions

SCond1 = 
      Table[(∫_0^1) ST1[[k]] (Subscript[Φ, k])[x] 
                \[DifferentialD]x == 0, {k, 2}] /. modeshapes

SCond1Rule=Solve[SCond1,{(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]

whose complex conjugate is

ccSCond1Rule = SCond1Rule /. conjugateRule

Substituting SCond1Rule and ccSCond1Rule into order2Eq and using modeshapes, we have

\(order2Eqm = 
      \(\(\(\(order2Eq /. SCond1Rule\) /. ccSCond1Rule\) /. modeshapes\) /. 
          int -> Integrate\) /. Integrate -> int; 

    order2Eqm /. displayRule\)

The associated boundary conditions are

order2BC = eqEps[[2]] // Rest

We use the method of undetermined coefficients to determine a particular solution of order2Eqm and order2BC. To accomplish this, we first determine the forms of the nonhomogeneous terms (i.e., the terms on the right-hand side of order2Eqm) as

rhsTerms = 
      Cases[order2Eqm[[2]], #] & /@
             {a_ f_[b_ x] :> 
                f[b x] /; FreeQ[a, (Subscript[T, 0])], _ (E^(a_ (Subscript[T, 0]) + b_. )) 
                  f_[c_ x] -> (E^(a (Subscript[T, 0]) + b)) f[c x]} // Flatten // 
        Union

which all satisfy the order2BC. Since we have only even spatial derivatives on the left-hand side of order2Eqm, we seek a particular solution as a linear combination of the rhsTerms as follows:

symbolList = Table[(Subscript[α, i]), {i, Length[rhsTerms]}]

sol2Form = 
      (Subscript[u, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
          symbolList . rhsTerms // Evaluate];

Substituting sol2Form into order2Eqm, equating coefficients of like terms, and solving for the symbolList, we have

symbolRule = 
      Solve[Coefficient[
                  Subtract @@ order2Eqm /. sol2Form /. intRule1
                         //. intRule2 /. int -> Integrate, rhsTerms] == 
                0 /. Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Thread, symbolList][[1]]

Substituting for the parameter values in symbolRule yields

symbolRuleN = 
      symbolRule /. Ω -> 2  (Subscript[ω, 2]) /. frequencies
             /. bucklingLoad /. values1

Substituting symbolRuleN into sol2Form yields

sol2 = sol2Form /. symbolRuleN

Solvability Conditions

Substituting sol1 and sol2 into eqEps[[3,1]] and using intRule1 and intRule2 yields

order3Eq = 
        (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 3])) == 
          ((eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 3])) - 
                            eqEps[[3, 1, 1]] /. sol1 /. sol2 /. 
                      intRule1 //. intRule2 /. modeshapes /. 
                int -> Integrate // Expand);  // Timing

Collecting the terms that may lead to secular terms, the terms proportional to (E^(I (Subscript[ω, i]) (Subscript[T, 0]))), we have

\(ST2 = 
      \(Coefficient[order3Eq[\([2]\)] /. expRule1[#], 
            E\^\(I\ ω\_#\ T\_0\)]\ &\)\ /@\ {1, 2}; 

    ST2 /. displayRule\)

Demanding that ST2[[i]] be orthogonal to the (Subscript[Φ, i])[x], we obtain the solvability conditions

\(SCond2 = 
      \(\(\(\(\(Table[
                    int[ST2[\([k]\)]\ Φ\_k[x], {x, 0, 1}] == 0, {
                      k, 2}] /. intRule1\) //. intRule2\) /. modeshapes\) /. 
            values1\) /. int -> Integrate\) /. D[SCond1Rule, T\_1]; 

    SCond2 /. displayRule\)

8.1.3 The Function MMSDirect11

Collecting the commands in the preceding section, we build a function MMSDirect11, specifically for the nonlinear transverse vibrations of a hinged-hinged beam around its first buckled mode (EOM and BC1), to determine the solvability conditions for different combinations of frequencies, modeshapes, and ResonanceConds in which the first mode might be involved in either a one-to-one or a three-to-one internal resonance with one of the other modes.

MMSDirect11[frequencies_, modeshapes_, ResonanceConds_]:=
	Module[{},
		multiScales = {u[x_, t] -> u[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
                Derivative[m_, n_][u_][x_, t] -> 
                  dt[n][D[u[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], {x, m}]], t -> (Subscript[T, 0])};
     solRule = 
              u -> (Evaluate[
                      Sum[(ϵ^j) (Subscript[u, j])[#1, #2, #3, #4], {j, 3}]] &)
                  ;
		scaling = {c -> 2  (ϵ^2) μ , 
                f -> (ϵ^2) f};
     eqs = (
                  Join[{EOM}, BC1] /. multiScales /. solRule /. 
                            scaling /. intRule1 //. intRule2 // TrigToExp
                       // ExpandAll) /. 
                (ϵ^(n_ /; n > 3)) -> 0;
     eqEps = 
              Thread[CoefficientList[Subtract @@ #, ϵ] == 0] &
                       /@ eqs // Transpose // Rest;
     \(sol1 = {
                u\_1 -> 
                  Function[{x, T\_0, T\_1, T\_2}, 
                    Sum[A\_i[T\_1, T\_2]\ Φ\_i[x]\ 
                            Exp[I\ ω\_i\ T\_0] + 
                          A\&_\_i[T\_1, T\_2]\ Φ\_i[x]\ 
                            Exp[\(-I\)\ ω\_i\ T\_0], {i, 2}] // 
                      Evaluate]}\);
     order2Eq = 
              (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 2])) == 
                        (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 2])) - 
                          eqEps[[2, 1, 1]] /. sol1 /. intRule1 //. 
                  intRule2 // ExpandAll;
		omgList = {(Subscript[ω, 1]), (Subscript[ω, 2])};
     OmgRule = 
              Solve[ResonanceConds, 
                      Complement[omgList, {#}]~Join~{Ω}][
                    [1]] & /@ omgList;
     expRule1[i_] := 
              Exp[arg_] :> 
                Exp[Expand[arg /. OmgRule[[i]]] /. 
                    (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])];
     \(ST1 = \(Coefficient[order2Eq[\([2]\)] /. expRule1[#], 
                    E\^\(I\ ω\_#\ T\_0\)]\ &\)\ /@\ {1, 2}\);
     SCond1 = 
              Table[(∫_0^1) ST1[[k]] (Subscript[Φ, k])[x] 
                        \[DifferentialD]x == 0, {k, 2}] /. modeshapes;
     SCond1Rule=Solve[SCond1,{(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]];
     ccSCond1Rule = SCond1Rule /. conjugateRule;
     order2Eqm = 
              order2Eq /. SCond1Rule /. ccSCond1Rule /. modeshapes
                     /. int -> Integrate /. Integrate -> int;
     rhsTerms = 
              Cases[order2Eqm[[2]], #] & /@
                     {a_ f_[b_ x] :> 
                        f[b x] /; FreeQ[a, (Subscript[T, 0])], _ (E^(a_ (Subscript[T, 0]) + b_. )) 
                          f_[c_ x] -> (E^(a (Subscript[T, 0]) + b)) f[c x]} // Flatten
                   // Union;
     symbolList = Table[(Subscript[α, i]), {i, Length[rhsTerms]}];
     sol2Form = 
              (Subscript[u, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
                  symbolList . rhsTerms // Evaluate];
     symbolRule = 
              Solve[Coefficient[
                          Subtract @@ order2Eqm /. sol2Form /. 
                                intRule1 //. intRule2 /. 
                            int -> Integrate, rhsTerms] == 0 /. 
                      Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Thread, symbolList][
                [1]];
 	 values1 = 
              Solve[Select[ResonanceConds, FreeQ[#, Ω] &] /. 
                      ϵ -> 0 /. frequencies, b] // Flatten;
     symbolRuleN = 
              symbolRule /. Ω -> 2  (Subscript[ω, 2]) /. 
                    frequencies /. bucklingLoad /. values1;
     sol2 = sol2Form /. symbolRuleN;
     order3Eq = 
              (eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 3])) == 
                ((eqEps[[1, 1, 1]] /. (Subscript[u, 1]) -> (Subscript[u, 3])) - 
                                  eqEps[[3, 1, 1]] /. sol1 /. sol2 /. 
                            intRule1 //. intRule2 /. modeshapes /. 
                      int -> Integrate // Expand);
     \(ST2 = \(Coefficient[order3Eq[\([2]\)] /. expRule1[#], 
                    E\^\(I\ ω\_#\ T\_0\)]\ &\)\ /@\ {1, 2}\);
     SCond2 = 
              Table[
                            int[ST2[[k]] (Subscript[Φ, k])[x], {x, 0, 
                                  1}] == 0, {k, 2}] /. intRule1 //. intRule2
                         /. modeshapes /. values1 /. 
                  int -> Integrate /. D[SCond1Rule, (Subscript[T, 1])];
 	 Return[
              SCond2 /. (Subscript[A, i])_ -> ((Subscript[A, i])[#2] &) /. 
                f_[(Subscript[T, 1]), (Subscript[T, 2])] -> f[(Subscript[T, 2])]]
 ]

8.1.4 Three-to-One Internal Resonances Between the First Two Modes

In this section, we consider the case (Subscript[ω, 2]) ≈ 3  (Subscript[ω, 1]) and Ω ≈ 2  (Subscript[ω, 2]). To describe the nearness of the resonances, we introduce the detuning parameters (Subscript[σ, i]) defined by

ResonanceConds1 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies, modeshapes, ResonanceConds1] // Timing

which indicates that this three-to-one internal resonance is not activated.

Next, we consider the case in which the frequency of the first mode is three times that of the second mode; that is, the buckling load is above the first crossover value. To accomplish this, we redefine the frequencies and denote the frequency of the first mode by (Subscript[ω, 2]) and that of the second mode by (Subscript[ω, 1]) as follows:

frequencies2 = {(Subscript[ω, 2]) -> (ω /. ωSol1), 
        (Subscript[ω, 1]) -> (ω /. ωSol2 /. n -> 2)}

Moreover, we redefine the modeshapes according to

modeshapes2 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape2 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape1 // Evaluate]}];

To describe the nearness of the resonances, we introduce the two detuning parameters (Subscript[σ, i]) defined by

ResonanceConds2 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies2, modeshapes2, ResonanceConds2] // Timing

Again, this three-to-one internal resonance is not activated.

8.1.5 One-to-One and Three-to-One Internal Resonances Between the First and Third Modes

Variations of the first and third natural frequencies with the buckling level are plotted as follows:

It follows from this figure that the frequencies of the first and third modes are equal when b (Sqrt[ α]) = (b^(1 : 1)), the frequency of the third mode is three times that of the first mode when b (Sqrt[ α]) = (b^(1 : 3)),  and the frequency of the first mode is three times that of the third mode when b (Sqrt[ α]) =   (b^(3 : 1)).  Next, we use the function MMSDirect11 to treat these three cases.

One-to-One Internal Resonance

We denote the frequency of the third mode by (Subscript[ω, 2]) and that of the first mode by (Subscript[ω, 1]) as follows:

frequencies3 = {(Subscript[ω, 1]) -> (ω /. ωSol1), 
        (Subscript[ω, 2]) -> (ω /. ωSol2 /. n -> 3)}

Moreover, we define (Subscript[ϕ, k]) and the first and third eigenmodes of the buckled beam as

modeshapes3 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape1 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape3 // Evaluate]}];

We consider a combination of a one-to-one internal resonance between the first and third modes and a principal parametric resonance of the third mode; that is,

ResonanceConds3 = {
        (Subscript[ω, 2]) == (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies3, modeshapes3, ResonanceConds3] // Timing

which indicates that the one-to-one internal resonance between the first and third modes is activated.

Three-to-One Internal Resonances

To treat the case in which the natural frequency of the third mode is approximately three times that of the first mode using the function MMSDirect11, we denote these frequencies by (Subscript[ω, 2]) and (Subscript[ω, 1]), respectively, according to

frequencies4 = {(Subscript[ω, 1]) -> (ω /. ωSol1), 
        (Subscript[ω, 2]) -> (ω /. ωSol2 /. n -> 3)}

Moreover, we define (Subscript[ϕ, k]) and the first and third eigenmodes of the buckled beam as

modeshapes4 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape1 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape3 // Evaluate]}];

We consider a combination of a three-to-one internal resonance between the first and third modes and a principal parametric resonance of the third mode; that is,

ResonanceConds4 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies4, modeshapes4, ResonanceConds4] // Timing

which indicates that this three-to-one internal resonance between the first and third modes is not activated.

When the frequency of first mode is approximately three times that of the third mode, we define the frequencies as

frequencies5 = {(Subscript[ω, 2]) -> (ω /. ωSol1), 
        (Subscript[ω, 1]) -> (ω /. ωSol2 /. n -> 3)}

Moreover, we define (Subscript[ϕ, k]) and the third and first eigenmodes of the buckled beam as

modeshapes5 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape3 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape1 // Evaluate]}];

We consider a combination of a three-to-one internal resonance between the first and third modes and a principal parametric resonance of the third mode; that is,

ResonanceConds5 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies5, modeshapes5, ResonanceConds5] // Timing

which indicates that this three-to-one internal resonance between the first and third modes is not activated.

8.1.6 One-to-One and Three-to-One Internal Resonances Between the First and Fourth Modes

Variations of the first and fourth natural frequencies with the buckling level are plotted as follows:

Again, there exist buckling levels for which the frequencies of the first and fourth modes are in the ratio of 1:1, 3:1, and 1:3. We use the function MMSDirect11 to treat these cases.

One-to-One Internal Resonance

When the frequencies of the first and fourth modes are approximately equal, we denote them by (Subscript[ω, 1]) and (Subscript[ω, 2]), respectively, according to

frequencies6 = {(Subscript[ω, 1]) -> (ω /. ωSol1), 
        (Subscript[ω, 2]) -> (ω /. ωSol2 /. n -> 4)}

Moreover, we define (Subscript[ϕ, k]) and the first and fourth eigenmodes of the buckled beam as

modeshapes6 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape1 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape4 // Evaluate]}];

We consider a combination of a one-to-one internal resonance between the first and fourth modes and a principal parametric resonance of the fourth mode; that is,

ResonanceConds6 = {
        (Subscript[ω, 2]) == (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies6, modeshapes6, ResonanceConds6] // Timing

which indicates that this one-to-one internal resonance between the first and fourth modes is activated.

Three-to-One Internal Resonances

When the frequency of the fourth mode is approximately three times that of the first mode, we denote these frequencies by (Subscript[ω, 2]) and (Subscript[ω, 1]), respectively, as

frequencies7 = {(Subscript[ω, 1]) -> (ω /. ωSol1), 
        (Subscript[ω, 2]) -> (ω /. ωSol2 /. n -> 4)}

Moreover, we define (Subscript[ϕ, k]) and the first and fourth eigenmodes of the buckled beam as

modeshapes7 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape1 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape4 // Evaluate]}];

We consider a combination of a three-to-one internal resonance between the first and fourth modes and a principal parametric resonance of the fourth mode; that is,

ResonanceConds7 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies7, modeshapes7, ResonanceConds7] // Timing

which indicates that this three-to-one internal resonance between the first and fourth modes is not activated.

When the frequency of the first mode is approximately three times that of the fourth mode, we denote these frequencies by (Subscript[ω, 2]) and (Subscript[ω, 1]), respectively, as

frequencies8 = {(Subscript[ω, 2]) -> (ω /. ωSol1), 
        (Subscript[ω, 1]) -> (ω /. ωSol2 /. n -> 4)}

Moreover, we define (Subscript[ϕ, k]) and the fourth and first eigenmodes of the buckled beam as

modeshapes8 = 
      Join[modeshape, {(Subscript[Φ, 1]) -> Function[x, shape4 // Evaluate], 
          (Subscript[Φ, 2]) -> Function[x, shape1 // Evaluate]}];

We consider a combination of a three-to-one internal resonance between the first and fourth modes and a principal parametric resonance of the fourth mode; that is,

ResonanceConds8 = {
        (Subscript[ω, 2]) == 3  (Subscript[ω, 1]) + (ϵ^2) (Subscript[σ, 1]), 
        Ω == 2  (Subscript[ω, 2]) + (ϵ^2) (Subscript[σ, 2])};

Using MMSDirect11, we obtain the solvability conditions

MMSDirect11[frequencies8, modeshapes8, ResonanceConds8] // Timing

which indicates again that this three-to-one internal resonance between the first and fourth modes is not activated.

8.2 Circular Cylindrical Shells

Following McIvor (1966) and Raouf and Nayfeh (1990), we write the equations of motion governing the dynamic response of an infinitely long cylindrical shell to a primary-resonance excitation of one of its two orthogonal flexural modes as

EOM = {
          w - (Subscript[ψ, θ]) + (Subscript[w, t, t]) + 
              (α^2) 
                (w + 2 (Subscript[w, θ, θ]) + 
                    (Subscript[w, θ, θ, θ, θ])) == 
            -((Subscript[w, θ]^2)/2) - w (Subscript[w, θ]^2) - (Subscript[ψ, t]^2) + 
              w (Subscript[ψ, t]^2) - 2 w (Subscript[ψ, θ]) + 
              (1/2) (Subscript[w, θ]^2) (Subscript[ψ, θ]) + (Subscript[ψ, θ]^2) - 
              w (Subscript[ψ, θ]^2) + 
              f (ϵ^3) Cos[n θ] Cos[Ω t] 
                (1 - w + (Subscript[ψ, θ])) - 
              w (Subscript[w, θ, θ]) - 
              (w^2) (Subscript[w, θ, θ]) + 
              (3/2) (Subscript[w, θ]^2) (Subscript[w, θ, θ]) + 
              (Subscript[ψ, θ]) (Subscript[w, θ, θ]) + 
              w (Subscript[ψ, θ]) (Subscript[w, θ, θ]) - 
              (Subscript[ψ, θ]^2) (Subscript[w, θ, θ]) + 
              (Subscript[w, θ]) (Subscript[ψ, θ, θ]) + 
              w (Subscript[w, θ]) (Subscript[ψ, θ, θ]) - 
              2 (Subscript[w, θ]) (Subscript[ψ, θ]) 
                (Subscript[ψ, θ, θ]) + 
              (α^2) 
                (-3 (w^2) - ((11 (Subscript[w, θ]^2))/2) - 
                    11 w (Subscript[w, θ, θ]) + 
                    4 (Subscript[ψ, θ]) (Subscript[w, θ, θ]) - 
                    6 (Subscript[w, θ, θ]^2) + 
                    4 (Subscript[w, θ]) (Subscript[ψ, θ, θ]) - 
                    8 (Subscript[w, θ]) (Subscript[w, θ, θ, θ]) + 
                    8 (Subscript[ψ, θ, θ]) 
                      (Subscript[w, θ, θ, θ]) + 
                    w (Subscript[ψ, θ, θ, θ]) + 
                    5 (Subscript[w, θ, θ]) 
                      (Subscript[ψ, θ, θ, θ]) - 
                    4 w (Subscript[w, θ, θ, θ, θ]) + 
                    4 (Subscript[ψ, θ]) 
                      (Subscript[w, θ, θ, θ, θ]) + 
                    (Subscript[w, θ]) 
                      (Subscript[ψ, θ, θ, θ, θ])), 
          (Subscript[w, θ]) + (Subscript[ψ, t, t]) - (Subscript[ψ, θ, θ]) == 
            2 w (Subscript[w, θ]) + ((Subscript[w, θ]^3)/2) + 2 (Subscript[w, t]) (Subscript[ψ, t]) - 
              2 w (Subscript[w, t]) (Subscript[ψ, t]) - 2 (Subscript[w, θ]) (Subscript[ψ, θ]) + 
              2 w (Subscript[w, θ]) (Subscript[ψ, θ]) + 
              (Subscript[w, θ]) (Subscript[w, θ, θ]) + 
              w (Subscript[w, θ]) (Subscript[w, θ, θ]) - 
              2 (Subscript[w, θ]) (Subscript[ψ, θ]) (Subscript[w, θ, θ]) + 
              2 w (Subscript[ψ, t, t]) - (w^2) (Subscript[ψ, t, t]) - 
              2 w (Subscript[ψ, θ, θ]) + 
              (w^2) (Subscript[ψ, θ, θ]) - 
              (Subscript[w, θ]^2) (Subscript[ψ, θ, θ]) + 
              (α^2) 
                ((Subscript[w, θ]) (Subscript[w, θ, θ]) - 
                    w (Subscript[w, θ, θ, θ]) - 
                    (Subscript[w, θ, θ]) 
                      (Subscript[w, θ, θ, θ]) + 
                    (Subscript[w, θ]) 
                      (Subscript[w, θ, θ, θ, θ]))} /. {
          w -> w[θ, t], 
          (Subscript[((s : w | ψ)), m])__ :> (∂_m) s[θ, t]};

We seek a second-order uniform asymptotic expansion of the solution of EOM in the form

solRule = {
        w -> (Evaluate[Sum[(ϵ^j) (Subscript[w, j])[#1, #2, #3, #4], {j, 3}]] &)
            , ψ -> 
          (Evaluate[Sum[(ϵ^j) (Subscript[ψ, j])[#1, #2, #3, #4], {j, 3}]]
                 &)}

where #1 stands for θ and #2, #3, and #4 stand for (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), respectively. Transforming the time derivatives in EOM in terms of the scales (Subscript[T, 0]), (Subscript[T, 1]), and (Subscript[T, 2]), substituting solRule into EOM, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq82a = 
      (EOM /. {w[x_, t] -> w[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
                  ψ[x_, t] -> ψ[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
                  Derivative[m_, n_][w_][x_, t] -> 
                    dt[n][D[w[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], {x, m}]], t -> (Subscript[T, 0])}
                   /. solRule // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we have

eqEps = 
      Rest[Thread[CoefficientList[Subtract @@ #, ϵ] == 0]] &
             /@ eq82a // Transpose; eqEps /. displayRule

8.2.1 First-Order Solution

To determine the solution of the first-order problem, eqEps[[1]], we use the method of separation of variables and let

svRule = {
        (Subscript[w, 1]) -> (A (E^(I m #1)) (E^(I (Subscript[ω, m]) #2)) &), 
        (Subscript[ψ, 1]) -> (B (E^(I m #1)) (E^(I (Subscript[ω, m]) #2)) &)};

Substituting the svRule into eqEps[[1]], we obtain the characteristic equations

eqEps1lhs = 
      Expand[(E^(-I m θ - I (Subscript[ω, m]) (Subscript[T, 0]))) #[
                  [1]] /. svRule] & /@ eqEps[[1]]

Setting the determinant of the coefficient matrix in eqEps1lhs equal to zero yields the following equation governing the natural frequencies (Subscript[ω, m]) of the shell:

(Outer[D, eqEps1lhs, {A, B}] // Det // Collect[#, (Subscript[ω, m])] &)
         == 0

It follows from eqEps1lhs that the ratio (Subscript[Γ, m]) = I  (B/A) of the amplitudes of (Subscript[ψ, 1]) and (Subscript[w, 1]) is given by

ΓRule = 
      (Subscript[Γ, m]) -> 
        Factor[I  (B/A) /. Solve[eqEps1lhs[[1]] == 0, B][[1]]]

Hence, the solution of eqEps[[1]] can be expressed in terms of the linear free-vibration modes

modeShapes = 
      ComplexExpand[{Re[#], Im[#]}] & /@
         {Exp[I n θ], 
          -I (Subscript[Γ, n]) Exp[I n θ]}

as

\(sol1Form = 
      modeShapes . 
          Table[A\_\(i, n\)[T\_1, T\_2]\ Exp[I\ ω\_n\ T\_0] + 
              A\&_\_\(i, n\)[T\_1, T\_2]\ Exp[\(-I\)\ ω\_n\ T\_0], {i, 
              2}] // Expand\)

It follows from the sol1Form that all of the flexural modes are degenerate because two linearly independent eigenfunctions, namely Cos[n θ] and Sin[n θ], correspond to the same frequency (Subscript[ω, n]) for each n.

We consider the case of a primary-resonance excitation of the Cos[n θ] mode corresponding to the frequency (Subscript[ω, n]); it is called the driven mode. Because of the degeneracy, this mode is involved in a one-to-one internal resonance with the Sin[n θ] mode, which is called the companion mode. We assume that neither of these two modes is involved in an internal resonance with any other mode. Therefore, in the presence of viscous damping, out of the infinite number of modes present in (Subscript[w, 1]) and (Subscript[ψ, 1]), only the driven mode Cos[n θ] and its companion mode Sin[n θ] will contribute to the steady-state response (Nayfeh and Mook, 1979). Consequently, we take the solution of the first-order problem to consist only of these two modes and express it as

sol1 = {
        (Subscript[w, 1]) -> Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1Form[[1]] // Evaluate], 
        (Subscript[ψ, 1]) -> 
          Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol1Form[[2]] // Evaluate]};

To describe the nearness of the primary resonance, we introduce the detuning parameter σ defined by

ResonanceCond = {
        Ω == (Subscript[ω, n]) + (ϵ^2) σ};

where the detuning is ordered at (ϵ^2) because secular terms appear first at order (ϵ^3).

8.2.2 Second-Order Solution

Substituting the sol1 into the second-order problem, eqEps[[2]], we obtain

order2Eq = 
      Table[(eqEps[[1, i, 1]] /. {(Subscript[w, 1]) -> (Subscript[w, 2]), (Subscript[ψ, 1]) -> (Subscript[ψ, 2])})
             == (
              eqEps[[1, i, 1]] /. {(Subscript[w, 1]) -> (Subscript[w, 2]), (Subscript[ψ, 1]) -> (Subscript[ψ, 2])})
               - (Subtract @@ eqEps[[2, i]] /. sol1 // Expand), {
          i, 2}];

Collecting the terms that may lead to secular terms, the terms proportional to (E^(I (Subscript[ω, n]) (Subscript[T, 0]))), we have

ST1 = 
      Coefficient[#[[2]], (E^(I (Subscript[ω, n]) (Subscript[T, 0])))] & /@ 
        order2Eq

Because the homogeneous parts of order2Eq have a nontrivial solution, the corresponding nonhomogeneous equations have a solution only if solvability conditions are satisfied. These conditions demand that ST1 be orthogonal to every solution of the adjoint homogeneous equations. In this case, these equations are self-adjoint and hence solutions of the adjoint equations are given by the transpose of the modeShapes. Demanding that ST1 be orthogonal to the transpose of the modeShapes, we obtain the solvability conditions

SCond1=Solve[(int[#, {θ, 0, 2  π}]& /@ 
                          (Transpose[modeShapes] . ST1) /. intRule1 //. 
                      intRule2 /. {int[Cos[s_] Sin[s_], _] -> 0, 
                      int[f_[n θ]^2, arg_] -> π}) == 0 // Thread
              ,{(Subscript[A, 1, n]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])],(Subscript[A, 2, n]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]}][[1]]

Having determined the solvability conditions of the second-order equations, order2Eq, we use a combination of the superposition principle and the method of undetermined coefficients to determine a particular solution of them. To this end, we identify the forms of all of the possible terms on the right-hand sides of order2Eq as follows:

\(collectForm = 
      Outer[Times, {Cos[n\ θ], Sin[n\ θ]}, 
          Table[{A\_\(i, n\)[T\_1, T\_2]\ E\^\(I\ ω\_n\ T\_0\), 
                A\&_\_\(i, n\)[T\_1, T\_2]\ 
                  E\^\(\(-I\)\ ω\_n\ T\_0\)}, {i, 2}] // Flatten] // 
        Flatten\)

\(possibleTerms = 
      Flatten[\(If[Head[#] === Plus, List\ @@\ #, #]\ &\)\ /@\ 
              \((\(\(\(Outer[Times, collectForm, collectForm] // Flatten\) // 
                      Union\) // TrigReduce\) // Expand)\)] /. 
          c_?NumberQ\ form_ -> form // Union; 

    possibleTerms /. displayRule\)

possibleTerms // Length

We note that the complex conjugate of each term on the right-hand sides of order2Eq is also a term on their right-hand sides. Therefore, we use this fact to reduce the algebra and calculate only the particular solutions corresponding to almost half of the terms because the particular solutions corresponding to the remaining terms are the complex conjugates of those already calculated. Then, we express the right-hand side of order2Eq[[i]] as the sum of two terms, tempList[i] and its complex conjugate, for i =  1 and 2, according to

\(Do[eq2rhsList = List\ @@\ order2Eq[\([i, 2]\)]; 
      tempList[i] = eq2rhsList; 
	\ \ 
      \(If[term = # /. conjugateRule; MemberQ[tempList[i], term], 
            tempList[i] = Complement[tempList[i], {#}]; 
            If[# === term, AppendTo[tempList[i], 1\/2\ #]]]\ &\)\ /@\ 
        Reverse[eq2rhsList], {i, 2}]\)

Next, we determine the coefficients of the possibleTerms in each of the tempList[i] according to

Do[rhs = Plus @@ tempList[i] /. SCond1 // TrigReduce // Expand; 
      coeffs[i] = 
        Coefficient[rhs, possibleTerms] /. 
          (Cos | Sin)[2  n θ] -> 0, {i, 2}]

and denote them by (Subscript[c, i,  j]) as

symbolList1 = 
      MapIndexed[If[#1 =!= 0, (Subscript[c, Sequence @@ #2]), 0] &, 
        Array[coeffs, 2], {2}]

symbolList2[s_] = Plus @@ symbolList1 /. c -> s;

where (Subscript[c, i,  j]) =  coeffs[i][[j]]. Moreover, we associate with the jth term in the symbolList1[[i]] the undetermined coefficient (Subscript[d, i,  j]) in (Subscript[w, 2]) and the undetermined coefficient (Subscript[e, i,  j]) in (Subscript[ψ, 2]) by replacing the (Subscript[c, i,  j]) in symbolList1[[i]] with (Subscript[d, i,  j]) and (Subscript[e, i,  j]) according to

sol2temp = {
        (Subscript[w, 2]) -> Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            symbolList2[d] . possibleTerms // Evaluate], 
        (Subscript[ψ, 2]) -> 
          Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            symbolList2[e] . possibleTerms // Evaluate]};

Substituting the sol2temp into the left-hand sides of the second-order equations, order2Eq, setting the coefficient of the jth term in possibleTerms in the ith equation equal to (Subscript[c, i,  j]), and solving the resulting algebraic equations, we obtain

deRule = 
      Solve[Flatten[
                  Table[Coefficient[order2Eq[[i, 1]] /. sol2temp, 
                        possibleTerms] /. f_[2  n θ] -> 0, {i, 
                      2}]] == Flatten[symbolList1] // Thread // Union, 
          {symbolList2[d], symbolList2[e]} // Flatten // Union // 
            Rest][[1]]

To simplify the expressions in the second-order solution, we introduce the notation

paperNotationRule = 
      Thread[(symbolList2[c] // Union // Rest) -> {(Subscript[c, 2]), (Subscript[c, 1]), 
              2  (Subscript[c, 1]), (Subscript[c, 2]), -(Subscript[c, 1]), -(Subscript[c, 4]), -(1/2) (Subscript[c, 3]), 
              -(Subscript[c, 3]), -(Subscript[c, 4]), (1/2) (Subscript[c, 3]), (Subscript[d, 1]), -2 (Subscript[d, 1]), 
              -(Subscript[d, 1]), -(1/2) (Subscript[d, 2]), (Subscript[d, 2]), (1/2) (Subscript[d, 2])}]~
        Join~{(
              1 + (α^2) - 8 (n^2) (α^2) + 
                16 (n^4) (α^2) - 4 (Subscript[ω, n]^2)) -> 
            (((Subscript[Δ, n]) + 4  (n^2))/(4 (n^2) - 
                4 (Subscript[ω, n]^2))), 
          Power[-1 + 4 (n^2), -2] Power[α, -2] -> 
            ((4 (n^2))/(Subscript[Θ, n]))}

Then, we express the second-order solution as

\(w2Sol = 
      \(symbolList2[d] . possibleTerms /. deRule\) //. paperNotationRule; 

    w2Sol /. displayRule\)

\(ψ2Sol = 
      \(\(symbolList2[e] . possibleTerms /. deRule\) //. paperNotationRule
          \) /. e\_i__ -> 0; 
ψ2Sol /. displayRule\)

We express this solution in pure function form as

sol2 = {
        (Subscript[w, 2]) -> Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            w2Sol + (w2Sol /. conjugateRule) // Evaluate], 
        (Subscript[ψ, 2]) -> 
          Function[{θ, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            ψ2Sol + (ψ2Sol /. conjugateRule) // Evaluate]};

8.2.3 Solvability Conditions

Substituting the sol1 and sol2 into eqEps[[3]] yields

order3Eq = 
        Table[(
              eqEps[[1, i, 1]] /. {(Subscript[w, 1]) -> (Subscript[w, 3]), (Subscript[ψ, 1]) -> (Subscript[ψ, 3])})
               == (
                eqEps[[1, i, 1]] /. {(Subscript[w, 1]) -> (Subscript[w, 3]), 
                    (Subscript[ψ, 1]) -> (Subscript[ψ, 3])}) - 
              (Subtract @@ eqEps[[3, i]] /. sol1 /. sol2 // 
                    TrigToExp // Expand), {i, 2}];  // Timing

To convert the terms that lead to small-divisor terms in the third-order equations, order3Eq, into terms that lead to secular terms, we introduce the rule

expRule1 = 
      Exp[arg_] :> 
        Exp[Expand[arg /. (ResonanceCond /. Equal -> Rule)] /. 
            (ϵ^2) (Subscript[T, 0]) -> (Subscript[T, 2])];

Then, the sources of secular terms in order3Eq, the terms proportional to \(TraditionalForm
      \`Exp[\(±\ I\)\ n\ θ + I\ ω\_n\ T\_0]\), are

ST2 = 
        Expand[
              Coefficient[
                  #[[2]] /. 
                      Exp[c1_.  + 
                            Complex[0, s_ /; Abs[s] != 1] n θ] -> 0
                         /. expRule1, Exp[I (Subscript[ω, n]) (Subscript[T, 0])]] /. 
                Exp[c1_.  +  Complex[0, m_] n θ] -> 
                  (Cos[n θ] + I m Sin[n θ]) Exp[c1]] &
             /@ order3Eq;  // Timing

These terms consist of linear and cubic terms in the complex-valued amplitudes (Subscript[A, i]). The cubic terms are proportional to

\(cubicTerms = 
      \(\(\(Flatten[
                  \(If[Head[#] === Plus, List\ @@\ #, #]\ &\)\ /@\ 
                    \((\(\(\(Nest[Outer[Times, collectForm, #]\ &, 
                                collectForm, 2] // Flatten\) // Union\) // 
                          TrigReduce\) // Expand)\)] /. {
                  f_[s_\ n\ θ] -> 0, 
                  Exp[c1_.  + Complex[0, s_ /; s != 1]\ ω\_n\ T\_0] -> 
                    0}\) /. Exp[_] -> 1\) /. c_?NumberQ\ form_ -> form // 
          Union\) // Rest; 
cubicTerms /. displayRule\)

And, their coefficients in ST2 are

coef2 = Outer[Coefficient, ST2, cubicTerms];

which we denote by

symbolList3 = 
      MapIndexed[If[#1 =!= 0, (Subscript[g, Sequence @@ #2]), 0] &, coef2, {2}];

notationRule = 
      (symbolList3 // Flatten // Union // Rest) -> {(Subscript[E, 2]), (Subscript[E, 1]), 
            (Subscript[E, 3]), (Subscript[E, 3]), (Subscript[E, 1]), (Subscript[E, 2]), -(Subscript[G, 2]), (Subscript[G, 1]), -(Subscript[G, 3]), (Subscript[G, 3]), 
            -(Subscript[G, 1]), (Subscript[G, 2])} // Thread

In terms of this notation, ST2 can be rewritten as

\(ST2New = 
      \(\((\(# /. Thread[cubicTerms -> 0]\ &\)\ /@\ ST2)\) + 
            symbolList3 . cubicTerms /. notationRule\) /. D[SCond1, T\_1]; 

    ST2New /. displayRule\)

The solvability conditions of order3Eq demand that the vector ST2New be orthogonal to the transpose of the modeShapes; that is,

\(SCond2 = 
      Expand[\(\(\(1\/π\ int[#, {θ, 0, 2  π}]\ &\)\ /@\ 
                    \((Transpose[modeShapes] . ST2New)\) /. intRule1\) //. 
                intRule2\) /. {int[Cos[s_]\ Sin[s_], _] -> 0, 
                int[f_[n\ θ]^2, _] -> π}] == 0 // Thread; 

    SCond2 /. displayRule\)

8.3 Near-Square Plates

We consider the nonlinear transverse vibrations of a near-square simply supported plate that is initially flat. We introduce a Cartesian coordinate system x y z, with the x and y axes lying in the undisturbed plane of the middle surface of the plate, and let u, v, and w denote the displacements of the middle surface in the x, y, and z directions, respectively. Then, the equations of motion are (Chu and Herrmann, 1956)

Needs["Calculus`VectorAnalysis`"]

BH[w] = Biharmonic[w[x, y, t], Cartesian[x, y, z]]

\(\(EOM = {
		
          u\_\(x, x\) + w\_x\ w\_\(x, x\) + 
              ν\ \((v\_\(x, y\) + w\_y\ w\_\(x, y\))\) + 
              \(1\/2\) \((1 - ν)\)\ 
                \((u\_\(y, y\) + v\_\(x, y\) + w\_x\ w\_\(y, y\) + 
                    w\_y\ w\_\(x, y\))\) == c\_p\%\(-2\)\ u\_\(t, t\), 
		
          v\_\(y, y\) + w\_y\ w\_\(y, y\) + 
              ν\ \((u\_\(x, y\) + w\_x\ w\_\(x, y\))\) + 
              \(1\/2\) \((1 - ν)\) 
                \((v\_\(x, x\) + u\_\(x, y\) + w\_y\ w\_\(x, x\) + 
                    w\_x\ w\_\(x, y\))\) == c\_p\%\(-2\)\ v\_\(t, t\), 
		
          \(1\/12\) h\^2\ BH[w] - N\_1\%i\ w\_\(x, x\) - 
              2  N\_12\%i\ w\_\(x, y\) - N\_2\%i\ w\_\(y, y\) - 
              u\_x\ w\_\(x, x\) - \(1\/2\) w\_x\%2\ w\_\(x, x\) - 
              v\_y\ w\_\(y, y\) - \(1\/2\) w\_y\%2\ w\_\(y, y\) - 
              ν \((
                  v\_y\ w\_\(x, x\) + \(1\/2\) w\_y\%2\ w\_\(x, x\) + 
                    u\_x\ w\_\(y, y\) + \(1\/2\) w\_x\%2\ w\_\(y, y\))\) - 
              \((1 - ν)\) \((u\_y + v\_x + w\_x\ w\_y)\) w\_\(x, y\) == 
            \(-c\_p\%\(-2\)\) 
                \((w\_x\ u\_\(t, t\) + w\_y\ v\_\(t, t\) + w\_\(t, t\))\) + 
              \(F \((1 - ν\^2)\)\)\/\(E\_s\ h\)
	} /. 
        \((s : u | v | w)\)\_m__ :> ∂\_m s[x, y, t]; \)\)

and the boundary conditions are

BC = {{u[x, y, t] == 0, w[x, y, t] == 0, 
              (∂_(x, x))w[x, y, t] == 0} /. {{x -> 0}, {x -> a}}, {
              v[x, y, t] == 0, w[x, y, t] == 0, 
              (∂_(y, y))w[x, y, t] == 0} /. {{y -> 0}, {
                y -> b}}} // Flatten

where (Subscript[c, p]^2) = (Subscript[E, s])/(ρ(1 - (ν^2))), (Subscript[E, s]) is Young's modulus, ν is Poisson's ratio, ρ and h are the density and thickness of the plate, respectively, F is the transverse load, the (Subscript[N, k]^i) are the uniformly distributed in-plane applied edge loads, and a and b are the edge lengths, with b ≈ a. EOM are usually called the dynamic analog of the von Karman equations because they reduce to von Karman's equations (1910) in the absence of the time derivatives.

For thin plates, h/a is very small, which we denote by (Sqrt[ϵ]). Next, we introduce the aspect ratio c = a/b, define the rule

scaleRule1 = {(h^2) -> ϵ (a^2), a -> c b};

and nondimensionalize the variables as

charT = (Sqrt[12/ϵ]) a/(Subscript[c, p]);

nondimRule = {
        (dv : u | v | w) -> (a dv[#1/a, #2/b, #3/charT] &), 
        F -> F (Subscript[E, s]) h/(12  a (1 - (ν^2))), x -> a x, y -> b y, 
        t -> charT t};

Substituting nondimRule and scaleRule1 into EOM and BC and considering the case of no in-plane edge loads, we obtain

\(\(EOM1 = 
      \(\((\(\(a \((Subtract\ @@\ #)\) /. N\_k_\%i -> 0\) /. nondimRule\) //. 
                  scaleRule1 // Expand)\) == 0\ &\)\ /@\ EOM; \)\)

BC1 = BC /. nondimRule /. c_ h_[a__, t] :> h[a, t]

To write EOM1 in a more concise way, we introduce the display rule

display1[expr_] := 
      expr /. {Derivative[a_, b_, c_][w_][x, y, t] :> 
            Subscript[w, 
              Sequence @@ 
                Join[Table[x, {a}], Table[y, {b}], Table[t, {c}]]]}

and obtain

EOM1[[1]] // display1EOM1[[2]] // display1EOM1[[3]] // display1

We use the method of multiple scales to directly attack EOM1 and BC1. To transform the time derivatives in EOM1 in terms of the scales (Subscript[T, 0]) and (Subscript[T, 1]), we let

multiScales = {
        (dv : u | v | w)[x_, y_, t] -> dv[x, y, (Subscript[T, 0]), (Subscript[T, 1])], 
        Derivative[m_, n_, o_][u_][x_, y_, t] -> 
          dt[o][D[u[x, y, (Subscript[T, 0]), (Subscript[T, 1])], {x, m}, {y, n}]]};

Then, we seek a second-order approximate solution in the form

solRule = {
        (dv : u | v) -> 
          (Evaluate[
                Sum[(ϵ^(j + 1)) (Subscript[dv, j])[#1, #2, #3, #4], {j, 2}]] &)
            , w -> 
          (Evaluate[Sum[(ϵ^j) (Subscript[w, j])[#1, #2, #3, #4], {j, 3}]] &)}

where #1, #2, #3, and #4 stand for x, y, (Subscript[T, 0]),  and (Subscript[T, 1]), respectively.

We introduce the detuning (Subscript[σ, 1]) to describe the nearness of the edge lengths a and b, consider the case of primary resonance, and hence let

scaleRule2 = {c -> 1 + ϵ (Subscript[σ, 1]), 
        F -> (ϵ^3) F[x, y] Cos[Ω (Subscript[T, 0])]};

Substituting multiScales, solRule, and scaleRule2 into EOM1 and BC1, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq83a = 
      (EOM1 /. multiScales /. solRule /. scaleRule2 // TrigToExp
               // ExpandAll) /. (ϵ^(n_ /; n > 3)) -> 0;

bc83a = BC1 /. multiScales /. solRule;

Equating coefficients of like powers of ϵ in eq83a and bc83a yields

eqEps = 
      Thread[CoefficientList[Subtract @@ #, ϵ] == 0] & /@
               eq83a // Transpose // Rest // Rest;

bcEpsuv = 
      Thread[CoefficientList[Subtract @@ #, ϵ] == 0] & /@
               Select[bc83a, FreeQ[#, w] &] // Transpose // Rest // Rest

bcEpsw = 
      Thread[CoefficientList[Subtract @@ #, ϵ] == 0] & /@ 
            Select[bc83a, ! FreeQ[#, w] &] // Transpose // Rest

Next, we introduce the display rule

\(\(displayRule1 = {
        \(\(Derivative[a_, b_, c__]\)[w_\_i_]\)[x_, y_, T\_0, __] :> 
          SequenceForm[
            If[arg1 = 
                Times\ @@\ MapIndexed[D\_\(#2[\([1]\)] - 1\)\%#1\ &, {c}]; 
              arg1 =!= 1, arg1, ""], 
            Subscript[w, 
              Sequence\ @@\ Join[{i}, Table[x, {a}], Table[y, {b}]]]], 
        \(\(Derivative[a__]\)[A\_i__]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\_i], \(\(Derivative[a__]\)[A\&_\_i__]\)[__] :> 
          SequenceForm[Times\ @@\ MapIndexed[D\_\(#2[\([1]\)]\)\%#1\ &, {a}], 
            A\&_\_i], w_\_i_[__, T\_0, __] -> w\_i, A\_i__[__] -> A\_i, 
        A\&_\_i__[__] -> A\&_\_i}; \)\)

and write eqEps in a more concise way as

(eqEpsm = {eqEps[[1, 3]], eqEps[[1, {1, 2}]], 
            eqEps[[2, 3]]}) /. displayRule1

8.3.1 First-Order Solution

Linear vibrations are governed by eqEpsm[[1]] and bcEpsw[[1]]; that is,

(order1Eq = eqEpsm[[1]]) /. displayRule1

order1BC = bcEpsw[[1]]

It follows from order1Eq and order1BC that the linear modeshapes and associated frequencies are given by

\(modeshapes = {
        ϕ\_\(m_, n_\) :> \((Sin[m\ π\ #1]\ Sin[n\ π\ #2]\ &)\)}; 
    
frequencies = {ω\_\(m_, n_\) :> \((m\^2 + n\^2)\) π\^2}; 
    \)

To investigate the case of one-to-one internal resonance between the mnth and nmth modes (m ≠ n), we let

\(\(w1sol = {
        w\_1 -> Function[{x, y, T\_0, T\_1}, 
            A\_1[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\)\ 
                ϕ\_\(m, n\)[x, y] + 
              A\&_\_1[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\)\ 
                ϕ\_\(m, n\)[x, y] + 
              A\_2[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\)\ 
                ϕ\_\(n, m\)[x, y] + 
              A\&_\_2[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\)\ 
                ϕ\_\(n, m\)[x, y]]}; \)\)

8.3.2 Second-Order Solution

Substituting the w1sol into eqEpsm[[2]], we obtain

order2Equ = 
        (eqEpsm[[2, 1, 1]] /. (Subscript[w, 1]) -> (0 &)) == 
          ((eqEpsm[[2, 1, 1]] /. (Subscript[w, 1]) -> (0 &)) - 
                      eqEpsm[[2, 1, 1]] /. w1sol /. modeshapes // 
                TrigReduce // Expand);  // Timing

order2Eqv = 
        (eqEpsm[[2, 2, 1]] /. (Subscript[w, 1]) -> (0 &)) == 
          ((eqEpsm[[2, 2, 1]] /. (Subscript[w, 1]) -> (0 &)) - 
                      eqEpsm[[2, 2, 1]] /. w1sol /. modeshapes // 
                TrigReduce // Expand);  // Timing

The associated boundary conditions are

order2BC = bcEpsuv[[1]]

To obtain particular solutions for (Subscript[u, 1]) and (Subscript[v, 1]), we use the method of undetermined coefficients. To accomplish this, we first look for all possible terms that appear on the right-hand sides of order2Equ and order2Eqv as follows:

possibleTerms1 = 
      Cases[{order2Equ[[2]], order2Eqv[[2]]}, _ Sin[a_] (E^b)_ :> 
            Sin[a] (E^b), ∞] // Union;

possibleTerms2 = possibleTerms1 /. Exp[_] -> 1 // Union;

possibleTerms = Join[possibleTerms1, possibleTerms2]

which satisfy order2BC. We then define the following symbol lists to represent the undetermined coefficients.

usymbolList = Table[(Subscript[c, i]), {i, Length[possibleTerms]}]

vsymbolList = usymbolList /. c -> d

Then, we express the solution as

\(u1sol = 
      u\_1 -> Function[{x, y, T\_0, T\_1}, 
          usymbolList . possibleTerms // Evaluate]; 

    v1sol = v\_1 -> 
        Function[{x, y, T\_0, T\_1}, 
          vsymbolList . possibleTerms // Evaluate]; \)

Substituting u1sol and v1sol into order2Equ and order2Eqv, collecting coefficients of possibleTerms, and solving for the usymbolList and vsymbolList, we obtain

algEq1 = 
        Coefficient[Subtract @@ order2Equ /. u1sol /. v1sol, 
                possibleTerms] == 0 /. Exp[_ (Subscript[T, 0])] -> 0 // Thread;  // 
      Timing

algEq2 = 
        Coefficient[Subtract @@ order2Eqv /. u1sol /. v1sol, 
                possibleTerms] == 0 /. Exp[_ (Subscript[T, 0])] -> 0 // Thread;  // 
      Timing

symbolRule1 = 
      Solve[Join[algEq1, algEq2], Join[usymbolList, vsymbolList]][
          [1]] // Simplify

8.3.3 Solvability Conditions

Substituting the w1sol, u1sol, and v1sol into eqEpsm[[3]], we obtain

order3Eq = 
        (order1Eq[[1]] /. (Subscript[w, 1]) -> (Subscript[w, 2])) == 
          ((order1Eq[[1]] /. (Subscript[w, 1]) -> (Subscript[w, 2])) - 
                        eqEpsm[[3, 1]] /. w1sol /. u1sol /. v1sol
                   /. modeshapes // Expand);  // Timing

The associated boundary conditions are

order3BC = bcEpsw[[2]]

To describe the nearness of the primary resonance, we introduce the detuning parameter (Subscript[σ, 2]) defined by

ResonanceCond = {
        Ω -> (Subscript[ω, m, n]) + ϵ (Subscript[σ, 2])};

Because the homogeneous parts of order3Eq and order3BC have a nontrivial solution, the corresponding nonhomogeneous problem has a solution only if solvability conditions are satisfied. These conditions demand that the right-hand side of order3Eq be orthogonal to every solution of the adjoint homogeneous problem. In this case, the problem is self-adjoint.

To determine the solvability conditions, we collect the terms that may lead to secular terms, the terms proportional to (E^(I (Subscript[ω, m, n]) (Subscript[T, 0]))). To this end, we use the rule

expRule = 
      Exp[a_] :> Exp[Expand[a /. ResonanceCond] /. ϵ (Subscript[T, 0]) -> (Subscript[T, 1])];

and obtain

ST = 
        Coefficient[order3Eq[[2]] /. expRule, 
            (E^(I (Subscript[ω, m, n]) (Subscript[T, 0])))] // TrigReduce;  // Timing

Among all the spatial-dependent non-forcing terms in ST, only the following terms are not orthogonal to the modeshapes:

terms1 = 
      Cases[ST, Cos[m π x + a_.  n π y /; Abs[a] == 1], 
          ∞] // Union

terms2 = 
      Cases[ST, Cos[n π x + a_.  m π y /; Abs[a] == 1], 
          ∞] // Union

The coefficients of these terms in ST are

terms1Coef = Coefficient[ST, terms1];

terms2Coef = Coefficient[ST, terms2];

The projection of terms1 onto Sin[m π x] Sin[n π y] can be obtained as

trigRule = {Cos[2 _.  π] -> 1, Sin[_.  π] -> 0};

proj1Coef = 
      (∫_0^1)
                (∫_0^1) Sin[m π x] Sin[n π y] # 
                  \[DifferentialD]x \[DifferentialD]y & /@ terms1 /. 
        trigRule

and the projection of terms2 onto Sin[n π x] Sin[m π y] can be obtained as

proj2Coef = 
      (∫_0^1)
                (∫_0^1) Sin[n π x] Sin[m π y] # 
                  \[DifferentialD]x \[DifferentialD]y & /@ terms2 /. 
        trigRule

The forcing term in ST is

f0 = ST /. Cos[_] -> 0

Then, the solvability conditions demand that ST be orthogonal to every linear eigenfunction, modeshapes; that is,

intRule3 = {
        int[a_ fun_, a1_, b1_] :> 
          a int[fun, a1, b1] /; FreeQ[a, First[a1] | First[b1]]};

SCond1 = 
      Expand[24 
            (int[Sin[m π x] Sin[n π y] f0, {x, 0, 1}, {y, 
                          0, 1}] + proj1Coef . terms1Coef /. symbolRule1
                     //. intRule3 /. int -> Integrate)] == 0;

SCond2 = 
      Expand[24 
            (int[Sin[n π x] Sin[m π y] f0, {x, 0, 1}, {y, 
                          0, 1}] + proj2Coef . terms2Coef /. symbolRule1
                     //. intRule3 /. int -> Integrate)] == 0;

To express these solvability conditions in a more readable form, we define

\(\(basicTerms = {A\_1[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\), 
        A\&_\_1[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\), 
        A\_2[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\), 
        A\&_\_2[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\)}; \)\)

cubicTerms = 
      ((E^(-I (Subscript[ω, m, n]) (Subscript[T, 0]))) 
                  Nest[Outer[Times, basicTerms, #] &, basicTerms, 2] // 
                Flatten) /. Exp[_ (Subscript[T, 0])] -> 0 // Union // Rest

Then, we collect the coefficients of cubicTerms from SCond1 and SCond2 and obtain

coef1 = Coefficient[SCond1[[1]], cubicTerms] // Simplify

coef2 = Coefficient[SCond2[[1]], cubicTerms] // Simplify

Hence, SCond1 and SCond2 can be rewritten and combined as

\(SCond = {
        Collect[SCond1[\([1]\)] /. Thread[cubicTerms -> 0], 
              σ\_1\ A\_1[T\_1]] + coef1 . cubicTerms == 0, 
        Collect[SCond2[\([1]\)] /. Thread[cubicTerms -> 0], 
              σ\_1\ A\_2[T\_1]] + coef2 . cubicTerms == 0}; 

    SCond /. displayRule1\)

8.3.4 Mixed Approach

Alternatively, we can use a mixed approach, a combination of discretization and direct approaches, to attack the problem, EOM1 and BC1. We first define

scaleRule3 = {
        (dv : u | v) -> ((ϵ^2) dv[#1, #2, #3] &), 
        w -> (ϵ w[#1, #2, #3] &)};

Substituting scaleRule2 and scaleRule3 into EOM1, expanding the result for small ϵ, and discarding terms of order higher than (ϵ^3), we obtain

eq83b = 
      (EOM1 /. scaleRule2 /. scaleRule3 // TrigToExp // ExpandAll)
           /. (ϵ^(n_ /; n > 3)) -> 0;

To write eq83b in a more concise way, we introduce the dispaly rule

display2[expr_] := 
      expr /. {Derivative[a_, b_, c_][w_][x, y, t] :> 
            Subscript[w, 
              Sequence @@ 
                Join[Table[x, {a}], Table[y, {b}], Table[t, {c}]]]}

Hence,

(uEOM = Coefficient[eq83b[[1, 1]], (ϵ^2)] == 0) // 
      display2

(vEOM = Coefficient[eq83b[[2, 1]], (ϵ^2)] == 0) // 
      display2

(wEOM = 
          eq83b[[3]] // Thread[(ϵ^(-2)) #, Equal] & // 
            ExpandAll) // display2

Linear vibrations are governed by wEOM with ϵ = 0; that is,

(order1Eqw = wEOM /. ϵ -> 0) // display2

subject to the boundary conditions BC1. It follows from order1Eqw and BC1 that the linear modeshapes and associated frequencies are given by

(Subscript[ϕ, m, n])[x, y] /. modeshapes

(Subscript[ω, m, n]) /. frequencies

To investigate the case of one-to-one internal resonance between the mnth and nmth modes (m ≠ n), we let

w1sol = {
        w -> Function[{x, y, t}, 
            η[t] Sin[m π x] Sin[n π y] + 
              ζ[t] Sin[n π x] Sin[m π y]]};

Substituting the w1sol into uEOM and vEOM, we obtain

\(order2Equ = 
      \((uEOM[\([1]\)] /. w -> \((0\ &)\))\) == 
        \((\((uEOM[\([1]\)] /. w -> \((0\ &)\))\) - uEOM[\([1]\)] /. w1sol // 
            TrigReduce)\); 
order2Equ // display2\)

\(order2Eqv = 
      \((vEOM[\([1]\)] /. w -> \((0\ &)\))\) == 
        \((\((vEOM[\([1]\)] /. w -> \((0\ &)\))\) - vEOM[\([1]\)] /. w1sol // 
            TrigReduce)\); 
order2Eqv // display2\)

To obtain the particular solutions for u and v, we use the method of undetermined coefficients. To accomplish this, we first look for all possible terms that appear on the right-hand sides of order2Equ and order2Eqv as follows:

possibleTerms = 
      Cases[{order2Equ[[2]], order2Eqv[[2]]}, #, ∞] &
               /@ {_ Sin[a_] h1_[t]^2 -> 
                Sin[a] h1[t]^2, _ Sin[a_] h1_[t] h2_[t] -> 
                Sin[a] h1[t] h2[t]} // Flatten // Union

which satisfy BC1. We then define the following symbol lists to represent the undetermined coefficients.

usymbolList = Table[(Subscript[c, i]), {i, Length[possibleTerms]}]

vsymbolList = usymbolList /. c -> d

Then, we express the solution as

\(usol = 
      u -> Function[{x, y, t}, usymbolList . possibleTerms // Evaluate]; 

    vsol = v -> Function[{x, y, t}, vsymbolList . possibleTerms // Evaluate]; 
    \)

Substituting usol and vsol into order2Equ and order2Eqv, collecting coefficients of possibleTerms, and solving for usymbolList and vsymbolList, we obtain

algEq1 = 
      Coefficient[Subtract @@ order2Equ /. usol /. vsol, 
            possibleTerms] == 0 // Thread;

algEq2 = 
      Coefficient[Subtract @@ order2Eqv /. usol /. vsol, 
            possibleTerms] == 0 // Thread;

(symbolRule2 = 
          Solve[Join[algEq1, algEq2], Join[usymbolList, vsymbolList]][
              [1]] // Simplify) // Timing

Solvability Conditions

We introduce two time scales (Subscript[T, 0]) and (Subscript[T, 1]) and seek a first-order uniform asymptotic expansion of the solution of wEOM in the form

solRule = {
        w -> Function[{x, y, (Subscript[T, 0]), (Subscript[T, 1])}, 
            Sum[(ϵ^(j - 1)) (Subscript[w, j])[x, y, (Subscript[T, 0]), (Subscript[T, 1])], {j, 2}]]};

Transforming the time derivatives in wEOM in terms of the scales (Subscript[T, 0]) and (Subscript[T, 1]), substituting solRule into wEOM, expanding the result for small ϵ, and discarding terms of order higher than ϵ, we obtain

eq83c = 
      (wEOM /. {(s : u | v | w)[x, y, t] -> s[x, y, (Subscript[T, 0]), (Subscript[T, 1])], 
                  Derivative[a_, b_, c_][h_][x, y, t] -> 
                    dt[c][D[h[x, y, (Subscript[T, 0]), (Subscript[T, 1])], {x, a}, {y, b}]]} /. 
              solRule // ExpandAll) /. (ϵ^(n_ /; n > 1)) -> 0;

Equating coefficients of like powers of ϵ, we have

eqEps = CoefficientList[eq83c[[1]], ϵ] == 0 // Thread;

Accounting for the two time scales (Subscript[T, 0]) and (Subscript[T, 1]), we create a rule for η[t] and ζ[t]

\(\(ampRule = {
        η[t] -> 
          A\_1[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\) + 
            A\&_\_1[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\), 
        ζ[t] -> 
          A\_2[T\_1]\ E\^\(I\ ω\_\(m, n\)\ T\_0\) + 
            A\&_\_2[T\_1]\ E\^\(\(-I\)\ ω\_\(m, n\)\ T\_0\)}; \)\)

and rewrite w1sol, usol, and vsol as

sol1 = {
        (Subscript[w, 1]) -> Function[{x, y, (Subscript[T, 0]), (Subscript[T, 1])}, 
            w[x, y, t] /. w1sol /. ampRule // Evaluate], 
        u -> Function[{x, y, (Subscript[T, 0]), (Subscript[T, 1])}, 
            u[x, y, t] /. usol /. ampRule // Evaluate], 
        v -> Function[{x, y, (Subscript[T, 0]), (Subscript[T, 1])}, 
            v[x, y, t] /. vsol /. ampRule // Evaluate]};

Substituting the sol1 into the second-order problem, eqEps[[2]], we obtain

order2Eqw = 
        (eqEps[[1, 1]] /. (Subscript[w, 1]) -> (Subscript[w, 2])) == 
              (eqEps[[1, 1]] /. (Subscript[w, 1]) -> (Subscript[w, 2])) - eqEps[[2, 1]] /. 
            sol1 // ExpandAll;  // Timing

Collecting the terms that may lead to secular terms, the terms proportional to (E^(I (Subscript[ω, m, n]) (Subscript[T, 0]))), we have

ST = 
        Coefficient[order2Eqw[[2]] /. expRule, 
            (E^(I (Subscript[ω, m, n]) (Subscript[T, 0])))] // TrigReduce;  // Timing

Among all the spatial-dependent non-forcing terms in ST, only the following terms are not orthogonal to the modeshapes:

terms1 = 
      Cases[ST, Cos[m π x + a_.  n π y /; Abs[a] == 1], 
          ∞] // Union

terms2 = 
      Cases[ST, Cos[n π x + a_.  m π y /; Abs[a] == 1], 
          ∞] // Union

The coefficients of these terms in ST are

terms1Coef = Coefficient[ST, terms1];

terms2Coef = Coefficient[ST, terms2];

The projection of terms1 onto Sin[m π x] Sin[n π y] can be obtained as

p1Coef = 
      (∫_0^1)
                (∫_0^1) Sin[m π x] Sin[n π y] # 
                  \[DifferentialD]x \[DifferentialD]y & /@ terms1 /. 
        trigRule

and the projection of terms2 onto Sin[n π x] Sin[m π y] can be obtained as

p2Coef = 
      (∫_0^1)
                (∫_0^1) Sin[n π x] Sin[m π y] # 
                  \[DifferentialD]x \[DifferentialD]y & /@ terms2 /. 
        trigRule

The forcing term in ST is

f0 = ST /. Cos[_] -> 0

Then, the solvability conditions demand that ST be orthogonal to every linear eigenfunction, modeshapes; that is,

SC1 = 
      Expand[24 
            (int[Sin[m π x] Sin[n π y] f0, {x, 0, 1}, {y, 
                          0, 1}] + p1Coef . terms1Coef /. symbolRule2 //. 
                  intRule3 /. int -> Integrate)] == 0;

SC2 = 
      Expand[24 
            (int[Sin[n π x] Sin[m π y] f0, {x, 0, 1}, {y, 
                          0, 1}] + p2Coef . terms2Coef /. symbolRule2 //. 
                  intRule3 /. int -> Integrate)] == 0;

They agree with the results obtained in the previous section:

{SC1, SC2} === {SCond1, SCond2}