Chapter 9
Higher Approximations of Continuous Systems Having Two-to-One Internal Resonances

In this chapter, we use the method of multiple scales to determine second-order uniform asymptotic expansions of the solutions of continuous systems with quadratic and cubic nonlinearities possessing two-to-one internal resonances. Higher-order treatments of such systems lead to inconsistent results if the time derivatives in their governing equations are expressed in second-order rather than first-order form. Therefore, we express the time derivatives in the governing equations in first-order form before treating them with the method of multiple scales. We describe the methodology using two examples: two-mode interactions in hinged-hinged buckled beams and four-mode interactions in suspended cables.

Preliminaries

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

To determine second-order uniform asymptotic expansions of the solutions of  buckled beams and suspended cables by using the method of multiple scales, taking into account two-to-one internal resonances, 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

Needs["Utilities`Notation`"]

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}];

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]]}; \)

In the ensuing analysis, we express the solutions in terms of complex-variable functions and hence we introduce the complex-conjugate rule

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

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]}; \)\)

9.1 Two-Mode Interactions in Buckled Beams

In this section, we consider the nonlinear response of a hinged-hinged buckled beam possessing a two-to-one internal resonance to a principal parametric resonance of the higher mode. The analysis assumes a unimodal static buckled deflection. The nondimensional equation of motion is (Section 8.1.1)

EOM=(Subscript[u, t, t])+(Subscript[u, x, x, x, x])+(Subscript[P, k]) (Subscript[u, x, x])-2(b^2) α ((∫_0^1)(Subscript[u, x]) (Subscript[ϕ, k]^′)[x]\[DifferentialD]x) (Subscript[ϕ, k]^″)[x]==-c (Subscript[u, t])+b α ((∫_0^1) (Subscript[u, x]^2) \[DifferentialD]x) (Subscript[ϕ, k]^″)[x]+2b α ((∫_0^1)(Subscript[u, x]) (Subscript[ϕ, k]^′)[x]\[DifferentialD]x) (Subscript[u, x, x])+α 
                ((∫_0^1) (Subscript[u, x]^2) \[DifferentialD]x) 
                (Subscript[u, x, x])+Cos[t Ω] f (Subscript[u, x, x])+b Cos[t Ω] f (Subscript[ϕ, k]^″)[x]/.(Subscript[u, m])__ :> (∂_m) u[x, t];

and the associated boundary conditions are

BC={u[0, t] == 0,Derivative[2, 0][u][0, t]==0,u[1, t] == 0,Derivative[2, 0][u][1, t]==0};

where (Subscript[P, k]) is the kth critical Euler buckling load, (Subscript[ϕ, k]) is the kth buckling mode shape, and b is the nondimensional buckling level.

Because the two-to-one internal resonance is activated by the quadratic nonlinearities, as indicated in Chapter 5, second- and higher-order treatments of such systems lead to inconsistent results if the time derivatives in their governing equations are expressed in second-order rather than first-order form. Therefore, we express the time derivatives in first-order form by letting v[x, t] == (∂_t) u[x, t], adding it to EOM, and obtaining

EOM1 = {(∂_t) u[x, t] - v[x, t] == 0, 
        EOM /. (∂_(t, t))u[x, t] -> (∂_t) v[x, t]};

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

