Nonlinear chaotic responses of obliquely stiffened porous FG doubly curved shallow shells under primary and internal resonances

Abstract

This study investigates the nonlinear dynamic (ND) behavior of obliquely stiffened porous functionally graded (PFG) doubly curved shallow shells (DCSSs) exposed to external excitation, with particular emphasis on primary resonance (PR) and 1:2 internal resonance (IR). Both the shell and the oblique stiffeners (OSs) are made of PFG materials. Two PFG layouts are examined: (1) shells graded from ceramic-rich at the surfaces toward metal-rich, with uniform porosity through the thickness; and (2) the same grading is used, but the porosity varies through the thickness. In both layouts, the stiffeners are functionally graded (FG) members with uniform porosity through their thickness. A mechanical model is developed for PFG shells reinforced by two crossed families of OSs, where the stiffener angles can be identical or different. Based on first-order shear deformation theory (FSDT), von Kármán kinematics, and Hooke’s law, the stress–strain relationships are established, and the governing partial differential equations are obtained via Hamilton’s principle. These equations are then converted into a two-degree-of-freedom nonlinear ordinary differential system using the Galerkin method. After that, the method of multiple scales (MMSs) is used to derive a four-dimensional system of nonlinear averaged equations. The MMS-based analytical formulation is therefore established within a weakly nonlinear framework, in which small-amplitude responses, weak external excitation, and light damping are assumed. Finally, numerical simulations are conducted to explore key response features, such as phase portraits, Poincaré sections, and time histories, showing how material gradation patterns and stiffener orientation affect the shells’ nonlinear dynamics.

1. Introduction

Shallow shell structures are extensively adopted in engineering thanks to their lightweight nature and strong load-bearing capability. This advantageous feature makes them attractive for a broad range of applications, from aerospace components to large-span architectural domes. To enhance their mechanical performance even further, stiffening elements are often added to the shell structure. The effectiveness of these reinforced shells can be greatly improved by using advanced materials, especially FG materials (FGMs). Owing to their superior characteristics, FGMs have been widely adopted in many fields such as offshore structures, civil engineering, aerospace systems, marine and underwater vehicles, and petroleum facilities. As a result, applying FGMs in the design of stiffeners offers strong potential for improving the strength and stability of DCSSs. It should also be noted that porosity may form during the manufacturing process of these materials, which can adversely affect the structural integrity of FG-based systems.^1,2 This issue has attracted considerable research interest in recent years. Accordingly, numerous studies have been devoted to analyzing and improving the stability of such advanced structural systems.

FG shells and beams have attracted significant attention because they provide a smooth spatial variation of material properties, which can reduce interfacial stress concentrations and improve structural performance under severe mechanical, thermal, and dynamic loading conditions. In mechanical engineering applications, FG shell- and beam-type components are frequently encountered in aerospace panels, pressure vessels, marine and offshore structures, turbine components, biomedical implants, thermal protection systems, energy-related equipment, and lightweight structural members. In these systems, the gradual transition between ceramic-rich and metal-rich phases can be designed to improve stiffness, thermal resistance, vibration characteristics, and damage tolerance. For shell structures, curvature, transverse shear deformation, geometric nonlinearity, and possible thickness stretching may significantly influence the static and dynamic responses. Therefore, recent studies have examined the nonlinear static and dynamic behavior of FG shells by considering advanced deformation theories, porosity distributions, and thickness-stretching effects. In particular, the study on porosity effects on nonlinear static performance of FG shells considering thickness stretching has shown that porosity effects can considerably alter the structural response and should be carefully considered in the modeling of porous FG shell systems. ³ In parallel, FG beams have also been widely studied as fundamental structural elements in micro/nano-electromechanical systems, sensors, actuators, foundation-supported components, and lightweight mechanical structures. The vibration and stability of nonlinear nonlocal strain-gradient FG beams on visco-Pasternak foundations further demonstrate the importance of size effects, foundation interaction, and material gradation in predicting the response of advanced FG structural members. ⁴

A large body of research has been conducted on the vibration behavior of shells made from various types of FGMs, such as FG structures,^5–7 FG porous structures,^8–10 FG carbon nanotube-reinforced composites,^11–13 and PFG systems.^14–16 Among these studies, considerable attention has been given to DCSSs constructed from different FGMs.^17–19 Moreover, some investigations have addressed the vibration and resonance behavior of FG shells in a broader sense,^20–22 while others have focused specifically on FG DCSSs.^23–25 Recently, Foroutan et al. ²⁶ studied the ND response of laminated composite DCSSs under combined external and parametric excitations, taking both internal and primary resonances into account. From another perspective, porosity and size-dependent effects have been shown to play an important role in the mechanical response of FG structural components. For example, the nonlinear static performance of porous FG shells has been investigated by considering thickness-stretching effects, showing that both porosity distribution and transverse deformation assumptions can noticeably affect the predicted stiffness and load-deflection behavior. In addition, the vibration and stability characteristics of nonlinear nonlocal strain-gradient FG beams resting on visco-Pasternak foundations have been analyzed, demonstrating that nonlocality, strain-gradient effects, foundation parameters, and material gradation can significantly influence the dynamic stability and natural frequencies of FG beam-type structures. These studies provide important research foundations for the present work, because they confirm that accurate modeling of FG structures requires careful consideration of porosity, geometric nonlinearity, deformation assumptions, and interaction between material gradation and dynamic stability.

Stiffened structures are commonly employed in diverse engineering applications, yet their ND and vibration behavior have not been studied as extensively as their static or load-bearing performance. Much of the earlier work has mainly focused on improving structural efficiency through the use of stiffeners to increase load-carrying capacity.^27–30 Meanwhile, several researchers have investigated the resonance and vibration behavior of stiffened shell structures made of FGMs under different loading conditions.^31–33 Special emphasis has also been placed on stiffened FG-DCSSs, particularly with respect to their dynamic response.^34–36 More recently, Foroutan and Torabi ³⁷ examined the ND and combination resonance behavior of obliquely stiffened multilayer FG-DCSSs. In ³⁷ , only the shell was modeled as a PFG material, whereas the OSs are treated as FGMs, and only the combination resonance of obliquely stiffened multilayer FG-DCSSs with a single mode is investigated.

In addition to analytical and semi-analytical approaches, numerical methods have been widely used for the analysis of shell structures. Among them, the finite element (FE) method is generally regarded as one of the most powerful and flexible tools for modeling linear and nonlinear responses of complex shell systems. Recent developments have also introduced isogeometric finite element formulations, which provide a closer connection between geometric modeling and numerical analysis. For example, Milić et al. ³⁸ developed a Reissner–Mindlin-based isogeometric finite element formulation for piezoelectric active laminated shells, while Milić et al. ³⁹ further extended the isogeometric FE framework to geometrically nonlinear analysis of piezoelectric active laminated shells. In the context of stiffened shell dynamics, He and Mohammadian ⁴⁰ recently investigated high-accuracy nonlinear frequency estimation of stringer-stiffened shells using efficient analytical frequency-amplitude formulations. These studies demonstrate the continuing progress in advanced numerical and analytical modeling strategies for shell-type structures and provide useful context for positioning the present analytical/semi-analytical resonance formulation.

A review of the available literature shows that, although many works have investigated the ND and vibration behavior of obliquely stiffened FG-DCSSs, very limited attention has been given to such structures under external excitation, especially when PR and 1:2 IR occur simultaneously. Although the above-mentioned studies have significantly advanced the understanding of FG shells and beams, most of them have focused on static behavior, free vibration, foundation-supported beams, unstiffened shells, or single-resonance conditions. Recent studies on porosity effects on nonlinear static performances of functionally graded shells considering thickness stretching, as well as vibration and stability of nonlinear nonlocal strain-gradient FG beams on visco-Pasternak foundations, have further demonstrated that porosity distribution, deformation assumptions, size-dependent effects, foundation interaction, and material gradation can strongly influence the mechanical and dynamic responses of FG structural components. However, comparatively fewer studies have addressed the nonlinear chaotic response of obliquely stiffened porous FG doubly curved shallow shells under external excitation. Moreover, the simultaneous influence of porosity distribution, oblique stiffener orientation, material gradation in both the shell and stiffeners, and coupled primary/internal resonance has not been fully clarified. This motivates the present study, because such a coupled configuration combines several important mechanical features: the stiffness-to-weight efficiency of doubly curved shallow shells, the ability of oblique stiffeners to tailor stiffness distribution and modal coupling, the smooth variation of material properties provided by FG materials, and the influence of porosity on effective stiffness, inertia, and dynamic stability. Therefore, investigating this structure under simultaneous PR and IR conditions can provide a more realistic understanding of how geometry, reinforcement orientation, material gradation, and porosity jointly affect periodic, quasi-periodic, and chaotic dynamic responses. To the best of the authors’ knowledge, no previous ND investigation has provided a unified theoretical and numerical treatment of obliquely stiffened PFG-DCSSs that accounts simultaneously for PFG in both the DCSSs and the OSs, external excitations, and the coupled effects of PR and 1:2 IR. In addition, most earlier works on stiffened FG-DCSS structures have neither considered PFG in both the shell and the OSs nor systematically investigated how stiffener orientation and material gradation influence the transition between chaotic and periodic responses. A clear understanding of these complex vibration and resonance phenomena is important because of their relevance to advanced engineering fields such as aeronautics, petroleum systems, aerospace, and biomedical applications. Resonance and vibration-control concepts are also important in acoustic and flow-related engineering systems. For example, Zhao et al. ⁴¹ developed and experimentally validated a nonlinear model of a Helmholtz resonator with an oscillating volume, showing that active adjustment of the cavity volume can enhance acoustic damping and broaden the effective frequency range of the resonator. In light of these research gaps, the key contributions of the present work are listed as: (I) investigating the ND behavior of obliquely stiffened PFG-DCSSs under external excitation; (II) modeling both the DCSSs and the external OSs as PFG materials, with through-thickness material variation from metal-rich outer layers toward a ceramic-rich inner layer; (III) simultaneously addressing 1:2 IR and PR; and (IV) examining two PFG shell layouts with the same ceramic-rich-to-metal-rich gradation, considering uniform and nonuniform porosity distributions through the shell thickness, while FG stiffeners are assumed to have uniform porosity through their thickness.