multiScales = {u[x_, t] -> u[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], 
        v[x, t] -> v[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 of EOM1 and BC in the form

solRule = 
      h : u | v -> 
        (Evaluate[Sum[(ϵ^j) (Subscript[h, 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.

We consider the case of a two-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 == 2  ω\_1 + ϵ\ σ\_1, 
        Ω == 2  ω\_2 + ϵ\ σ\_2}; \)

where the (Subscript[σ, i]) are detuning parameters that describe the nearness of the resonances. Because the influence of the two-to-one internal resonance appears at O(ϵ), we scale the damping and forcing terms as

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

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

eq91a = 
      (EOM1 /. Integrate -> int /. multiScales /. solRule
                       /. scaling /. intRule1 //. intRule2 // TrigToExp
               // ExpandAll) /. (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ, we obtain

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

9.1.1 First-Order Solution

The homogeneous part of the first-order equations, eqEps[[1]], can be written as

\(linearSys = \(#[\([1]\)]\ &\)\ /@\ eqEps[\([1]\)] /. f -> 0; 

    linearSys /. displayRule\)

Because in the presence of damping all of the modes that are not directly or indirectly excited decay with time (Nayfeh and Mook, 1979), the solution of eqEps[[1]] is taken to consist of the lowest two modes; that is,

\(\(solForm = 
      Ψ[x]\ Exp[I\ Ω\ T\_0] + 
        Ψ[x]\ Exp[\(-I\)\ Ω\ T\_0] + 
        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}]; \)\)

sol1 = {(Subscript[u, 1]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, solForm // Evaluate], 
        (Subscript[v, 1]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, D[solForm, (Subscript[T, 0])] // Evaluate]}

where the (Subscript[Φ, i]) are the eigenmodes, which satisfy the orthonormality condition \(TraditionalForm
      \`∫\_0\%1 Φ\_n[x]\ Φ\_m[x]\ 
            dx = δ\_\(n, m\)\), and Ψ is uniquely determined by the following boundary-value problem:

\(eq91b = 
      Coefficient[
          \(\(eqEps[\([1, 2, 1]\)] /. 
                \((sol1 /. {A\_i_ -> \((0\ &)\), A\&_\_i_ -> \((0\ &)\)})
                  \)\) /. intRule1\) //. intRule2, 
          E\^\(I\ Ω\ T\_0\)] == 0\)

bc91b = BC /. u -> Ψ

In what follows, we consider nonlinear vibrations of a first-mode buckled beam; hence, the critical buckling load is

bucklingLoad = (Subscript[P, k]) -> (π^2);

We note that, as the buckling level increases from zero, the symmetric mode is the first mode. However, as the buckling level exceeds the first cross-over value, the symmetric mode becomes the second mode. Hence, the second natural frequency (Subscript[ω, 2]) depends on the buckling level. For buckling levels greater than the cross-over value, the first static buckled deflection and the first two normalized mode shapes and associated natural frequencies are

\(modeshapes = {ϕ\_k -> \((Sin[π\ #]\ &)\), 
        Φ\_1 -> \((\(\@2\) Sin[2  π\ #]\ &)\), 
        Φ\_2 -> \((\(\@2\) Sin[π\ #]\ &)\)}; 

    frequencies = {ω\_2 -> b\ \(π\^2\) \@α, 
        ω\_1 -> 2 \@ 3\ π\^2}; \)

The buckling level at which (Subscript[ω, 2]) = 2  (Subscript[ω, 1]) is

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

Substituting the bucklingLoad and modeshapes into eq91b yields

eq91c = eq91b /. bucklingLoad /. modeshapes

Inspection of eq91c and bc91b suggests that Ψ[x] has the form

Ψsol = 
      Ψ -> (Γ Sin[π #] &);

Substituting Ψsol into eq91c, collecting the coefficient of Sin[π x], solving the resulting equation for Γ, and simplifying it using frequencies and values1, we obtain

ΓRule = 
      Solve[
                Coefficient[
                    eq91c[[1]] /. Ψsol /. 
                      int -> Integrate, Sin[π x]] == 0, Γ]
                [[1]] /. Ω -> 2  (Subscript[ω, 2]) /. 
          frequencies /. values1

Because the buckled-beam problem is self-adjoint, we have

adjoint = {{-I (Subscript[ω, 1]) (Subscript[Φ, 1])[x], 
          (Subscript[Φ, 1])[x]}, {-I (Subscript[ω, 2]) (Subscript[Φ, 2])[x], 
          (Subscript[Φ, 2])[x]}};

whose complex conjugate is

adjointC = adjoint /. conjugateRule

9.1.2 Second-Order Solution

Substituting sol1, Ψsol, and some of the parameter values into eqEps[[2]] yields

\(order2Eq = 
      \(\((linearSys /. u_\_1 -> u\_2)\) == 
              \((\(\(\(\(\(\(\((linearSys /. u_\_1 -> u\_2)\) - 
                                  \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([2]\)])
                                    \) /. sol1\) /. Ψsol\) /. 
                            ΓRule\) /. intRule1\) //. intRule2
                        \) /. modeshapes\) /. int -> Integrate // Expand)
                \) /. bucklingLoad\) /. values1 // Thread; 

    order2Eq /. displayRule\)

In order to collect the terms that may lead to secular terms from the right-hand sides of order2Eq, we define the rules

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

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

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

\(ST11 = 
      Coefficient[\(#[\([2]\)]\ &\)\ /@\ order2Eq /. expRule1[1], 
        E\^\(I\ ω\_1\ T\_0\)]; 
ST11 /. displayRule\)

\(ST12 = 
      Coefficient[\(#[\([2]\)]\ &\)\ /@\ order2Eq /. expRule1[2], 
        E\^\(I\ ω\_2\ T\_0\)]; 
ST12 /. displayRule\)

Demanding that ST11 be orthogonal to the adjointC[[1]], we obtain the solvability condition

SCond11=Solve[int[adjointC[[1]] . ST11, {x, 0, 1}] == 0 /. intRule1
                       //. intRule2 /. modeshapes /. 
                int -> Integrate,(Subscript[A, 1]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll

Demanding that ST12 be orthogonal to the adjointC[[2]], we obtain the solvability condition

SCond12=Solve[int[adjointC[[2]] . ST12, {x, 0, 1}] == 0 /. intRule1
                       //. intRule2 /. modeshapes /. 
                int -> Integrate,(Subscript[A, 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]//ExpandAll

Hence, we have the solvability conditions

SCond1 = Join[SCond11, SCond12];

whose complex conjugates are

ccSCond1 = SCond1 /. conjugateRule;

To remove (Subscript[D, 1])(Subscript[A, 1]) and (Subscript[D, 1])(Subscript[A, 2]) from the right-hand sides of order2Eq, we first define the rules

sigRule = 
      Solve[ResonanceConds, {(Subscript[σ, 1]), (Subscript[σ, 2])}][[1]]

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

Substituting SCond1, ccSCond1, and expRule2 into order2Eq, we obtain

\(order2Eqm = 
      \(\(order2Eq /. SCond1\) /. ccSCond1\) /. expRule2 // ExpandAll; 

    order2Eqm /. displayRule\)

We use the method of undetermined coefficients to determine the solution of order2Eqm and their associated boundary conditions. To implement this, we first determine all possible forms of the terms on the right-hand sides of order2Eqm as follows:

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

where all the spatial dependent functions satisfy the hinged-hinged boundary conditions. We note that sol2Form consists of (a) the terms that are proportional to (E^(-I (Subscript[ω, 1]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 1]) (Subscript[T, 0]))), (E^(-I (Subscript[ω, 2]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 2]) (Subscript[T, 0]) ))and (b) the remainder of the terms. The positions of the former terms, which are orthogonal to the adjoint, in sol2Form are

pos[1] = 
      (Position[
                  # sol2Form /. expRule1[1] /. Exp[_.  + _ (Subscript[T, 0])] -> 0, 
                  a_ /; a =!= 0, 1] // Flatten // Rest) & /@
         {(E^(-I (Subscript[ω, 1]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 1]) (Subscript[T, 0])))}

pos[2] = 
      (Position[
                  # sol2Form /. expRule1[2] /. Exp[_.  + _ (Subscript[T, 0])] -> 0, 
                  a_ /; a =!= 0, 1] // Flatten // Rest) & /@
         {(E^(-I (Subscript[ω, 2]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 2]) (Subscript[T, 0])))}

Next, we seek the solution of  (Subscript[u, 2]) and (Subscript[v, 2]) in two parts. First, we consider the part of sol2Form that is not related to the secular terms

sol2Forma = 
      Delete[sol2Form, {#} & /@ Flatten[Array[pos, 2]]]

and define the undetermined coefficients as

\(uSymbola = Table[ψ\_i, {i, Length[sol2Forma]}]; 

    vSymbola = uSymbola /. ψ -> η; \)

The general solution in terms of the sol2Forma can be written as

sol2a = {
        (Subscript[u, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            uSymbola . sol2Forma // Evaluate], 
        (Subscript[v, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            vSymbola . sol2Forma // Evaluate]};

Second, we consider the part of sol2Form that is related to the secular terms

sol2Formb1 = sol2Form[[pos[1][[1]]]]

sol2Formb2 = sol2Form[[pos[2][[1]]]]

and define the undetermined coefficients as

\(uSymbolb1 = Table[ψs\_\(1, i\), {i, Length[sol2Formb1]}]; 

    uSymbolb2 = Table[ψs\_\(2, i\), {i, Length[sol2Formb2]}]; \)

where the undetermined coefficients for (Subscript[v, 2]) are related to uSymbolb1 and uSymbolb2 due to the orthogonality condition. Therefore, the general solution in terms of the sol2Formb1 and sol2Formb2 can be written as

sol2b = {
        (Subscript[u, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            uSymbolb1 . sol2Formb1 + uSymbolb2 . sol2Formb2 // Evaluate], 
        (Subscript[v, 2]) -> Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            -I (Subscript[ω, 1]) uSymbolb1 . sol2Formb1 - 
                I (Subscript[ω, 2]) uSymbolb2 . sol2Formb2 // Evaluate]};

The total solution is the sum of sol2a, sol2b, and the complex conjugate of sol2b.

Next, we solve for these undetermined coefficients. Substituting sol2a into order2Eqm and equating the coefficients of sol2Forma, we have

algEqa = 
      Flatten[Coefficient[
                    Subtract @@ # /. sol2a /. intRule1 //. 
                          intRule2 /. modeshapes /. int -> Integrate, 
                    sol2Forma] & /@ order2Eqm] == 0 /. 
          Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Thread;

Substituting some of the parameter values into algEqa and solving for uSymbola and vSymbola, we obtain

symbolaRule = 
      Solve[algEqa, uSymbola~Join~vSymbola][[1]] /. 
            Ω -> 2  (Subscript[ω, 2]) /. frequencies /. values1

Substituting sol2b into either of the two equations in order2Eqm, equating the coefficients of sol2Formb1 and sol2Formb2, and solving for uSymbolb1 and uSymbolb2, we have

symbolbRule = 
      Solve[
                Coefficient[Subtract @@ order2Eqm[[1]] /. sol2b, 
                        sol2Formb1~Join~sol2Formb2] == 0 /. 
                    Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Thread, 
                uSymbolb1~Join~uSymbolb2][[1]] /. 
            Ω -> 2  (Subscript[ω, 2]) /. frequencies /. values1

Substituting these symbol values into sol2a and sol2b yields

\(sol2aForm = 
      \({u\_2[x, T\_0, T\_1, T\_2], v\_2[x, T\_0, T\_1, T\_2]} /. sol2a\) /. 
        symbolaRule; 

    sol2bForm = 
      \({u\_2[x, T\_0, T\_1, T\_2], v\_2[x, T\_0, T\_1, T\_2]} /. sol2b\) /. 
        symbolbRule; \)

sol2 = {
        (Subscript[u, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2aForm[[1]] + sol2bForm[[1]] + 
                  (sol2bForm[[1]] /. conjugateRule) // Expand // 
              Evaluate], 
        (Subscript[v, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2aForm[[2]] + sol2bForm[[2]] + 
                  (sol2bForm[[2]] /. conjugateRule) // Expand // 
              Evaluate]}

9.1.3 Solvability Conditions

Substituting sol1, Ψsol, sol2, and some of the parameter values into eqEps[[3]] yields

\(\(order3Eq = 
        \((linearSys /. u_\_1 -> u\_3)\) == 
              \((\(\(\(\(\(\((linearSys /. u_\_1 -> u\_3)\) - 
                                \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([3]\)])
                                  \) /. sol1\) /. Ψsol\) /. 
                          ΓRule\) /. sol2\) /. intRule1\) //. 
                    intRule2 // Expand)\) /. values1 // Thread; \) // Timing
      \)

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

ST21 = 
      Coefficient[#[[2]] & /@ order3Eq /. expRule1[1], 
            (E^(I (Subscript[ω, 1]) (Subscript[T, 0])))] /. modeshapes /. int -> Integrate;

ST22 = 
      Coefficient[#[[2]] & /@ order3Eq /. expRule1[2], 
            (E^(I (Subscript[ω, 2]) (Subscript[T, 0])))] /. modeshapes /. int -> Integrate;

Demanding that ST21 be orthogonal to the adjointC[[1]], we obtain the solvability condition

SCond21=Solve[(int[adjointC[[1]] . ST21, {x, 0, 1}] /. 
                                  intRule1 //. intRule2) == 0,(Subscript[A, 1]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]]/.modeshapes/.int -> Integrate/.SCond1/.frequencies/.values1//ExpandAll;
SCond21 /. displayRule

Demanding that ST22 be orthogonal to the adjointC[[2]], we obtain the solvability condition

SCond22=Solve[(int[adjointC[[2]] . ST22, {x, 0, 1}] /. 
                                    intRule1 //. intRule2) == 0,(Subscript[A, 2]^((0, 1)))[(Subscript[T, 1]), 
                              (Subscript[T, 2])]][[1]]/.modeshapes/.int -> Integrate/.SCond1/.Ω -> 2  (Subscript[ω, 2])/.frequencies/.values1//ExpandAll;
SCond22 /. displayRule

Because in the absence of damping the system is conservative, the modulation equations, solvability conditions, must satisfy symmetry conditions. To check for these symmetries, we let

\(\(formList1a = {A\&_\_1[T\_1, T\_2]\ A\&_\_2[T\_1, T\_2], 
        A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_2[T\_1, T\_2]}; \)\)

a1coefs = {(Subscript[s, 11]), (Subscript[s, 12])} -> 
          Coefficient[(2 (Subscript[ω, 1]) I) SCond21[[1, 2]], 
            formList1a] // Thread

\(\(formList2a = {
        A\_1[T\_1, T\_2]\ A\_2[T\_1, T\_2]\ A\&_\_1[T\_1, T\_2], 
        A\&_\_1[T\_1, T\_2]\^2}; \)\)

a2coefs = {(Subscript[s, 21]), (Subscript[s, 22])} -> 
          Coefficient[(2 (Subscript[ω, 2]) I) SCond22[[1, 2]], 
            formList2a] // Thread

The symmetry condition demands that (Subscript[s, 12]) = (Subscript[s, 21]), which is true because

(Subscript[s, 12]) - (Subscript[s, 21]) /. a1coefs /. a2coefs /. frequencies /. 
      values1

9.2 Four-Mode Interactions in Suspended Cables

We determine a second-order uniform asymptotic expansion of the three-dimensional response of a suspended cable to a transverse harmonic distributed excitation of one of the first two vertical or horizontal modes when their frequencies are in the ratio of either two-to-one or one-to-one and no other internal resonance is activated. Benedettini, Rega, and Alaggio (1995) derived nondimensional equations governing three-dimensional finite-amplitude vibrations of a suspended homogeneous elastic cable excited by harmonic distributed vertical and horizontal loads. The time derivatives are expressed in second-order form. As discussed in Chapter 5, second- and higher-order treatments of such equations lead to inconsistent results unless the time derivatives are expressed in first-order rather than second-order form. Therefore, we rewrite these equations as

(Subscript[δ, i_Integer, j_Integer]) := If[i == j, 1, 0]

\(eq1 = Table[u\_\(j, t\) - v\_j == 0, {j, 2}]; 

    eq2 = Table[
        v\_\(j, t\) + 2\ μ\_j\ v\_j - u\_\(j, x, x\) - 
            α\ 
              \((b\ \(y\^″\)[x]\ δ\_\(1, j\) + 
                  u\_\(j, x, x\))\)\ 
              \(∫\_0\%1
                  \((b\ \(y\^′\)[x]\ u\_\(1, x\) + 
                      1\/2\ \((u\_\(1, x\)\%2 + u\_\(2, x\)\%2)\))\) 
                  \[DifferentialD]x\) == 
          P\_j[x]\ Cos[Ω\ t + τ\_j], {j, 2}]; \)

\(EOM = 
      Transpose[{eq1, eq2}] /. {u_\_\(s_, m__\) :> ∂\_m u\_s[x, t], 
            v\_s_ :> v\_s[x, t]} // Flatten\)

yRule = {y -> (4 # (1 - #) &)};

Moreover, we write the associated boundary conditions as

BC = {(Subscript[u, j])[0, t] == 0, (Subscript[u, j])[1, t] == 0};

Here, the index j = 1 refers to the vertical component of the displacement and j = 2 refers to the horizontal component of the displacement; the (Subscript[μ, j]) are the viscous damping coefficients; the (Subscript[P, j])[x] are the distributed amplitudes of the excitation with the frequency Ω;  (Subscript[δ, i,  j]) denotes the Kronecker delta; and b y[x], where y[x] = 4 x (1 - x), defined in the yRule, is the initial static configuration. The spatial and temporal variables were nondimensionalized using the span of the cable ℓ and the characteristic time ℓ (Sqrt[ ϱ/E]), where ϱ and E are the cable density and modulous of elasticity, respectively.

We seek an asymptotic expansion of the response of the cable when its span is such that the first vertical and horizontal natural frequencies are in the ratio of two-to-one; that is, (b^2) α≈(π^2)/16. When (b^2) α=(π^2)/16, the natural frequencies are

params = {(Subscript[ω, j_ /; j != 3]) :> 2 π, 
        (Subscript[ω, 3]) -> π, (b^2) α -> (π^2)/16};

and the corresponding normalized modeshapes are

modeshapes = {
        (Subscript[Φ, 1]) -> 
          ((Sqrt[(2/3)]) (1 - Cos[2  π #]) &), 
        (Subscript[Φ, 2]) -> ((Sqrt[2]) Sin[2  π #] &), 
        (Subscript[Φ, 3]) -> ((Sqrt[2]) Sin[π #] &), 
        (Subscript[Φ, 4]) -> ((Sqrt[2]) Sin[2  π #] &)};

We use the method of multiple scales and seek a second-order uniform asymptotic expansion of the response of the cable to a primary-resonance excitation of the first vertical mode in the form

\(\(multiScales = {u_\_j_[x, t] -> u\_j[x, T\_0, T\_1, T\_2], 
        \(\(Derivative[m_, n_]\)[u_]\)[x, t] -> 
          \(dt[n]\)[D[u[x, T\_0, T\_1, T\_2], {x, m}]], t -> T\_0}; \)\)

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

Taking into account the resonance conditions, we scale the damping and forcing terms as

scaling = {(Subscript[μ, j])_ -> (ϵ^2) (Subscript[μ, j]), 
        (Subscript[P, j])_[x] -> (ϵ^3) (Subscript[P, j])[x]};

where the forcing has been scaled at order (ϵ^3) so that its influence first appears at the same order as the nonlinear shift in the frequencies in the absence of the internal resonance.

We consider the case in which the orders of magnitude of the forcing levels for both the vertical and horizontal motions are the same. Otherwise, a primary resonance of the vertical modes (current study) could also initiate a subharmonic resonance of order one-half  of the first horizontal mode and a primary resonance of the second horizontal mode at different levels of approximation.

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

eq92a = 
      ((EOM /. Integrate -> int /. multiScales /. solRule /. 
                    scaling // ExpandAll) //. intRule2 // ExpandAll) /. 
        (ϵ^(n_ /; n > 3)) -> 0;

Equating coefficients of like powers of ϵ in eq92a, we obtain

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

9.2.1 First-Order Solution

The homogeneous parts of the first-order equations are

(linearSys = #[[1]] & /@ eqEps[[1]]) /. displayRule

The first-order problems are identical to the linear eigenvalue problems. To account for the interactions arising from the multiple internal resonances and the primary resonance, we include the lowest two modes in each plane;  the other modes will decay due to damping. Hence, the solution of eqEps[[1]] can be expressed as

\(\(sol1u = {
        u\_\(1, 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, 1, 2}] // Evaluate], 
        u\_\(2, 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, 3, 4}] // Evaluate]}; 
    \)\)

sol1v = 
      Table[(Subscript[v, k, 1]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            D[(Subscript[u, k, 1])[x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])], (Subscript[T, 0])] /. sol1u // Evaluate], {
          k, 2}];

sol1 = Join[sol1u, sol1v]

where (Subscript[Φ, 1])[x] and (Subscript[Φ, 2])[x] are the lowest symmetric and antisymmetric vertical eigenmodes, and (Subscript[Φ, 3])[x] and (Subscript[Φ, 4])[x] are the lowest symmetric and antisymmetric horizontal eigenmodes. The (Subscript[Φ, n])[x] are orthogonal; they are normalized so that they satisfy the orthonormality condition \(TraditionalForm
      \`∫\_0\%1 Φ\_n[x]\ Φ\_m[x]\ 
            dx = δ\_\(n, m\)\).

Because the problem is self-adjoint, we have

adjoint = 
      Table[{-I (Subscript[ω, k]) (Subscript[Φ, k])[x], (Subscript[Φ, k])[x]}, {
          k, 4}];

whose complex conjugate is

adjointC = adjoint /. conjugateRule

9.2.2 Second-Order Solution

Substituting sol1 into eqEps[[2]] yields

\(order2Eq = 
      \((linearSys /. u_\_\(k_, 1\) -> u\_\(k, 2\))\) == 
          \((\(\(\(\(\(\((linearSys /. u_\_\(k_, 1\) -> u\_\(k, 2\))\) - 
                            \((\(#[\([1]\)]\ &\)\ /@\ eqEps[\([2]\)])\) /. 
                          sol1\) /. intRule1\) //. intRule2\) /. modeshapes
                    \) /. yRule\) /. int -> Integrate // Expand)\) // Thread; 
    
order2Eq /. displayRule\)

We consider the case of  primary resonance of the first vertical mode, a two-to-one internal resonance between the first vertical and horizontal modes, a one-to-one internal resonance between the first and second vertical modes, and a one-to-one internal resonance between the first vertical and  second horizontal modes. In order to collect the terms that may lead to secular terms from the right-hand sides of order2Eq, we define the rules

\(omgList = Table[ω\_i, {i, 4}]; 

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

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

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

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

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

    ST1v /. displayRule\)

Collecting the terms that may lead to secular terms from the equations governing the horizontal modes, we have

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

    ST1h /. displayRule\)

Demanding that ST1v[[i]] be orthogonal to the adjointC[[i]], we obtain the solvability conditions

SCond1v=Solve[(int[adjointC[[#]] . ST1v[[#]], {x, 0, 
                                      1}] /. intRule1 //. intRule2 /. 
                              modeshapes /. params /. int -> Integrate)
                         == 0,(Subscript[A, #]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]] & /@ {1, 2}//Flatten;
SCond1v /. displayRule

Demanding that ST1h[[i]] be orthogonal to the adjointC[[i+2]], we obtain the solvability conditions

SCond1h=Solve[(int[adjointC[[# + 2]] . ST1h[[#]], {x, 
                                      0, 1}] /. intRule1 //. intRule2 /. 
                              modeshapes /. params /. int -> Integrate)
                         == 0,(Subscript[A, # + 2]^((1, 0)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]] & /@ {1, 2}//Flatten;
SCond1h /. displayRule

SCond1 = Join[SCond1v, SCond1h];

whose complex conjugates are

ccSCond1 = SCond1 /. conjugateRule;

Next, we use the solvability conditions to eliminate the (Subscript[D, 1]) (Subscript[A, i ])from the right-hand sides of order2Eq. To this end, we define the rules

sigRule = 
      Solve[ResonanceConds, Table[(Subscript[σ, i]), {i, 4}]][[1]]

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

Substituting SCond1, ccSCond1, and expRule2 into order2Eq, we obtain

order2Eqm = order2Eq /. SCond1 /. ccSCond1 /. expRule2;

To determine the solution of order2Eqm and the associated boundary conditions, we use the method of undetermined coefficients. To accomplish this, we first determine the forms of the terms on the right-hand sides of order2Eqm by picking up the x dependent or x and (Subscript[T, 0]) dependent terms by defining the rule

baits = {_ (E^(a_ (Subscript[T, 0]) + b_. )) f_[c_ x] :> 
          (E^(a (Subscript[T, 0]) + b)) f[c x], 
        a_ f_[b_ x] :> f[b x] /; FreeQ[a, (Subscript[T, 0])]};

Then, we collect the solution forms for the horizontal modes according to

sol2Formh = 
      Table[Cases[order2Eqm[[k, 2]], #] & /@ baits, {k, 3, 
              4}] // Flatten // Union

Each of these terms satisfies the boundary conditions. However, not all of the solution forms

sol2Formv = 
      Table[Cases[order2Eqm[[k, 2]], #] & /@ baits, {k, 2}] // 
          Flatten // Union

for the vertical modes satisfy the boundary conditions. Besides, there exist terms that are not functions of x on the right-hand sides of order2Eqm[[{1,2}]], namely

#[[2]] & /@ order2Eqm[[{1, 2}]] /. 
        Thread[sol2Formv -> 0] /. displayRule

Hence, we need to either consider some algebraic forms or combined algebraic and sinusoidal forms for x to manually adjust the forms in sol2Formv or use a more general approach to let Mathematica pick up the right forms. The latter approach is used here.

Vertical Modes

Instead of looking for the x dependent forms, we seek all of the (Subscript[T, 0]) dependent forms as

sol2Formv1 = 
      Join[{1}, 
        Table[Cases[order2Eqm[[k, 2]], Exp[_], ∞], {k, 2}] // 
            Flatten // Union]

where the first argument in sol2Formv1 is included to account for the (Subscript[T, 0]) independent terms.

We define the undetermined coefficients for (Subscript[u, 1,  2]) and (Subscript[v, 1,  2]) as

\(uSymbolv = Table[ψ\_\(1, i\)[x], {i, Length[sol2Formv1]}]; 

    vSymbolv = uSymbolv /. ψ -> η; \)

and write the general solution in the form

sol2v = {
        (Subscript[u, 1, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, uSymbolv . sol2Formv1 // Evaluate], 
        (Subscript[v, 1, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            vSymbolv . sol2Formv1 // Evaluate]};

Substituting sol2v into order2Eqm[[1]] and equating the coefficients of the sol2Formv1 on both sides yields

expr1 = Subtract @@ order2Eqm[[1]] /. sol2v;

eqv1 = 
      Join[{expr1 == 0 /. Exp[_] -> 0}, 
        Coefficient[expr1, Rest[sol2Formv1]] == 0 // Thread]

Solving for vSymbolv from eqv1, we obtain

eta1Sol = Solve[eqv1, vSymbolv][[1]]

Substituting sol2v and eta1Sol into order2Eqm[[2]] and equating the coefficients of the sol2Formv1 on both sides yields

expr2 = 
      Subtract @@ order2Eqm[[2]] /. sol2v /. eta1Sol /. 
              intRule1 //. intRule2 /. 
          int[(y^′)[x] ((Subscript[ψ, 1, k_])^′)[x], _] -> 
            (Subscript[β, k]) /. yRule;

eqv2 = 
      Join[{expr2 == 0 /. Exp[_] -> 0}, 
          Coefficient[expr2, Rest[sol2Formv1]] == 0 // Thread] /. params

Next, we determine the solutions of the boundary-value problems: eqv2 and the associated boundary conditions. First, we obtain the information about the solution forms uSymbolv according to

\(ψForm = 
      \(\(\(DSolve[{#[\([1]\)], #[\([2]\)] == 0 /. x -> 0, #[\([2]\)] == 0 /. 
                        x -> 1}, #[\([2]\)], x]\)[\([1]\)] /. C[i_] -> 0\ &
              \)\ /@\ Transpose[{eqv2, uSymbolv}] // Flatten\) // ExpandAll; 
    
ψForm /. displayRule\)

Then, according to ψForm, we pick up the spatial function forms for uSymbolv as

bait1 = {a_.  (x^n)_.  :> (x^n) /; FreeQ[a, x], 
        a_ f_[b_ x] :> f[b x] /; FreeQ[a, x], 
        a_ x f_[b_ x] :> x f[b x]};

funcList = 
      Table[Prepend[
          Cases[ψForm[[i, 2]], #] & /@ bait1 // Flatten
               // Union, 1], {i, Length[ψForm]}]

We assume the solution for (Subscript[ψ, 1,  i]) as a linear combination of funcList[[i]], substitute it into ψForm[[i]],  equate both sides, solve for the undetermined coefficients, and obtain the solution for (Subscript[ψ, 1,  i]) as

\(psi1Sol = 
      Table[
		form = funcList[\([k]\)]; 
		
        cList = Table[c\_i, {i, Length[form]}]; 
		sol = cList . form; 
        algExpr = 
          sol - ψForm[\([k, 2]\)] /. 
            β\_k_ -> 
              Integrate[
                Expand[\(y\^′\)[x]\ D[sol, x] /. yRule], {x, 0, 1}]; 

        		algEq = 
          Append[Coefficient[algExpr, Rest[form]], 
                algExpr /. Thread[Rest[form] -> 0]] == 0 // Thread; 
	\ \ 
        ψ\_\(1, k\) -> 
          Function[x, 
            \(sol /. \(Solve[algEq, cList]\)[\([1]\)] // Expand\) // 
              Evaluate], 
	\ {k, Length[ψForm]}]\)

We note that all of the spatial dependent functions satisfy the boundary conditions.

The sol2Formv1 consists of two parts: a part is related to the secular terms, which consists of terms proportional to (E^(±I (Subscript[ω, 1]) (Subscript[T, 0]))) and (E^(±I (Subscript[ω, 2]) (Subscript[T, 0]))) and the other part consists of terms that are not related to the secular terms. The positions of the former are

pos[1] = 
      (Position[
                  # sol2Formv1 /. expRule1[1] /. 
                    Exp[_.  + _ (Subscript[T, 0])] -> 0, a_ /; a =!= 0, 1] // Flatten
                 // Rest) & /@
         {(E^(-I (Subscript[ω, 1]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 1]) (Subscript[T, 0])))}

pos[2] = 
      (Position[
                  # sol2Formv1 /. expRule1[2] /. 
                    Exp[_.  + _ (Subscript[T, 0])] -> 0, a_ /; a =!= 0, 1] // Flatten
                 // Rest) & /@
         {(E^(-I (Subscript[ω, 2]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 2]) (Subscript[T, 0])))}

Since pos[1] and pos[2] are the same, we simply define

{pos1, pos1cc} = pos[1] // Flatten

Next, we augment (Subscript[ψ, 1, pos1]) and (Subscript[η, 1, pos1]) by a linear combination of all possible homogeneous solutions. The unknown coefficients can be determined by the orthogonality conditions. ((Subscript[ψ, 1, pos1cc]) and (Subscript[η, 1, pos1cc]) are simply the complex conjugates of (Subscript[ψ, 1, pos1]) and (Subscript[η, 1, pos1]), respectively.) First, we express their solution forms using psi1Sol and eta1Sol as

STFormv = 
      {(Subscript[ψ, 1, pos1])[x] + 
                      Sum[(Subscript[c, 1, j]) (Subscript[Φ, j])[x], {j, 2}], 
                    (Subscript[η, 1, pos1])[x] + 
                      Sum[(Subscript[c, 1, j]) I (Subscript[ω, j]) (Subscript[Φ, j])[x], {
                          j, 2}]} /. eta1Sol /. psi1Sol /. 
              (Subscript[c, i])_ -> 0 /. modeshapes /. params // Expand

Demanding that STFormv be orthogonal to the adjointC[[i]], we obtain the unknown coefficients

c1Rule = 
      Solve[int[adjointC[[#]] . STFormv, {x, 0, 1}] == 0 /. 
                        modeshapes /. params /. intRule1 //. intRule2
                   /. int -> Integrate, (Subscript[c, 1, #])] & /@ {1, 2} // 
        Flatten

Then, we adjust psi1Sol with

psi1Sol[[pos1]] = 
      (Subscript[ψ, 1, pos1]) -> 
        Function[x, STFormv[[1]] /. c1Rule // Evaluate]

psi1Sol[[pos1cc]] = 
      (Subscript[ψ, 1, pos1cc]) -> 
        Function[x, 
          STFormv[[1]] /. c1Rule /. conjugateRule // Evaluate]

Horizontal Modes

The sol2Formh consists of two parts: a part is related to the secular terms, which consists of terms proportional to (E^(±I (Subscript[ω, 3]) (Subscript[T, 0]))) and (E^(±I (Subscript[ω, 4]) (Subscript[T, 0]))) and the other part consists of terms that are not related to the secular terms. The positions of the former are

pos[3] = 
      (Position[
                  # sol2Formh /. expRule1[3] /. Exp[_.  + _ (Subscript[T, 0])] -> 0, 
                  a_ /; a =!= 0, 1] // Flatten // Rest) & /@
         {(E^(-I (Subscript[ω, 3]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 3]) (Subscript[T, 0])))}

pos[4] = 
      (Position[
                  # sol2Formh /. expRule1[4] /. Exp[_.  + _ (Subscript[T, 0])] -> 0, 
                  a_ /; a =!= 0, 1] // Flatten // Rest) & /@
         {(E^(-I (Subscript[ω, 4]) (Subscript[T, 0]))), (E^(I (Subscript[ω, 4]) (Subscript[T, 0])))}

Hence, we assume that the solution for (Subscript[u, 2,  2]) and (Subscript[v, 2,  2]) consists of two parts corresponding to the parts in sol2Formh. The part, which is not related to the secular terms, is

sol2Formha = Delete[sol2Formh, pos[3]]

We define the corresponding undetermined coefficients as

\(uSymbolha = Table[ψ\_\(2, i\), {i, Length[sol2Formha]}]; 

    vSymbolha = uSymbolha /. ψ -> η; \)

Hence, the general solution in terms of the sol2Formha can be written as

sol2ha = {
        (Subscript[u, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            uSymbolha . sol2Formha // Evaluate], 
        (Subscript[v, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            vSymbolha . sol2Formha // Evaluate]};

The part of sol2Formh, which is related to the secular terms, is

sol2Formhb = sol2Formh[[pos[3][[1]]]]

and define the undetermined coefficients as

uSymbolhb = Table[ψ(Subscript[s, 2, i]), {i, Length[sol2Formhb]}];

where the undetermined coefficients for (Subscript[v, 2,  2]) are related to the uSymbolhb due to the orthogonality condition. Therefore, the general solution in terms of the sol2Formhb can be written as

sol2hb = {
        (Subscript[u, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            uSymbolhb . sol2Formhb // Evaluate], 
        (Subscript[v, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            -I (Subscript[ω, 3]) uSymbolhb . sol2Formhb // Evaluate]};

The total solution is a combination of sol2ha, sol2hb, and the complex conjugate of sol2hb.

Substituting sol2ha into order2Eqm[[{3,4}]] and equating the coefficients of sol2Formha on both sides, we have

algEqa = 
      Flatten[Coefficient[Subtract @@ # /. sol2ha, sol2Formha] & /@ 
              order2Eqm[[{3, 4}]]] == 0 // Thread;

Solving for the uSymbolha and the vSymbolha yields

symbolaRule = 
      Solve[algEqa, uSymbolha~Join~vSymbolha][[1]] /. params

Substituting sol2hb into either of order2Eqm[[{3,4}]], equating the coefficients of sol2Formhb on both sides, and solving for the uSymbolhb, we have

symbolbRule = 
      Solve[Coefficient[Subtract @@ order2Eqm[[3]] /. sol2hb, 
                  sol2Formhb] == 0 // Thread, uSymbolhb][[1]] /. params

Substituting symbolaRule and symbolbRule into sol2ha and sol2hb yields

\(sol2haForm = 
      \({u\_\(2, 2\)[x, T\_0, T\_1, T\_2], 
            v\_\(2, 2\)[x, T\_0, T\_1, T\_2]} /. sol2ha\) /. symbolaRule; 

    sol2hbForm = 
      \({u\_\(2, 2\)[x, T\_0, T\_1, T\_2], 
            v\_\(2, 2\)[x, T\_0, T\_1, T\_2]} /. sol2hb\) /. symbolbRule; 
    \)

sol2h = {
        (Subscript[u, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2haForm[[1]] + sol2hbForm[[1]] + 
                  (sol2hbForm[[1]] /. conjugateRule) // Expand // 
              Evaluate], 
        (Subscript[v, 2, 2]) -> 
          Function[{x, (Subscript[T, 0]), (Subscript[T, 1]), (Subscript[T, 2])}, 
            sol2haForm[[2]] + sol2hbForm[[2]] + 
                  (sol2hbForm[[2]] /. conjugateRule) // Expand // 
              Evaluate]}

9.2.3 Solvability Conditions

Substituting sol1, sol2v, and sol2h into eqEps[[3]] yields

\(\(order3Eq = 
        \((linearSys /. u_\_\(k_, 1\) -> u\_\(k, 3\))\) == 
            \((\(\(\(\(\(\(\(\(\(\((linearSys /. u_\_\(k_, 1\) -> u\_\(k, 3\))
                                        \) - 
                                      \((\(#[\([1]\)]\ &\)\ /@\ 
                                        eqEps[\([3]\)])\) /. sol1\) /. sol2v
                                  \) /. eta1Sol\) /. psi1Sol\) /. sol2h\) /. 
                          intRule1\) //. intRule2\) /. modeshapes\) /. yRule
                    \) /. int -> Integrate // Expand)\) // Thread; \) // 
      Timing\)

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

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

Collecting the terms that may lead to secular terms from the equations governing the horizontal modes, we have

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

Demanding that ST2v[[i]] be orthogonal to the adjointC[[i]], we obtain the solvability conditions

SCond2v=Solve[(int[adjointC[[#]] . ST2v[[#]], {x, 
                                        0, 1}] /. intRule1 //. intRule2
                                   /. modeshapes /. params /. 
                            int -> Integrate) == 0,(Subscript[A, #]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]] & /@ {1, 2}//Flatten//ExpandAll;
SCond2v /. displayRule

Demanding that ST2h[[i]] be orthogonal to the adjointC[[i+2]], we obtain the solvability conditions

SCond2h=Solve[(int[adjointC[[# + 2]] . ST2h[[#]], {
                                        x, 0, 1}] /. intRule1 //. intRule2
                                   /. modeshapes /. params /. 
                            int -> Integrate) == 0,(Subscript[A, # + 2]^((0, 1)))[(Subscript[T, 1]), (Subscript[T, 2])]][[1]] & /@ {1, 2}//Flatten//ExpandAll;
SCond2h /. displayRule

We can rewrite SCond2v and SCond2h in compact form by collecting the coefficients of the same form as follows:

\(\(collectForm = 
      Table[{A\_i[T\_1, T\_2]\ E\^\(I\ ω\_i\ T\_0\), 
            A\&_\_i[T\_1, T\_2]\ E\^\(\(-I\)\ ω\_i\ T\_0\)}, {i, 
            4}] // Flatten; \)\)

cubicTerms = 
      Nest[Outer[Times, collectForm, #] &, collectForm, 2] // Flatten // 
        Union;

cubicST[i_] := 
      (E^(-I (Subscript[ω, i]) (Subscript[T, 0]))) cubicTerms /. expRule1[i] /. 
            Exp[_ (Subscript[T, 0]) + _. ] -> 0 // Union // Rest

\(moduEq = 
      MapIndexed[
        \((term = cubicST[#2[\([1]\)]]; 
            coef = Coefficient[#1[\([2]\)], term]; #1[\([1]\)] == 
              \((#1[\([2]\)] /. Thread[term -> 0])\) + coef . term)\)\ &, 
        Join[SCond2v, SCond2h]]; 
moduEq /. displayRule\)

These solvability conditions satisfy the symmetry properties because

\(formList = 
      Table[A\_j[T\_1, T\_2]\ 
          Table[A\_i[T\_1, T\_2]\ A\&_\_i[T\_1, T\_2], {i, 4}], {j, 4}]; 

    formList /. displayRule\)

notationList = Table[(Subscript[s, i, j]), {i, 4}, {j, 4}]

\(\(coefRule = 
      MapThread[
          Rule, {notationList, 
            \(Coefficient[\((2\ ω\_#\ I)\)\ moduEq[\([#, 2]\)], 
                  formList[\([#]\)]]\ &\)\ /@\ Range[4]}, 2] // Flatten; 
    \)\)

testList = 
      Table[(Subscript[s, i, j]) - (Subscript[s, j, i]), {i, 4}, {j, i, 4}] // Flatten // 
          Union // Rest

testList /. coefRule /. params // Factor