The nonlinear governing equations (NGEs) of obliquely stiffened PFG-DCSSs are developed by combining the FSDT with von Kármán nonlinear kinematic equations. Next, the PFG-OSs are incorporated into the PFG-DCSSs model through the smeared-stiffener method, allowing the stiffener effects to be represented within the shell’s equivalent stiffness formulation. Galerkin’s method is then employed to discretize the NGEs of the system, and MMSs is used to build a unified theoretical approach for studying the ND behavior of the structure under PR and IR. Based on this framework, the study systematically examines the impact of material gradation and stiffener orientations on the ND behavior of the shells, using phase portraits, waveform plots, and Poincaré maps to illustrate the results. In this way, the findings provide useful insights and practical guidance for engineers and researchers concerned with the analysis and design of obliquely stiffened PFG-DCSSs subjected to external excitations.

It should be emphasized that the analytical resonance formulation developed in this study is based on the method of multiple scales and is therefore restricted to weakly nonlinear response regimes. This assumption is consistent with the present focus on PR and 1:2 IR under small-amplitude vibration, light damping, and weak external excitation. For stronger nonlinear regimes, alternative semi-analytical techniques, such as the homotopy perturbation method (HPM), may be employed as an extension of the present framework. In addition, recent developments in fractal and fractional nonlinear oscillators have shown that fractal derivatives, fractional derivatives, two-scale fractal theory, fractal complex transformation, and fractional complex transformation can be useful for modeling nonlinear oscillatory systems with complex spatial or temporal characteristics. For example, Zheng et al. ⁴² investigated a fractal oscillator for a pendulum with a rolling wheel using He’s fractal derivative and a two-scale fractal-based solution strategy, while Jin et al. ⁴³ analyzed a fractional rotating pendulum oscillator based on He’s fractional derivative and fractional variational theory. Inspired by these studies, future work may extend the present shell formulation into fractal or fractional space to capture more complex multi-scale or memory-dependent effects in porous FG shell systems.

2. Obliquely stiffened PFG-DCSSs

Figure 1 presents a schematic view of the obliquely stiffened PFG-DCSSs. The key geometric parameters of the DCSSs are the thickness $h$ , and the curved lengths $a$ and $b$ in the $x$ and $y$ directions, respectively, and $R_x$ and $R_y$ indicate the corresponding radii of curvature. The OSs are characterized by their width $d$ , spacing $S$ , thickness $h_s$ , and the stiffeners’ angles $\alpha$ and $\beta$ . The present formulation is developed for a general doubly curved shallow shell with two principal radii of curvature, $R_x$ and $R_y$ . Therefore, different shallow-shell geometries can be represented by selecting appropriate curvature parameters. In particular, $R_x=R_y$ corresponds to a spherical/synclastic shallow shell, whereas $R_x\ne R_y$ represents a general doubly curved shallow shell. The schematic shown in Figure 1 illustrates the spherical/synclastic shallow-shell configuration used in the main numerical examples. Figure 2 further illustrates how the material distribution varies within the obliquely stiffened PFG-DCSSs. As shown, two PFG layouts are considered. In the first layout, the shell is functionally graded from ceramic-rich surfaces toward a more metal-rich through the thickness, and the porosity is assumed uniform (i.e., constant) across the thickness. In the second layout, the shell follows the same through-thickness grading, but the porosity is nonuniform, meaning it varies along the thickness according to a prescribed distribution. In both cases, the OSs are FG members with uniform porosity through their thickness. In the present model, the common/contact layer between the PFG-OSs and the PFG-DCSSs is ceramic-rich and is assumed to be continuously attached. Therefore, the stiffeners and shell are considered to be perfectly bonded, and no relative displacement, sliding, or separation is allowed at the shell-stiffener interface. This assumption enables compatible deformation between the stiffeners and the shell and allows the stiffener contribution to be incorporated into the equivalent stiffness formulation of the shell. Delamination, imperfect bonding, and interfacial damage are not explicitly considered in the current formulation.

Figure 1.

Configuration of obliquely stiffened PFG-DCSSs for the spherical/synclastic shallow-shell case.

Figure 2.

Material distribution of obliquely stiffened PFG-DCSSs.

The mechanical properties of the DCSSs and OSs, namely Young’s modulus ( $E$ ) and mass density ( $\rho$ ), for the two configurations shown in Figure 2 are determined based on the corresponding material volume fractions (VFs).

• Type 1:

$P_{OSs}=P_c+\left(P_c-P_m\right){\left(\frac{2z+h}{2h_s}\right)}^{N_{OSs}}-\frac{1}{2}\gamma \left(P_m+P_c\right);-\left(h/2+h_s\right)\le z\le -h/2$

$P_{DCSSs}=P_c+\left(P_m-P_c\right){\left(\frac{2z+h}{2h}\right)}^{N_{DCSSs}}-\frac{1}{2}\gamma \left(P_m+P_c\right);-h/2\le z\le h/2$

(1a)

• • Type 2:

$P_{OSs}=P_c+\left(P_c-P_m\right){\left(\frac{2z+h}{2h_s}\right)}^{N_{OSs}}-\frac{1}{2}\gamma \left(P_m+P_c\right);-\left(h/2+h_s\right)\le z\le -h/2$

$P_{DCSSs}=P_c+\left(P_m-P_c\right){\left(\frac{2z+h}{2h}\right)}^{N_{DCSSs}}-\frac{1}{2}\gamma \left(P_m+P_c\right)\left(1-\frac{2\left|z\right|}{h}\right);-h/2\le z\le h/2$

(1b)

In Eq. (1), $P$ stands for a material property (such as $E$ or $\rho$ ). Also, $N_{DCSSs}$ and ${\gamma }_{DCSSs}$ , and $N_{OSs}$ and ${\gamma }_{OSs}$ , denote the material power-law indices and the porosity VFs of the DCSSs and OSs, respectively. The subscripts $c$ and $m$ indicate the ceramic and metal. It should be noted that Eq. (1) represents an equivalent-continuum quantification of the porosity distribution in the PFG-DCSSs and FG-OSs. In this approach, the effect of porosity is incorporated through prescribed through-thickness material-property functions, allowing the global influence of uniform and nonuniform porosity distributions on stiffness and inertia to be examined within a tractable shell-theory framework. The two porosity patterns considered in Eq. (1) are selected as representative cases to examine the influence of the assumed porosity distribution on the ND response of the shell. These porosity distributions have been widely adopted in previous studies on porous FG beams, plates, and shellsa^2,14,15; therefore, they are not arbitrary assumptions, but commonly used porosity idealizations within the equivalent-continuum modeling framework. Although two-scale fractal dimensions can provide a more detailed description of pore morphology, pore connectivity, and multi-scale pore-size distribution, such microstructural features are not explicitly resolved in the present formulation.

3. Governing equations

As noted earlier, the FSDT is used to derive the NGEs for investigating the ND behavior of obliquely stiffened PFG-DCSSs. The strain-displacement relationships, expressed in terms of the displacements and rotation components, are as:

$\left\{\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{l}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}+z\left\{\begin{array}{l}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right\};\hspace{0.17em}\left\{\begin{array}{l}{\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right\}=\left\{\begin{array}{l}{\partial }_x\\ {\partial }_y\end{array}\right\}w+\left\{\begin{array}{l}{\varphi }_x\\ {\varphi }_y\end{array}\right\};\hspace{0.17em}\left\{\begin{array}{l}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right\}=\left[\begin{array}{ll}{\partial }_x& 0\\ 0& {\partial }_y\\ {\partial }_y& {\partial }_x\end{array}\right]\left\{\begin{array}{l}{\varphi }_x\\ {\varphi }_y\end{array}\right\}$

(2)

Here, the in-plane normal strains $\left({\epsilon }_x^0,{\epsilon }_y^0\right)$ and the in-plane shear strain $\left({\gamma }_{xy}^0\right)$ are defined on the shell mid-surface. The terms ${\kappa }_x$ , ${\kappa }_y$ , and ${\kappa }_{xy}$ denote the changes in curvature and twist of the shell. In addition, ${\varphi }_x$ and ${\varphi }_y$ denote the rotations of the normal to the mid-surface about $y$ - and $x$ -axes, respectively. Accordingly, for a shallow shell, the mid-surface strain components based on the von Kármán nonlinear strain–displacement relations are given as^44,45:

$\left\{\begin{array}{l}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\}=\underset{\text{linear}}{\underbrace{\left[\begin{array}{lll}{\partial }_x& 0& -1/R_x\\ 0& {\partial }_y& -1/R_y\\ {\partial }_y& {\partial }_x& 0\end{array}\right]\left\{\begin{array}{l}u\\ v\\ w\end{array}\right\}}}+\underset{\text{nonlinear}}{\underbrace{\left\{\begin{array}{l}{\partial }_x^2/2\\ {\partial }_y^2/2\\ {\partial }_x{\partial }_y\end{array}\right\}w}}$

(3)

where, $u$ , $v$ , and $w$ represent the displacement components along the x-, y-, and z-axes, respectively. The stress-strain relations for the PFG-DCSSs and FG-OSs are expressed below:

• • PFG-DCSSs

${\left\{\begin{array}{l}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}}_{DCSSs}=\left[\begin{array}{lll}{S}_{11}& S_{12}& 0\\ S_{21}& S_{22}& 0\\ 0& 0& S_{66}\end{array}\right]\left\{\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}$

$S_{11}=S_{22}=\frac{E_{DCSSs}}{1-{\nu }^2};\hspace{0.17em}{S}_{12}=S_{21}=\frac{{\nu E}_{DCSSs}}{1-{\nu }^2};\hspace{0.17em}{S}_{66}=\frac{E_{DCSSs}}{2\left(1+\nu \right)}$

(4a)

• • FG-OSs

${\left\{\begin{array}{l}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}}_{OSs}=\left[\begin{array}{lll}{Q}_{11}& Q_{12}& 0\\ Q_{21}& Q_{22}& 0\\ 0& 0& Q_{66}\end{array}\right]\left\{\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}$

$Q_{11}=Q_T\frac{{\mathrm{C}\alpha }^3+{\mathrm{C}\theta }^3}{\mathrm{S}\alpha +\mathrm{S}\theta };\hspace{0.17em}{Q}_{12}=Q_T\frac{{\mathrm{S}\alpha }^2\mathrm{C}\alpha +{\mathrm{S}\theta }^2\mathrm{C}\theta }{\mathrm{C}\alpha +\mathrm{C}\theta };\hspace{0.17em}{Q}_{21}=Q_T\frac{\mathrm{S}\alpha {\mathrm{C}\alpha }^2+\mathrm{S}\theta {\mathrm{C}\theta }^2}{\mathrm{S}\alpha +\mathrm{S}\theta }$

$Q_{22}=Q_T\frac{{\mathrm{S}\alpha }^3+{\mathrm{S}\theta }^3}{\mathrm{C}\alpha +\mathrm{C}\theta };\hspace{0.17em}{Q}_{66}=Q_T\frac{\mathrm{S}\alpha \mathrm{C}\alpha +\mathrm{S}\theta \mathrm{C}\theta }{2};\hspace{0.17em}{Q}_T=\frac{d{E}_{OSs}\mathrm{S}\left(\alpha +\theta \right)}{S}$

(4b)

Here, $\nu$ denotes Poisson’s ratio. By applying Lekhnitskii’s smeared stiffener technique, we examine how the FG-OSs influence the PFG-DCSSs. The smeared stiffener technique represents the stiffeners through equivalent stiffness contributions distributed over the shell surface. This approach is suitable for regularly spaced stiffeners and for evaluating the global dynamic response of stiffened shell structures, including natural frequencies, modal coupling, and resonance behavior. However, it does not explicitly capture local effects associated with discrete stiffener-shell junctions, stress concentrations, local deformation around individual stiffeners, or possible local buckling of the stiffeners. Therefore, the present model should be interpreted as a global-response formulation rather than a detailed local stress analysis of the stiffener-shell connection regions. With this approach, and by integrating Eq. (3) through the thickness of the shell, we obtain the following governing expressions for resultant forces and moments:

$\left\{\begin{array}{l}N\\ M\end{array}\right\}=\left[\begin{array}{ll}A& B\\ B& D\end{array}\right]\left\{\begin{array}{l}{\epsilon }^0\\ \kappa \end{array}\right\};\hspace{0.17em}\left\{\begin{array}{l}{Q}_y\\ Q_x\end{array}\right\}=K\overline{A}\gamma$

(5)

where $K=5/6$ denotes the shear correction factor,^46,47 and

$N=\left\{\begin{array}{l}{N}_x\\ N_y\\ N_{xy}\end{array}\right\};\hspace{0.17em}M=\left\{\begin{array}{l}{M}_x\\ M_y\\ M_{xy}\end{array}\right\};\hspace{0.17em}{\epsilon }^0=\left\{\begin{array}{l}{\epsilon }_x^0\\ {\epsilon }_y^0\\ {\gamma }_{xy}^0\end{array}\right\};\hspace{0.17em}\kappa =\left\{\begin{array}{l}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right\};\hspace{0.17em}\mathrm{\Lambda }=\left[\begin{array}{lll}{\mathrm{\Lambda }}_{11}& {\mathrm{\Lambda }}_{12}& 0\\ {\mathrm{\Lambda }}_{21}& {\mathrm{\Lambda }}_{22}& 0\\ 0& 0& {\mathrm{\Lambda }}_{66}\end{array}\right];\hspace{0.17em}\mathrm{\Lambda }=A,B,D$

$\overline{A}=\left[\begin{array}{ll}{A}_{44}& 0\\ 0& A_{55}\end{array}\right];\hspace{0.17em}\gamma =\left\{\begin{array}{l}{\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}$

$\left(A_{ij},B_{ij},D_{ij}\hspace{0.17em}\right)={\int }_{-\left(h/2+h_s\right)}^{-h/2}{Q}_{ij}\left(1,z,z^2\right)dz+{\int }_{-h/2}^{h/2}{S}_{ij}\left(1,z,z^2\right)dz;\hspace{0.5em}i,j=1,2,6$

$A_{44}=A_{55}=A_{66}$

(6)

4. Nonlinear dynamic equations

Hamilton’s principle is employed to derive the governing partial differential equations for the obliquely stiffened PFG-DCSSs, as ⁴⁸ :

${\int }_0^T\left(\delta K-\delta U-\delta V\right)dt=0$

(7)

In Eq. (7), $\delta K$ , $\delta U$ , and $\delta V$ represent the virtual variations of the kinetic energy, strain energy, and the potential associated with the external excitation, respectively, and can be expressed as follows:

$\delta K={\int }_0^a\underset{-\left(h\backslash /2+h_s\right)}{\overset{h/2}{{\displaystyle \int }}}{\int }_0^b[\left(\rho \delta {\dot{\mathrm{d}}}^T\dot{\mathrm{d}}\right)\hspace{0.17em}dy\hspace{0.17em}dz\hspace{0.17em}dx+\rho z\left({\dot{\psi }}_x\delta \dot{u}+{\dot{\psi }}_y\delta \dot{v}+\dot{u}\delta {\dot{\psi }}_x+\dot{v}\delta {\dot{\psi }}_y\right)+\rho z^2\left({\dot{\psi }}_x\delta {\dot{\psi }}_x+{\dot{\psi }}_y\delta {\dot{\psi }}_y\right)]dy\hspace{0.17em}dz\hspace{0.17em}dx\hspace{0.17em}\delta U={\int }_0^a\underset{-\left(h/2+h_s\right)}{\overset{h/2}{{\displaystyle \int }}}{\int }_0^b\left(\delta {ϵ}^T\sigma \right)\hspace{0.17em}dy\hspace{0.17em}dz\hspace{0.17em}dx;\hspace{0.17em}\delta V={\int }_0^a\underset{-\left(h/2+h_s\right)}{\overset{h/2}{{\displaystyle \int }}}{\int }_0^b\left(\delta {\mathrm{d}}^T{q}^e\right)\hspace{0.17em}dy\hspace{0.17em}dz\hspace{0.17em}dx$

$\mathrm{d}={\left\{\begin{array}{lll}u& v& w\end{array}\right\}}^T;\hspace{0.17em}\sigma ={\left\{\begin{array}{lllll}{\sigma }_x& {\sigma }_y& {\tau }_{yz}& {\tau }_{xz}& {\tau }_{xy}\end{array}\right\}}^T;\hspace{0.17em}ϵ={\left\{\begin{array}{lllll}{\epsilon }_x& {\epsilon }_y& {\gamma }_{yz}& {\gamma }_{xz}& {\gamma }_{xy}\end{array}\right\}}^T$

(8)

where $\delta$ denotes the variation, and $q^e$ is the external load vector in three directions. Using Eq. (8), together with the fundamental lemma of variational calculus in terms of virtual displacements and integration by parts, the Lagrange equations of motion of the obliquely stiffened PFG-DCSSs are derived below:

$\delta u:{N_x}_{,x}+{N_{xy}}_{,y}={\rho }_1{u}_{,tt}+{\rho }_2{{\varphi }_x}_{,tt}$

(9a)

$\delta v:{N_{xy}}_{,x}+{N_y}_{,y}={\rho }_1{v}_{,tt}+{\rho }_2{{\varphi }_y}_{,tt}$

(9b)

$\delta w:{Q_x}_{,x}+{Q_y}_{,y}+N_x{w}_{,xx}+2N_{xy}{w}_{,xy}+N_y{w}_{,yy}+\frac{N_x}{R_x}+\frac{N_y}{R_y}+q={\rho }_1{w}_{,tt}+2{\rho }_{1}ct$

(9c)

$\delta {\varphi }_x:{M_x}_{,x}+{M_{xy}}_{,y}-Q_x={\rho }_2{u}_{,tt}+{\rho }_3{{\varphi }_x}_{,tt}$

(9d)

$\delta {\varphi }_y:{M_{xy}}_{,x}+{M_y}_{,y}-Q_y={\rho }_2{v}_{,tt}+{\rho }_3{{\varphi }_y}_{,tt}$

(9e)

Using Eq. (3), the compatibility equation for the obliquely stiffened PFG-DCSSs can be written as follows:

$\left[\begin{array}{lll}{\partial }_{yy}& {\partial }_{xx}& -{\partial }_{xy}\end{array}\right]{\epsilon }^0=\left({\partial }_{xy}^2-{\partial }_{xx}{\partial }_{yy}-\frac{{\partial }_{xx}}{R_y}-\frac{{\partial }_{yy}}{R_x}\right)w$

(10)

The stress function $\psi$ , which satisfies the required conditions, is introduced below:

$N^T=\left[\begin{array}{lll}{\partial }_{yy}& {\partial }_{xx}& -{\partial }_{xy}\end{array}\right]\psi$

(11)

Based on Eqs. (11), (9a), and (9b), it follows that:

$u_{,tt}=-\frac{{\rho }_2}{{\rho }_1}{{\varphi }_x}_{,tt};\hspace{0.5em}{v}_{,tt}=-\frac{{\rho }_2}{{\rho }_1}{{\varphi }_y}_{,tt}$

(12)

5. System’s solution

Assume that the obliquely stiffened PFG-DCSSs are simply supported. The simply supported boundary condition is adopted because it provides admissible trigonometric mode functions and allows the Galerkin reduction to be performed in an analytically tractable form. It should be noted that boundary conditions can affect the natural frequencies, mode shapes, modal coupling coefficients, and nonlinear resonance response of shell structures. Therefore, the results presented in this study are specific to the simply supported case. For other boundary conditions, the same general formulation may be extended by employing the corresponding admissible functions or numerical mode shapes in the Galerkin discretization. Based on these boundary conditions, the approximate solutions can be expressed as follows^46,49:

$w=\left[\begin{array}{ll}{W}_1\left(t\right)& W_2\left(t\right)\end{array}\right]\left[\begin{array}{l}\sin \frac{\pi x}{a}\sin \frac{3\pi y}{b}\\ \sin \frac{3\pi x}{a}\sin \frac{\pi y}{b}\end{array}\right];\hspace{0.17em}{\varphi }_x=\left[\begin{array}{ll}{\mathrm{\Phi }}_{x1}\left(t\right)& {\mathrm{\Phi }}_{x2}\left(t\right)\end{array}\right]\left[\begin{array}{l}\cos \frac{\pi x}{a}\sin \frac{3\pi y}{b}\\ \cos \frac{3\pi x}{a}\sin \frac{\pi y}{b}\end{array}\right]$

${\varphi }_y=\left[\begin{array}{ll}{\mathrm{\Phi }}_{y1}& {\mathrm{\Phi }}_{y2}\end{array}\right]\left[\begin{array}{l}\sin \frac{\pi x}{a}\cos \frac{3\pi y}{b}\\ \sin \frac{3\pi x}{a}\cos \frac{\pi y}{b}\end{array}\right];q=\left[\begin{array}{ll}{q}_1\left(t\right)& q_2\left(t\right)\end{array}\right]\left[\begin{array}{l}\sin \frac{\pi x}{a}\sin \frac{3\pi y}{b}\\ \sin \frac{3\pi x}{a}\sin \frac{\pi y}{b}\end{array}\right]$

(13)

The Galerkin method is employed in the present study because it provides an efficient and analytically tractable procedure for reducing the nonlinear governing partial differential equations of the obliquely stiffened PFG-DCSSs to a finite-dimensional system of ordinary differential equations. For the simply supported shell considered here, the selected trigonometric admissible functions satisfy the geometric boundary conditions and allow the governing equations to be projected onto the dominant vibration modes. This is especially suitable for the present resonance analysis, because the PR and 1:2 IR are mainly governed by the interaction of selected lower-frequency modes. In addition, the Galerkin procedure preserves the effects of material gradation, porosity distribution, shell curvature, inertia, and oblique stiffener orientation through the resulting modal coefficients. Therefore, the method provides a compact reduced-order model that is appropriate for the subsequent MMS-based nonlinear resonance analysis. In the present study, the lower-frequency vibration modes are considered to be more influential in the global ND response than higher-frequency modes. Therefore, two specific admissible modes are retained in Eq. (13) to represent the dominant interacting modal components associated with the considered PR and 1:2 IR conditions. This two-mode truncation allows the main coupled resonance mechanism to be captured while maintaining an analytically tractable reduced-order system for the subsequent MMS-based analysis.

By substituting Eqs. (5) and (13) into Eq. (10) together with Eq. (11), and then solving the resulting equation, $\psi$ can be expressed as follows:

$\psi ={\psi }_1\hspace{0.17em}\cos \left(\frac{\pi \left(ay-3bx\right)}{ab}\right)+{\psi }_2\hspace{0.17em}\cos \left(\frac{\pi \left(ay+3bx\right)}{ab}\right)+{\psi }_3\hspace{0.17em}\cos \left(\frac{2\pi \left(ay-2bx\right)}{ab}\right)$

$+{\psi }_4\hspace{0.17em}\cos \left(\frac{2\pi \left(ay-bx\right)}{ab}\right)+{\psi }_5\hspace{0.17em}\cos \left(\frac{2\pi \left(ay+bx\right)}{ab}\right)+{\psi }_6\hspace{0.17em}\cos \left(\frac{2\pi \left(ay+2bx\right)}{ab}\right)$

$+{\psi }_7\hspace{0.17em}\cos \left(\frac{\pi \left(3ay-bx\right)}{ab}\right)+{\psi }_8\hspace{0.17em}\cos \left(\frac{3\pi ay+\pi bx}{ab}\right)+{\psi }_9\hspace{0.17em}\cos \left(\frac{4\pi \left(ay-bx\right)}{ab}\right)$

$+{\psi }_{10}\hspace{0.17em}\cos \left(\frac{2\pi \left(2ay-bx\right)}{ab}\right)+{\psi }_{11}\hspace{0.17em}\cos \left(\frac{2\pi \left(2ay+bx\right)}{ab}\right)+{\psi }_{12}\hspace{0.17em}\cos \left(\frac{4\pi \left(ay+bx\right)}{ab}\right)$

$+{\psi }_{13}\hspace{0.17em}\cos \left(\frac{2\pi x}{a}\right)+{\psi }_{14}\hspace{0.17em}\cos \left(\frac{6\pi x}{a}\right)+{\psi }_{15}\hspace{0.17em}\cos \left(\frac{2\pi y}{b}\right)+{\psi }_{16}\hspace{0.17em}\cos \left(\frac{6\pi y}{b}\right)$

(14)

in which, the coefficients ${\psi }_i\left(i=1,2,\dots ,16\right)$ are obtained by substituting the assumed modal functions into the compatibility equation and collecting the corresponding harmonic terms. These coefficients depend on the geometric parameters $\pi$ , $a$ , and $b$ , the stiffness coefficients $A_{ij}$ and $B_{ij}$ , and the generalized modal amplitudes $W_i\left(i=1,2\right)$ . Because of their lengthy algebraic forms, the complete expressions are provided in the Supplementary Material, while the derivation procedure is summarized here to maintain readability. By substituting the expressions for the resultant forces and moments (written in terms of the strain relations and Eqs. (12), and (13)) into Eqs. (9c)-(9e), and subsequently applying Galerkin’s approach over $0\le x\le a$ and $0\le y\le b$ , the nonlinear discretized equations of motion (EM) for obliquely stiffened PFG-DCSSs are derived as follows:

$\begin{array}{l}{\ddot{W}}_1+2c{\dot{W}}_1+F_{11}{W}_1+F_{12}{W}_1{W}_2+F_{13}{W}_1^2+F_{14}{W}_2^2+F_{15}{W}_2^2{W}_1+F_{16}{W}_1^3+F_{17}{W}_1{\mathrm{\Phi }}_{x1}\\ +F_{18}{W}_1{\mathrm{\Phi }}_{x2}+F_{19}{W}_2{\mathrm{\Phi }}_{x1}+F_{110}{W}_2{\mathrm{\Phi }}_{x2}+F_{111}{W}_1{\mathrm{\Phi }}_{y1}+F_{112}{W}_1{\mathrm{\Phi }}_{y2}+F_{113}{W}_2{\mathrm{\Phi }}_{y1}+F_{114}{W}_2{\mathrm{\Phi }}_{y2}+F_{115}{\mathrm{\Phi }}_{x1}+F_{116}{\mathrm{\Phi }}_{y1}+F_{117}{q}_1=0\end{array}$

(15a)

$\begin{array}{l}{\ddot{W}}_2+2c{\dot{W}}_2+F_{21}{W}_2+F_{22}{W}_1{W}_2+F_{23}{W}_1^2+F_{24}{W}_2^2+F_{25}{W}_1^2{W}_2+F_{26}{W}_2^3+F_{27}{W}_1{\mathrm{\Phi }}_{x1}\\ +F_{28}{W}_1{\mathrm{\Phi }}_{x2}+F_{29}{W}_2{\mathrm{\Phi }}_{x1}+F_{210}{W}_2{\mathrm{\Phi }}_{x2}+F_{211}{W}_1{\mathrm{\Phi }}_{y1}+F_{212}{W}_1{\mathrm{\Phi }}_{y2}+F_{213}{W}_2{\mathrm{\Phi }}_{y1}+F_{214}{W}_2{\mathrm{\Phi }}_{y2}+F_{215}{\mathrm{\Phi }}_{x2}+F_{216}{\mathrm{\Phi }}_{y2}+F_{217}{q}_2=0\end{array}$

(15b)

$F_{31}{W}_1+F_{32}{W}_1{W}_2+F_{33}{W}_1^2+F_{34}{W}_2^2+F_{35}{\mathrm{\Phi }}_{x1}+F_{36}{\mathrm{\Phi }}_{y1}={\overline{\rho }}_2{\ddot{\mathrm{\Phi }}}_{x1}$

(15c)

$F_{41}{W}_2+F_{42}{W}_1{W}_2+F_{43}{W}_1^2+F_{44}{W}_2^2+F_{45}{\mathrm{\Phi }}_{x2}+F_{46}{\mathrm{\Phi }}_{y2}={\overline{\rho }}_2{\ddot{\mathrm{\Phi }}}_{x2}$

(15d)

$F_{51}{W}_1+F_{52}{W}_1{W}_2+F_{53}{W}_1^2+F_{54}{W}_2^2+F_{55}{\mathrm{\Phi }}_{x1}+F_{56}{\mathrm{\Phi }}_{y1}={\overline{\rho }}_3{\ddot{\mathrm{\Phi }}}_{y1}$

(15e)

$F_{61}{W}_2+F_{62}{W}_1{W}_2+F_{63}{W}_1^2+F_{64}{W}_2^2+F_{65}{\mathrm{\Phi }}_{x2}+F_{66}{\mathrm{\Phi }}_{y2}={\overline{\rho }}_3{\ddot{\mathrm{\Phi }}}_{y2}$

(15f)

in which, the coefficients $F_{ij}\hspace{0.17em}\left(i=1,2;j=1,2,\dots ,17\right)$ and $F_{ij}\hspace{0.17em}\left(i=3,4,\dots ,6;j=1,2,\dots ,6\right)$ appearing in Eq. (15) are obtained by substituting the stress resultants, moments, and assumed modal functions into Eqs. (9c)-(9e), followed by Galerkin integration over $0\le x\le a$ and $0\le y\le b$ . These coefficients are functions of the geometric parameters, stiffness terms $A_{ij}$ , $B_{ij}$ , and $D_{ij}$ , inertia terms, and curvature parameters of the shell. The subsequent reduction procedure used to obtain the coefficients $f_{ij}$ is described in Section 5.1, and the explicit expressions of the reduced coefficients are provided in the Appendix.

5.1. Model simplification assumptions

To obtain a tractable reduced-order model for the subsequent nonlinear resonance analysis, the following model simplification assumptions are introduced. First, the transverse modal amplitudes $W_1\left(t\right)$ and $W_2\left(t\right)$ are retained as the dominant generalized coordinates because they represent the main interacting modes associated with the considered PR and 1:2 IR conditions. Second, the rotational variables ${\mathrm{\Phi }}_{x1}\left(t\right)$ , ${\mathrm{\Phi }}_{x2}\left(t\right)$ , ${\mathrm{\Phi }}_{y1}\left(t\right)$ , and ${\mathrm{\Phi }}_{y2}\left(t\right)$ are treated as dependent variables associated with transverse shear deformation rather than independent inertia-dominated generalized coordinates. Therefore, the inertial effects associated with these rotational motions are assumed to be small compared with the transverse inertia of the shell, and the angular acceleration terms ${\ddot{\mathrm{\Phi }}}_{x1}$ , ${\ddot{\mathrm{\Phi }}}_{x2}$ , ${\ddot{\mathrm{\Phi }}}_{y1}$ , and ${\ddot{\mathrm{\Phi }}}_{y2}$ in Eqs. (15c)-(15f) are neglected. As a result, Eqs. (15c)-(15f) are reduced to a set of algebraic equations for the rotational variables. Solving these equations gives ${\mathrm{\Phi }}_{x1}\left(t\right)$ , ${\mathrm{\Phi }}_{x2}\left(t\right)$ , ${\mathrm{\Phi }}_{y1}\left(t\right)$ , and ${\mathrm{\Phi }}_{y2}\left(t\right)$ as functions of the transverse modal amplitudes $W_1\left(t\right)$ and $W_2\left(t\right)$ . The reduced coefficients $f_{ij}$ are then obtained by substituting these expressions into Eqs. (15a) and (15b) and collecting the corresponding linear, quadratic, and cubic nonlinear terms. Therefore, the coefficients $f_{ij}$ are not independent fitting parameters; they are analytically determined from the Galerkin coefficients $F_{ij}$ , stiffness coefficients, inertia terms, geometric parameters, curvature terms, and the assumed mode shapes. Accordingly, substituting these expressions into Eqs. (15a) and (15b) leads to the following two coupled nonlinear vibration equations for the obliquely stiffened PFG-DCSSs under external excitations:

${\ddot{W}}_1+2c{\dot{W}}_1+{\omega }_1^2{W}_1+f_{12}{W}_1^2{W}_2+f_{13}{W}_2^3+f_{14}{W}_1{W}_2^2+f_{15}{W}_1^3+f_{16}{W}_2^2+f_{17}{W}_1{W}_2+f_{18}{W}_1^2+f_{110}{q}_1=0$

(16a)

${\ddot{W}}_2+2c{\dot{W}}_2+{\omega }_2^2{W}_2+f_{22}{W}_1^2{W}_2+f_{23}{W}_2^3+f_{24}{W}_1{W}_2^2+f_{25}{W}_1^3+f_{26}{W}_2^2+f_{27}{W}_1{W}_2+f_{28}{W}_1^2+f_{210}{q}_2=0$

(16b)

The coefficients $f_{ij}$ appearing in Eq. (16) are reduced nonlinear coefficients obtained after eliminating the rotational variables ${\mathrm{\Phi }}_{x1}$ , ${\mathrm{\Phi }}_{x2}$ , ${\mathrm{\Phi }}_{y1}$ , and ${\mathrm{\Phi }}_{y2}$ from Eqs. (15c)-(15f) and substituting the resulting expressions into Eqs. (15a) and (15b). These coefficients are not independent fitting parameters; rather, they are analytically determined by the Galerkin coefficients $F_{ij}$ , the stiffness coefficients, inertia terms, geometric parameters, curvature terms, and the assumed mode shapes. To improve the reproducibility of the reduced-order nonlinear formulation, the explicit expressions of the coefficients $f_{ij}$ and ${\omega }_i\hspace{0.17em}\left(i=1,2\right)$ used in Eq. (16) are provided in the Appendix. Also, $q_i=Q_i\hspace{0.17em}\cos \mathrm{\Omega }t\left(i=1,2\right)$ , where $Q_i\hspace{0.17em}\left(i=1,2\right)$ denote the amplitudes of external excitation. To formulate the dimensionless EM for the obliquely stiffened PFG-DCSSs, relevant parameters and variables appearing in Eqs. (16a) and (16b) are defined below:

$\tau ={\omega }_{1}t;\hspace{0.17em}{W}_i=X_{i}h\hspace{0.17em}\left(i=1,2\right);\hspace{0.17em}\overline{\mathrm{\Omega }}=\frac{\mathrm{\Omega }}{{\omega }_1};\hspace{0.17em}\overline{C}=\frac{c}{{\omega }_1};\hspace{0.17em}{\overline{\omega }}_1=\frac{{\omega }_1}{{\omega }_1};{\overline{\omega }}_2=\frac{{\omega }_2}{{\omega }_1};\hspace{0.17em}{\overline{Q}}_1=\frac{G_{110}{Q}_1}{h{\omega }_1^2}$

${\overline{Q}}_2=\frac{G_{210}{Q}_2}{h{\omega }_1^2};\hspace{0.17em}{\overline{f}}_{ij}=\frac{f_{ij}h}{{\omega }_1^2}\hspace{0.17em}\left(i=1,2;j=6,7,8\right);\hspace{0.17em}{\overline{G}}_{ij}=\frac{G_{ij}{h}^2}{{\omega }_1^2}\hspace{0.17em}\left(i=1,2;j=2,3,\dots ,5\right)$

(17)

Using Eq. (17), Eqs. (16a) and (16b) are rewritten in the dimensionless form below:

${\ddot{X}}_1+2\overline{C}{\dot{X}}_1+{\overline{\omega }}_1^2{X}_1+{\overline{f}}_{12}{X}_1^2{X}_2+{\overline{f}}_{13}{X}_2^3+{\overline{f}}_{14}{X}_1{X}_2^2+{\overline{f}}_{15}{X}_1^3+{\overline{f}}_{16}{X}_2^2+{\overline{f}}_{17}{X}_1{X}_2+{\overline{f}}_{18}{X}_1^2+{\overline{Q}}_1\hspace{0.17em}\cos \hspace{0.17em}{\overline{\mathrm{\Omega }}}_2\tau =0$

(18a)

${\ddot{X}}_2+2\overline{C}{\dot{X}}_2+{\overline{\omega }}_2^2{X}_2+{\overline{f}}_{22}{X}_1^2{X}_2+{\overline{f}}_{23}{X}_2^3+{\overline{f}}_{24}{X}_1{X}_2^2+{\overline{f}}_{25}{X}_1^3+{\overline{f}}_{26}{X}_2^2+{\overline{f}}_{27}{X}_1{X}_2+{\overline{f}}_{28}{X}_1^2+{\overline{Q}}_2\hspace{0.17em}\cos \hspace{0.17em}{\overline{\mathrm{\Omega }}}_2\tau =0$

(18b)

Here, (^.) and (^..) indicate the first- and second-order derivatives in relation to $\tau$ , respectively.

6. Perturbation technique

To analyze the PR and 1:2 IR responses of the obliquely stiffened PFG-DCSSs, the MMSs is employed within a weakly nonlinear framework. In this approach, a small perturbation parameter $ϵ\ll 1$ is introduced to scale the displacement amplitudes, damping, and external excitation. Therefore, the present perturbation solution is valid for small-amplitude resonance responses, light damping, and weak forcing conditions. This assumption is consistent with the objective of the present work, which focuses on PR and 1:2 IR. The MMS is particularly appropriate for this purpose because it introduces separate fast and slow time scales and enables the secular terms to be identified and eliminated systematically. As a result, the method directly yields averaged modulation equations that describe the evolution of modal amplitudes and phases under detuning, damping, excitation, and nonlinear modal coupling. This feature is especially useful for the present study, since the main objective is to clarify the coupled PR and 1:2 IR mechanisms rather than only to obtain an approximate time-domain solution of Eqs. (18a) and (18b). It should also be noted that the MMS formulation is not intended to describe arbitrary strong nonlinear responses. For strongly nonlinear regimes, complementary approaches such as the HPM may be more suitable and can be considered as an extension of the present analytical framework. Accordingly, the quantities $\overline{C}$ , $X_i$ , and ${\overline{Q}}_i\hspace{0.17em}\left(i=1,2\right)$ are defined as follows:

$\stackrel{\sim }{C}={ϵ}^2\overline{C};\hspace{0.17em}{X}_i=ϵ{\stackrel{\sim }{X}}_i;\hspace{0.17em}{\overline{Q}}_i={ϵ}^2{\stackrel{\sim }{Q}}_i\hspace{0.17em}\left(i=1,2\right)$

(19)

Upon substituting Eq. (19) into Eqs. (18a) and (18b) and applying the MMSs, the approximate solutions can be expressed in the following form:

${\stackrel{\sim }{X}}_i\left(\tau ,\epsilon \right)=Z_{i0}\left(T_0,T_1,T_2\right)+ϵZ_{i1}\left(T_0,T_1,T_2\right)+{ϵ}^2{Z}_{i2}\left(T_0,T_1,T_2\right)+\dots ;\hspace{0.17em}\left(i=1,2\right)$

(20)

in which, $T_i={ϵ}^i\tau \hspace{0.17em}\left(i=0,1,2\right)$ . Since the MMS solution is constructed as an asymptotic expansion, the accuracy of the approximation depends on the smallness of the perturbation parameter $ϵ$ and on the consistent truncation of higher-order terms. In the present formulation, terms are retained up to the order required to capture the PR and 1:2 IR interactions. The neglected higher-order terms represent the perturbation truncation error and are assumed to remain small under the weakly nonlinear conditions adopted in this study. Therefore, the obtained averaged equations are valid for small-amplitude responses, light damping, and weak external excitation, which are consistent with the scaling assumptions introduced in Eq. (19). Also, using Eqs. (18a) and (18b), the following differential operators are obtained:

$\frac{d}{dt}=\frac{\partial }{\partial T_0}\frac{\partial T_0}{\partial \tau }+\frac{\partial }{\partial T_1}\frac{\partial T_1}{\partial \tau }+\frac{\partial }{\partial T_2}\frac{\partial T_2}{\partial \tau }+\dots =D_0+\epsilon D_1+{\epsilon }^2{D}_2+\dots$

$\frac{d^2}{{dt}^2}={\left(D_0+\epsilon D_1+{\epsilon }^2{D}_2+\dots \right)}^2=D_0^2+2\epsilon D_0{D}_1+{\epsilon }^2\left(D_1^2+2D_0{D}_2\right)+\dots$

(21)

This study focuses on PR and 1:2 IR. It is also worth noting that, in the present work, PR occurs when $\bar{\mathrm{\Omega }}$ is close to ${\bar{\omega }}_2$ . Accordingly, the resonance conditions can be written as follows:

${\overline{\omega }}_1^2=\frac{{\overline{\mathrm{\Omega }}}^2}{4}-{ϵ}^2{\sigma }_1;\hspace{0.17em}{\overline{\omega }}_2^2={\overline{\mathrm{\Omega }}}^2-{ϵ}^2{\sigma }_2;\hspace{0.17em}\overline{\mathrm{\Omega }}=2$

(22)

in which, ${\sigma }_i\hspace{0.17em}\left(i=1,2\right)$ represent the detuning parameters. Substituting Eqs. (19)-(22) into Eqs. (18a) and (18b), and then equating the coefficients of ${ϵ}^i\hspace{0.17em}\left(i=0,1,2\right)$ to zero, gives:

• • Order ${ϵ}^0$ :

$D_0^2{Z}_{10}+Z_{10}=0$

(23a)

$D_0^2{Z}_{20}+4Z_{20}=0$

(23b)

• • Order $ϵ$ :

$D_0^2{Z}_{11}+Z_{11}=-2D_0{D}_1{Z}_{10}-{\overline{f}}_{18}{Z}_{10}^2-{\overline{f}}_{17}{Z}_{10}{Z}_{20}-{\overline{f}}_{16}{Z}_{20}^2-{\stackrel{\sim }{Q}}_1\hspace{0.17em}\cos \left(2\tau \right)$

(24a)

$D_0^2{Z}_{21}+4Z_{21}=-2D_0{D}_1{Z}_{20}-{\overline{f}}_{28}{Z}_{10}^2-{\overline{f}}_{27}{Z}_{10}{Z}_{20}-{\overline{f}}_{26}{Z}_{20}^2-{\stackrel{\sim }{Q}}_2\hspace{0.17em}\cos \left(2\tau \right)$

(24b)

• • Order ${ϵ}^2$ :

$D_0^2{Z}_{12}+Z_{12}=-2D_0{D}_1{Z}_{11}-D_1^2{Z}_{10}-2D_0{D}_2{Z}_{10}-2\overline{C}{D}_0{Z}_{10}-{\sigma }_1{Z}_{10}-{\overline{f}}_{15}{Z}_{10}^3-{\overline{f}}_{12}{Z}_{10}^2{Z}_{20}-2{\overline{f}}_{18}{Z}_{10}{Z}_{11}-{\overline{f}}_{17}\left(Z_{10}{Z}_{21}+Z_{11}{Z}_{20}\right)-2{\overline{f}}_{16}{Z}_{20}{Z}_{21}-{\overline{f}}_{14}{Z}_{20}^2{Z}_{10}-{\overline{f}}_{13}{Z}_{20}^3$

(25a)

$D_0^2{Z}_{22}+4Z_{22}=-2D_0{D}_1{Z}_{21}-D_1^2{Z}_{20}-2D_0{D}_2{Z}_{20}-2\overline{C}{D}_0{Z}_{20}-{\sigma }_1{Z}_{20}-{\overline{f}}_{25}{Z}_{10}^3-{\overline{f}}_{22}{Z}_{10}^2{Z}_{20}-2{\overline{f}}_{28}{Z}_{10}{Z}_{11}-{\overline{f}}_{27}\left(Z_{10}{Z}_{21}+Z_{11}{Z}_{20}\right)-2{\overline{f}}_{26}{Z}_{20}{Z}_{21}-{\overline{f}}_{24}{Z}_{20}^2{Z}_{10}-{\overline{f}}_{23}{Z}_{20}^3$

(25b)

The solution form of Eq. (23) is:

$Z_{j0}=A_1{e}^{ji{T}_0}+{\overline{A}}_1{e}^{-ji{T}_0};\hspace{0.17em}\left(j=1,2\right)$

(26)

Inserting Eq. (26) into Eq. (23) gives:

$D_0^2{Z}_{11}+Z_{11}=-\left(2i{D}_1{A}_1+{\overline{f}}_{17}{A}_2{\overline{A}}_1\right)e^{i{T}_0}+c.c.+NST$

$D_0^2{Z}_{21}+{4Z}_{21}=-\left(4i{D}_1{A}_2+{\overline{f}}_{28}{A}_1^2+\frac{1}{2}{\stackrel{\sim }{Q}}_2\right)e^{2i{T}_0}+c.c.+NST$

(27)

Here, c.c. denotes complex conjugate, while NST stands for the non-secular terms. By setting secular terms in Eq. (27) to zero, the following complex-form equations are obtained:

$D_1{A}_1=\frac{1}{2}{\overline{f}}_{17}{A}_2{\overline{A}}_1$

(28a)

$D_1{A}_2=\frac{1}{4}{\overline{f}}_{28}{A}_1^2+\frac{1}{8}{\stackrel{\sim }{Q}}_2$

(28b)

Accordingly, the solution of Eq. (27) is given by:

$Z_{11}=\frac{1}{15}{\overline{f}}_{16}{A}_2^2{e}^{4i{T}_0}+\frac{1}{8}{\overline{f}}_{17}{A}_1{A}_2{e}^{3i{T}_0}+\left(\frac{1}{3}{\overline{f}}_{18}{A}_1^2-\frac{1}{6}{\stackrel{\sim }{Q}}_1\right)e^{2i{T}_0}-2{\overline{f}}_{18}{A}_1{\overline{A}}_1-2{\overline{f}}_{16}{A}_2{\overline{A}}_2+c.c.$

(29a)

$Z_{21}={\frac{1}{12}{\overline{f}}_{26}{A}_2^{2}e}^{4i{T}_0}+\frac{1}{5}{\overline{f}}_{27}{A}_1{A}_2{e}^{3i{T}_0}-\frac{1}{3}{\overline{f}}_{27}{\overline{A}}_1{A}_2{e}^{i{T}_0}-\frac{1}{2}{{\overline{f}}_{28}\overline{A}}_1{A}_1-\frac{1}{2}{{\overline{f}}_{26}\overline{A}}_2{A}_2+c.c.$

(29b)

Inserting Eqs. (26) and (29) into Eq. (25) and eliminating the secular terms by setting them to zero gives:

$D_2{A}_1=\left({\overline{f}}_{14}+\frac{1}{16}{\overline{f}}_{17}^2-\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{26}-2{\overline{f}}_{18}{\overline{f}}_{26}-\frac{2}{15}{\overline{f}}_{16}{\overline{f}}_{17}\right)i{A_1{\overline{A}}_{2}A}_2+\left(\frac{3}{2}{\overline{f}}_{15}-\frac{5}{3}{\overline{f}}_{18}^2-\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{28}\right)i{A}_1^2{\overline{A}}_1+\left(-\frac{1}{2}i{\sigma }_1-\overline{C}\right)A_1-\frac{1}{6}{\overline{f}}_{18}{\stackrel{\sim }{Q}}_1{\overline{A}}_1$

(30a)

$D_2{A}_2=\left(\frac{1}{2}{\overline{f}}_{12}-\frac{1}{30}{\overline{f}}_{27}^2-\frac{1}{4}{\overline{f}}_{28}{\overline{f}}_{26}+\frac{1}{16}{\overline{f}}_{17}{\overline{f}}_{28}-\frac{1}{2}{\overline{f}}_{18}{\overline{f}}_{27}\right)i{A_1{\overline{A}}_{1}A}_2-\frac{5}{48}{f}_{26}{\overline{f}}_{29}i{A}_2{\overline{A}}_2+\left(\frac{3}{4}{\overline{f}}_{23}-\frac{5}{24}{\overline{f}}_{26}^2-\frac{29}{60}{\overline{f}}_{26}{\overline{f}}_{16}\right)i{A}_2^2{\overline{A}}_2-\left(\frac{1}{4}i{\sigma }_2+\overline{C}\right)A_2$

(30b)

The amplitude functions can be expressed in Cartesian form ( $A_{i+1}=x_{2i+1}+i{x}_{2i+2}\hspace{0.17em};\hspace{0.17em}\left(i=0,1\right)$ ). Accordingly, we adopt this representation and compute the response for each state. As a result, decomposing the governing equation into its real and imaginary components yields:

$\frac{d{x}_1}{d{T}_2}=\frac{-\overline{C}{x}_1}{2}+\frac{{\sigma }_1{x}_2}{2}-\frac{1}{6}{\overline{f}}_{18}{\stackrel{\sim }{Q}}_1{x}_2-\frac{1}{8}{\overline{f}}_{17}\left(\frac{1}{2}{\overline{f}}_{19}-{\overline{f}}_{29}\right)x_2{x}_3+\left(2{\overline{f}}_{18}{\overline{f}}_{16}-\frac{1}{16}{\overline{f}}_{17}^2+\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{26}+\frac{2}{15}{\overline{f}}_{16}{\overline{f}}_{27}\right)x_2\left(x_3^2+x_4^2\right)+\left(\frac{5}{3}{\overline{f}}_{18}^2+\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{28}-\frac{3}{2}{\overline{f}}_{15}\right)\left(x_2{x}_1^2+x_2^3\right)$

(31a)

$\frac{d{x}_2}{d{T}_2}=\frac{-\overline{C}{x}_2}{2}-\frac{{\sigma }_1{x}_1}{2}-\frac{1}{6}{\overline{f}}_{18}{\stackrel{\sim }{Q}}_1{x}_1+\frac{1}{8}{\overline{f}}_{17}\left(\frac{1}{2}{\overline{f}}_{19}-{\overline{f}}_{29}\right)x_1{x}_3-\left(2{\overline{f}}_{18}{\overline{f}}_{16}-\frac{1}{16}{\overline{f}}_{17}^2+\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{26}+\frac{2}{15}{\overline{f}}_{16}{\overline{f}}_{27}\right)x_1\left(x_3^2+x_4^2\right)-\left(\frac{5}{3}{\overline{f}}_{18}^2+\frac{1}{4}{\overline{f}}_{17}{\overline{f}}_{28}-\frac{3}{2}{\overline{f}}_{15}\right)\left(x_1{x}_2^2+x_1^3\right)$

(31b)

$\frac{d{x}_3}{d{T}_2}=\frac{-\overline{C}{x}_3}{2}+\frac{{\sigma }_2{x}_4}{4}+\left(\frac{1}{2}{\overline{f}}_{18}{\overline{f}}_{27}-\frac{1}{16}{\overline{f}}_{17}{\overline{f}}_{28}+\frac{1}{4}{\overline{f}}_{28}{\overline{f}}_{26}+\frac{1}{30}{\overline{f}}_{27}^2-\frac{1}{6}{\overline{f}}_{22}\right)x_4\left(x_1^2+x_2^2\right)+\left(\frac{5}{24}{\overline{f}}_{26}^2+\frac{29}{60}{\overline{f}}_{16}{\overline{f}}_{27}-\frac{3}{4}{\overline{f}}_{23}\right)\left(x_4{x}_3^2+x_4^3\right)$

(31c)

$\frac{d{x}_4}{d{T}_2}=\frac{-\overline{C}{x}_4}{2}-\frac{{\sigma }_2{x}_3}{4}-\left(\frac{1}{2}{\overline{f}}_{18}{\overline{G}}_{27}-\frac{1}{16}{\overline{f}}_{17}{\overline{f}}_{28}+\frac{1}{4}{\overline{f}}_{28}{\overline{f}}_{26}+\frac{1}{30}{\overline{f}}_{27}^2-\frac{1}{2}{\overline{f}}_{22}\right)x_3\left(x_1^2+x_2^2\right)-\left(\frac{5}{24}{\overline{f}}_{26}^2+\frac{29}{60}{\overline{f}}_{16}{\overline{f}}_{27}-\frac{3}{4}{\overline{f}}_{23}\right)\left(x_3{x}_4^2+x_3^3\right)$

(31d)

7. Numerical results

This section examines the ND response of obliquely stiffened PFG-DCSSs through waveform characteristics, phase portraits, and Poincaré maps, highlighting how changes in stiffener angle and material properties affect the shells’ nonlinear behavior. The materials considered are Aluminum ( ${\rho }_m=2702\hspace{0.17em}\text{kg}/{\mathrm{m}}^3$ , $E_m=70\hspace{0.17em}\text{GPa}$ ) and Alumina ( ${\rho }_c=3800\hspace{0.17em}\text{kg}/{\mathrm{m}}^3$ , $E_c=380\hspace{0.17em}\text{GPa}$ ). The geometric specifications of the obliquely stiffened PFG-DCSSs are summarized in Table 1. Since $R_x/R_y=1$ , the numerical shell considered in the main parametric study corresponds to a spherical shallow shell. Nevertheless, the theoretical formulation is expressed in terms of the two principal radii of curvature, $R_x$ and $R_y$ , and therefore remains applicable to general doubly curved shallow shells by selecting appropriate curvature ratios. Before presenting the numerical results, all geometric parameters used in the considered cases are defined here for clarity and reproducibility. The shell thickness is $h=0.002$ m, and the in-plane dimensions are $a=b=1$ m. The radii of curvature are $R_x=R_y=1.5$ m, satisfying $R_x/R_y=1$ and $R_x/a=1.5$ . The oblique stiffeners have thickness $h_s=0.01$ m and width $d=2.5$ mm. The number of stiffeners is $N_s=50$ , and the corresponding stiffener spacing is $S=\left(a-N_{s}d\right)/\left(N_s+1\right)$ for each stiffener family. Unless otherwise stated, these geometric parameters are used for all numerical cases, while the stiffener angles $\alpha$ and $\beta$ are varied to examine their influence on the ND response. All results reported in Subsection 7.1 correspond to obliquely stiffened PFG-DCSSs exhibiting these mode shapes.

Table 1.

Geometric characteristics of OSs and DCSSs.

$h\left(\mathrm{m}\right)$	$a/b$	$R_x/R_y$	$R_x/a_y$	$\nu$	$h_s\left(\mathrm{m}\right)$	$d\left(\text{mm}\right)$	$N_{DCSSs}$	$N_{OSs}$	$N_S$ *
0.002	1	$1$	$1.5$	0.3	0.01	2.5	1	1	50

^*Number of PFG-OSs.

The four-dimensional averaged equations given in Eq. (31) were solved numerically using the classical fourth-order Runge-Kutta method. Since Eq. (31) is formulated in terms of the MMS slow time scale, the numerical integration was performed with respect to the dimensionless slow time variable $T_2$ , where $T_2={ϵ}^2\tau$ and $\tau ={\omega }_{1}t$ . A fixed time step of $\mathrm{\Delta }{T}_2=0.001$ was used was used in the numerical integration. It should be noted that $T_2$ represents the dimensionless slow time and should not be interpreted as the original dimensional physical time $t$ . The system was integrated over the interval $T_2=0$ to $T_2=2000$ . To remove the influence of initial transient behavior, the initial portion of the response was discarded, and only the steady-state interval $T_2=1900$ to $T_2=2000$ was used to generate the time histories, phase portraits, and Poincaré maps. Accordingly, the horizontal axis in the waveform diagrams represents the dimensionless slow time $T_2$ . Unless otherwise stated, the initial conditions were taken as $x_1\left(0\right)=0$ , $x_2\left(0\right)=0$ , $x_3\left(0\right)=0$ , and $x_4\left(0\right)=0$ .

7.1. Verification, numerical accuracy, and convergence assessment

The geometric and material parameters used in the numerical analysis are summarized in Table 1, while their definitions are provided in Section 2. The mode shapes used in the present reduced-order formulation are defined by the admissible functions introduced in Eq. (13). In Figure 3, $m$ and $n$ denote the half-wave numbers in the $x$ - and $y$ -directions, respectively. The first retained mode corresponds to $m=1$ and $n=3$ , and the second retained mode corresponds to $m=3$ and $n=1$ . These two modes are selected to represent the dominant interacting modal components associated with the considered PR and 1:2 IR conditions.

Figure 3.

PFG-DCSS mode shapes used in the present analysis.

To assess the reliability of the proposed approach, the outcomes of this study were compared with those from previous works. In particular, the non-dimensional natural frequency (NDNF) of cylindrical shells (CSs) was evaluated in comparison with the results of Song and Li, ⁵⁰ Zhang et al., ⁵¹ and Yang et al., ⁵² as summarized in Table 2. As shown in Table 2, most of the present results are in very close agreement with the published benchmark solutions. The maximum discrepancy of 6.2% occurs only for the relatively thicker benchmark shell case ( $h/R=0.05$ , $n=1$ ). This case is less representative of the main configuration considered in the present study, which corresponds to a thin-shell structure. For the thin-shell benchmark case ( $h/R=0.002$ ), the discrepancies are extremely small, confirming the accuracy of the present formulation in the thickness regime most relevant to this work. The relatively larger discrepancy observed for the thicker shell case may be attributed to the stronger effects of transverse shear deformation, rotary inertia, curvature terms, shear correction assumptions, and differences in the shell theories or numerical procedures adopted in the reference studies. Likewise, as an additional limiting-case verification, the NDNF of a rectangular plate was checked against the results of Malik and Bert, ⁵³ as shown in Table 3. It should be emphasized that the rectangular plate comparison is included only as a benchmark limiting case of the present formulation and does not represent the main geometry investigated in the ND analyses. A rectangular plate can be recovered from the doubly curved shallow-shell formulation by neglecting the curvature terms, i.e., by taking $R_x,R_y\to \infty$ . In addition, Table 4 provides a comparison of the natural frequency (NF) of FG-DCSSs with that reported by Bich et al. ⁵⁴ Since the present study focuses on an analytical/semi-analytical reduced-order resonance formulation, a separate nonlinear FE model is not included here; however, the benchmark comparisons in Tables 2–4 support the accuracy of the proposed formulation for the considered shell-type structures. Therefore, the validation results are considered acceptable and support the reliability of the present formulation for the vibration analysis of thin shell-type structures.

Table 2.

Comparisons of NDNF of CSs ( $m=1,\nu =0.3,L/R=20$ ).

​	Present	Song et al. 50		Yang et al. 52		Zhang et al. 51
		Discrepancy (%)		Discrepancy (%)		Discrepancy (%)
$h/R=0.002$
$n=1$	0.0161008	0.0161011	0.002	0.0161003	0.003	0.0161011	0.002
$n=3$	0.0050421	0.0050424	0.006	0.0050423	0.004	0.0050418	0.006
$h/R=0.05$
$n=1$	0.0151266	0.0161299	6.2	0.0151185	0.05	0.0161065	6.1
$n=3$	0.1097431	0.1097653	0.02	0.1096523	0.08	0.1098113	0.06

Table 3.

Comparisons of the NDNF for rectangular plate ( $m=1,n=3,a/b=0.5,h/b=0.1,\nu =0.3$ ).

Present	Malik and Bert 53
	Discrepancy (%)
4.39602	4.27211	2.818

Table 4.

Comparisons of the NF (Hz) for FG-DCSSs ( $m=1,n=3,a/b=1,b/h=20$ ).

Present	Bich et al. 54
	Discrepancy (%)
1633.89	1633.93	0.002

To further assess the accuracy and reliability of the proposed formulation, additional validation studies have been included by comparing the present results with available benchmark results from the literature. These comparisons are performed for relevant limiting cases, since direct experimental data for obliquely stiffened PFG-DCSSs under simultaneous PR and 1:2 IR are not currently available in the open literature. The good agreement observed in the validation results confirms the correctness of the theoretical formulation and the numerical implementation adopted in the present study. Nevertheless, future experimental investigations on porous FG stiffened shell structures under controlled harmonic excitation would be valuable for further evaluating the predicted periodic, quasi-periodic, and chaotic responses.

In addition to the above validation studies, the numerical accuracy and convergence of the MMS-based solution were also examined. Since the MMS formulation is based on an asymptotic perturbation expansion, the perturbation series was consistently truncated up to the retained order, and the analysis was restricted to the weakly nonlinear regime defined in Eq. (19). In this regime, the response amplitudes, damping, and external excitation are scaled by the small perturbation parameter. To check the numerical convergence of the averaged equations, the integration time step was progressively reduced, and the resulting time histories, phase portraits, and Poincaré maps were compared. The qualitative response characteristics, including periodic, quasi-periodic, and chaotic motions, remained unchanged under time-step refinement. This confirms that the reported ND responses are not artifacts of numerical discretization and supports the reliability of the MMS-based numerical results.

7.2. Results on chaotic and periodic behavior

The current study investigates the ND analyses of obliquely stiffened PFG-DCSSs, with particular attention to the strong impact of stiffener angle variations and material composition on their behavior. The analysis concentrates on the presence of 1:2 IR together with the coupling associated with the primary resonance response in obliquely stiffened PFG-DCSSs. Examining these phenomena is important for assessing the dynamic stability of shell structures in real operating conditions. Changes in stiffener orientation alter the interaction among vibration modes, which directly affects the change from periodic to chaotic motion. In a similar way, changes in material properties, namely mass density and Young’s modulus, modify the vibration characteristics of the structure and underline the important contribution of material gradation to its overall dynamic performance. To present the ND behavior of obliquely stiffened PFG-DCSSs subjected to external excitation, particularly under primary and 1:2 internal resonance conditions, several graphical tools were used, including Poincaré maps, 3D and two-dimensional (2D) phase portraits, and time-history waveforms. It should be noted that the identification of quasi-periodic and chaotic responses in this study is based on the combined interpretation of standard nonlinear-dynamics indicators, including time-history responses, phase portraits, and Poincaré maps. To provide a clearer criterion, quasi-periodic motion is identified when the response remains bounded and non-repeating, the phase trajectory forms an organized closed-band or torus-like structure, and the Poincaré map consists of ordered points arranged along smooth closed curves or invariant sets. In contrast, chaotic motion is identified when the bounded response exhibits irregular aperiodic oscillations, the phase trajectory fills a broader region of the phase space without forming a regular closed structure, and the Poincaré map shows scattered points distributed irregularly over a finite area. Therefore, the response classification is not based on a single waveform or phase portrait alone, but on the consistent behavior observed from all three diagnostic tools. Quantitative measures, such as the largest Lyapunov exponent, correlation dimension, and bifurcation diagrams, can provide additional confirmation of chaotic behavior and may be incorporated in future studies. The computational results clearly show that obliquely stiffened PFG-DCSSs can exhibit both periodic and chaotic responses. Moreover, the findings indicate that modifications in stiffener orientations and material properties may drive the system toward different types of motion, including quasi-periodic and chaotic behavior.

The chaotic responses reported in this section are obtained from the MMS-based averaged equations and should therefore be interpreted as weakly nonlinear resonance-induced chaotic modulation responses. These responses result from the interaction of modal coupling, detuning, damping, and external excitation within the assumed small-amplitude framework. Thus, the term “chaotic” is used here to describe the qualitative dynamic features observed in the reduced slow-flow system, rather than to imply large-amplitude strongly nonlinear structural motion.

In the following waveform diagrams, the horizontal axis represents the dimensionless slow time $T_2$ , not the original dimensional time $t$ . The plotted range $T_2=1900$ to $2000$ is selected to show the post-transient steady-state response of the system.

Figures 4–11 present the ND response of obliquely stiffened PFG-DCSSs with type 1 for different stiffener angles. In this case, the shell is functionally graded from ceramic-rich outer surfaces toward a more metal-rich composition through the thickness, while the porosity is considered to be uniformly distributed across the thickness. The stiffeners are also considered functionally graded members with uniform porosity throughout their thickness. The results clearly show that both stiffener orientation and material gradation have a strong influence on the shell’s dynamic behavior. More specifically, the obliquely stiffened PFG-DCSSs display quasi-periodic motion for the angle combinations $\alpha =30°$ and $\beta =0°$ or $90°$ , and $\alpha =60°$ and $\beta =60°$ or $90°$ . By contrast, chaotic motion is observed for $\alpha =0°$ and $\beta =60°$ or $90°$ , and $\alpha =30°$ and $\beta =30°$ , $60°$ , or $90°$ .

Figure 4.

Quasi-periodic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =0°$ and $\beta =30°$ .

Figure 5.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =0°$ and $\beta =60°$ .

Figure 6.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =0°$ and $\beta =90°$ .

Figure 7.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =30°$ and $\beta =30°$ .

Figure 8.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =30°$ and $\beta =60°$ .

Figure 9.

Quasi-periodic motion occurs in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =30°$ and $\beta =90°$ .

Figure 10.

Quasi-periodic motion appears in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =60°$ and $\beta =60°$ .

Figure 11.

Quasi-periodic motion appears in the obliquely stiffened PFG-DCSSs with type 1 when $\alpha =60°$ and $\beta =90°$ .

Figures 4–11 provide a comparative illustration of the ND responses of the obliquely stiffened PFG-DCSSs with type 1 porosity distribution under different stiffener angle combinations. In these cases, the shell has a uniform porosity distribution through the thickness, while the stiffener orientations are varied to examine their influence on the resonance response. The comparison is performed using time-history waveforms, two-dimensional phase portraits, three-dimensional phase trajectories, and Poincaré maps. The results show that changing the stiffener angles significantly modifies the modal interaction and the resulting nonlinear response. In quasi-periodic cases, the response remains bounded with organized phase trajectories and structured Poincaré points, whereas chaotic cases are characterized by irregular waveforms, more complex phase portraits, and scattered Poincaré points. Therefore, Figures 4–11 clearly demonstrate that, even for the same uniform porosity distribution, the stiffener orientation can shift the system from quasi-periodic motion to chaotic motion or vice versa.

Figures 12–19 illustrate the ND response of obliquely stiffened PFG-DCSSs with type 2 for various stiffener angles. In this configuration, the shell follows the same through-thickness material grading, but the porosity is distributed nonuniformly, meaning that it changes along the thickness according to a specified pattern. The stiffeners are also assumed to be functionally graded members with uniform porosity through their thickness. The results clearly indicate that changes in stiffener orientation and material gradation significantly influence the shell’s dynamic behavior. In particular, the obliquely stiffened PFG-DCSSs show quasi-periodic motion for the angle combinations $\alpha =0°$ and $\beta =90°$ , and $\alpha =30°$ and $\beta =60°$ . In contrast, chaotic motion appears for $\alpha =0°$ and $\beta =30°$ or $60°$ , and $\alpha =30°$ and $\beta =30°$ or $90°$ , and $\alpha =60°$ and $\beta =60°$ or $90°$ .

Figure 12.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =0°$ and $\beta =30°$ .

Figure 13.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =0°$ and $\beta =60°$ .

Figure 14.

Quasi-periodic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =0°$ and $\beta =90°$ .

Figure 15.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =30°$ and $\beta =30°$ .

Figure 16.

Quasi-periodic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =30°$ and $\beta =60°$ .

Figure 17.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =30°$ and $\beta =90°$ .

Figure 18.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =60°$ and $\beta =60°$ .

Figure 19.

Chaotic motion occurs in the obliquely stiffened PFG-DCSSs with type 2 when $\alpha =60°$ and $\beta =90°$ .

Figures 12–19 present the corresponding ND responses of the obliquely stiffened PFG-DCSSs with type 2 porosity distribution, where the porosity varies nonuniformly through the shell thickness. Similar to the previous set of results, different stiffener angle combinations are considered while the other material, geometric, and excitation parameters are kept unchanged. The numerical illustrations show that the nonuniform porosity distribution changes the stiffness and inertia characteristics of the shell, which in turn affects the modal coupling and resonance behavior. Compared with the uniform porosity case, the regions of quasi-periodic and chaotic motion are redistributed, indicating that porosity distribution plays an important role in the nonlinear stability characteristics of the system. In the quasi-periodic responses, the Poincaré maps exhibit organized patterns, while in the chaotic responses, the points are scattered and the phase portraits become more irregular. Thus, Figures 12–19 confirm that both stiffener orientation and porosity distribution strongly influence the transition between quasi-periodic and chaotic responses in obliquely stiffened PFG-DCSSs.

8. Conclusions

In this work, the ND response of obliquely stiffened PFG-DCSSs under external excitation was examined. The governing nonlinear equations of motion were first formulated based on FSDT together with the von Kármán relations for PFG-DCSSs reinforced by two intersecting oblique stiffeners, whose angles could be identical or different. After that, Galerkin’s approach was employed to reduce the system to a two-degree-of-freedom ND model. The analysis considered PR and 1:2 IR. Subsequently, MMSs was used to derive four averaged nonlinear equations. It should be noted that the perturbation formulation used in this study is based on weakly nonlinear assumptions. Therefore, the analytical results are mainly applicable to small-amplitude resonance responses with light damping and weak external excitation. Extending the present formulation to strongly nonlinear regimes would require additional analytical or semi-analytical treatments. In this regard, the HPM represents a promising complementary approach for future investigations of strongly nonlinear responses in obliquely stiffened PFG-DCSSs. Based on the findings of the current study, the following main conclusions may be drawn:

• For the case of uniform porosity distribution, where the shell is graded from ceramic-rich outer surfaces toward a more metal-rich core, the response strongly depends on the selected stiffener orientations. The system exhibits quasi-periodic motion for the angle combinations $\alpha =30°$ and $\beta =0°$ or $90°$ , and $\alpha =60°$ and $\beta =60°$ or $90°$ . In contrast, chaotic motion is observed for $\alpha =0°$ and $\beta =60°$ or $90°$ , as well as $\alpha =30°$ and $\beta =30°$ , $60°$ , or $90°$ . These results show that even under uniform porosity conditions, changing the stiffener arrangement can significantly alter the vibration regime of the shell.

• For the case of nonuniform porosity distribution, the shell follows the same through-thickness material grading, but the porosity varies along the thickness according to a prescribed distribution pattern. Under this condition, the nonlinear response again changes noticeably with the stiffener angles. The results indicate quasi-periodic motion for $\alpha =0°$ and $\beta =90°$ , and $\alpha =30°$ and $\beta =60°$ . On the other hand, chaotic motion occurs for $\alpha =0°$ and $\beta =30°$ or $60°$ , and $\alpha =30°$ and $\beta =30°$ or $90°$ , and $\alpha =60°$ and $\beta =60°$ or $90°$ . This demonstrates that the nonuniform distribution of porosity modifies the stability characteristics and can shift the regions where chaotic behavior appears.

• Overall, the results clearly reveal that stiffener angles and porosity distribution are key parameters in controlling the ND behavior of obliquely stiffened PFG-DCSSs. By selecting suitable angles for the stiffeners, it is possible to guide the system toward more stable quasi-periodic responses or, conversely, into chaotic regimes. This provides useful insight into the design and optimization of advanced functionally graded shell structures exposed to complex dynamic environments.

The main limitations of the proposed method should also be acknowledged. The present analytical framework is based on a reduced-order Galerkin model combined with the MMS; therefore, it is mainly suitable for weakly nonlinear resonance responses under small-amplitude vibration, weak external excitation, and light damping conditions. In addition, only selected dominant modes are retained to capture the PR and 1:2 IR interactions, while higher-mode interactions are not explicitly included. The analysis is also limited to simply supported boundary conditions, perfect bonding between the shell and stiffeners, an equivalent-continuum porosity model, and a smeared stiffener representation. Therefore, the proposed method is appropriate for global ND and resonance analysis, but it does not explicitly capture local stiffener–shell junction effects, stress concentrations, local buckling, interfacial damage, or strongly nonlinear large-amplitude responses.

It should also be noted that the chaotic responses reported in this study are obtained within the weakly nonlinear MMS-based averaged formulation and should therefore be interpreted as resonance-induced chaotic modulation responses under small-amplitude vibration, weak excitation, and light damping conditions. Moreover, the present analysis is carried out for simply supported boundary conditions; therefore, the quantitative resonance responses may change under clamped, free, or elastically restrained edges, and investigating the sensitivity of the ND response to different boundary conditions is recommended as a future extension of this work. In addition, the oblique stiffeners were modeled using a smeared stiffener technique, which is appropriate for global ND analysis of regularly distributed stiffened shells but does not explicitly resolve local stiffener-shell junction effects, stress concentrations, or local buckling of individual stiffeners.

Furthermore, the present study adopts a widely used equivalent-continuum porosity model to quantify the porosity distribution of FG doubly curved shallow shells. In this model, uniform and nonuniform porosity distributions are represented through prescribed through-thickness material-property functions, which enables the global influence of porosity on the ND response to be examined. However, the detailed pore morphology, pore connectivity, and two-scale fractal characteristics of the porous medium are not explicitly considered. Future studies may extend the present formulation by introducing two-scale fractal dimensions into the effective material-property model to provide a more refined description of porous microstructures and their influence on structural mechanical behavior.

Finally, the present model assumes a continuously attached, ceramic-rich common/contact layer between the PFG-OSs and the PFG-DCSSs. Accordingly, the stiffeners and shell are treated as perfectly bonded, and no relative displacement, sliding, separation, interfacial damage, or delamination is considered at the shell-stiffener interface. Since imperfect bonding or delamination may reduce the effective stiffness transfer between the stiffeners and the shell and consequently influence the ND and resonance responses, future studies may extend the present formulation by incorporating imperfect interface conditions or delamination models.