Constrained dynamic optimization of bidirectional porous functionally graded beams using Touratier higher-order shear deformation theory

3.1. Displacement field

According to Touratier’s HSDT [26], as adapted to porous FG beam dynamics in recent higher-order formulations [13, 36], the displacement field is written as:

where $u_0$ is the mid-plane axial displacement, $w_b$ is the bending component of the transverse displacement, and $w_s$ is the shear component of the transverse displacement. The Touratier shear shape function is taken as:

and its derivative gives:

The axial and transverse shear strains are therefore:

3.2. Energy expressions and stiffness coefficients

Hamilton’s principle, following the variational derivations commonly adopted in porous FG wave-propagation models [27, 29, 30], is used to obtain the governing equations:

The strain energy is:

where $A$ is the cross-section and $G\left(y,z\right)=E\left(y,z\right)/\left[2\left(1+\nu \right)\right]$ assuming constant Poisson’s ratio $\nu$. The kinetic energy is:

The stiffness coefficients entering the HSDT formulation are obtained from cross-sectional integrations:

The mass coefficients are similarly defined as:

Eqs. (15-19) explicitly define the origin of the stiffness and mass matrices used later in the eigenvalue problem. These coefficients are recomputed for each design vector during optimization because $E\left(y,z\right)$ and $\rho \left(y,z\right)$ depend on $n$, $k$, and $\alpha$.

3.3. Eigenvalue problem and phase velocity

Using the generalized displacement vector:

the governing equations can be written in matrix form as:

where ${\mathbf{x}}_d={\left[n,k,\alpha ,b/h\right]}^T$ is the design vector, $\mathbf{K}$ is assembled from the stiffness coefficients in Eqs. (15-17), and $\mathbf{M}$ is assembled from the mass coefficients in Eqs. (18-19). Non-trivial solutions require:

For wave propagation analysis, the phase velocity of mode $i$ is evaluated from the dispersion relation in the standard form used for porous FG beam and plate models [27, 28, 36]:

where $\lambda$ is the wave number and ${\omega }_i$ is the corresponding circular frequency.

4.1. Design variables, objective function, and constraints

The design vector is defined as:

The four variables represent complementary mechanisms governing the dynamic response. The indices $n$ and $k$ independently control the through-thickness and width-wise material distributions, respectively; $\alpha$ modifies the local stiffness and inertia fields through porosity; and $b/h$ controls the cross-sectional geometry. Treating $n$ and $k$ independently is essential for a bidirectional FG beam because the two spatial directions enter the stiffness and mass integrals differently.

The objective is to increase the fundamental natural frequency while reducing structural mass. The normalized scalar objective is:

where $f_1={\omega }_1/\left(2\pi \right)$, $M_s$ is the structural mass, $f_1^{ref}$ and $M_s^{ref}$ are baseline values, and $w_f+w_m=1$. The weights are preference parameters rather than material constants; changing them changes the relative value assigned to dynamic performance and mass. A documentation audit performed during revision showed that the archived numerical output does not preserve the exact per-run values of $w_f$ and $w_m$. Rather than reconstructing unsupported metadata, the revised manuscript removes any claim that the archived response cases correspond to a specific unrecorded weight pair and reports an explicit a-posteriori sensitivity analysis for $w_f=$ 0.20, 0.30, 0.40, 0.50, and 0.80 in Section 6.8. This makes the weight dependence reproducible and transparent while avoiding retrospective fabrication.

The structural mass is computed from:

The optimization is subject to:

The ranges $0\le n$, $k\le 10$ span weak-to-steep power-law gradations, $0\le \alpha \le 0.5$ defines the admissible range of the analytical porosity model adopted here, and $5\le b/h\le 20$ spans moderately thick to slender cross-sectional configurations. These bounds define the design domain of the present study rather than universal material limits. For the archived response cases, the traceable values of $n$, $k$, and $\alpha$ are all interior to their prescribed intervals; their normalized distances to the nearest bound are reported in Section 6.7. The archived geometric entry associated with $b/h$ is not reproduced in the revised response table because the surviving historical tabulation is inconsistent with the stated $5\le b/h\le 20$ design domain, and no corrected value is inferred. In addition, the archived $n=k$ trade-off sweep is re-evaluated over reduced domains ($0\le n=k\le 8$ and $0\le n=k\le 10$) to show explicitly when a mass-dominant scalarization becomes bound-active.

The axial normal stress associated with the HSDT strain field is evaluated as:

and the maximum stress is defined as:

where the maximum is evaluated over the longitudinal stations and cross-sectional quadrature points used in the numerical implementation. The generalized eigenproblem is solved with LAPACK DSYGV using the standard type-1 form $\mathbf{K}\mathbf{\varphi }={\omega }^2\mathbf{M}\mathbf{\varphi }$. For this problem type, the returned eigenvectors are mass-normalized so that ${\mathbf{\varphi }}_i^T\mathbf{M}{\mathbf{\varphi }}_j={\delta }_{ij}$ [81]. The first mass-normalized eigenvector is therefore used consistently to reconstruct ${\epsilon }_x$ and ${\sigma }_x$, removing the otherwise arbitrary modal amplitude from the numerical definition. The stress check is expressed through the normalized ratio $\widehat{\sigma }={\sigma }_x/{\sigma }_{allow}$. The admissibility constraint is therefore:

This formulation makes the modal stress indicator unique and reproducible for the adopted eigenproblem normalization. It remains a modal screening criterion, not a substitute for a service-load strength assessment under a prescribed external load. The normalized stress field is recomputed for every trial design before acceptance.

4.2. BFGS solution strategy

The BFGS algorithm approximates the Hessian of the objective without explicitly computing second-order derivatives. Unlike population-based searches previously applied to FG frequency optimization [62, 63, 65], the present quasi-Newton procedure uses local gradient and curvature information while enforcing feasibility. The positive-definite rank-two update follows the classical variable-metric developments of Fletcher and Goldfarb [70, 71]. At iteration $j$, the search direction ${\mathbf{p}}_j$ is obtained from:

where ${\mathbf{B}}_j$ is the Hessian approximation.

The objective gradient is evaluated by a central finite-difference approximation:

with the nominal perturbation $\mathrm{\Delta }{x}_i=$ 10^-4. Finite differences can become expensive for high-dimensional problems; however, the present design space contains only four variables. A central-difference gradient therefore requires eight perturbed objective evaluations per iteration, in addition to the unperturbed evaluation. Central differences are retained because their leading truncation error is smaller than that of a one-sided forward approximation. The archived BFGS trace in Fig. 10 is retained only as a descent diagnostic and is not used to certify strict first-order convergence. A separate a-posteriori perturbation audit is reported in Section 6.9. It applies the same central-difference formula at three locations of a shape-preserving interpolation of the archived normalized frequency response ($n=2$, 5, and 8) and demonstrates systematic convergence of the derivative estimate as the perturbation is reduced from 10^-1 to 10^-6. This isolates finite-difference accuracy from missing historical optimizer metadata and provides a reproducible sensitivity check for the adopted 10^-4 perturbation.

The trial design is:

where the step length ${\eta }_j$ is determined using Wolfe conditions. The sufficient-decrease and curvature parameters are $c_1=\mathrm{ }$10^-4 and $c_2=$ 0.9, respectively. After each trial step, the design vector is projected onto the admissible bound domain and the stress condition $g_{\sigma }\le 0$ is checked before acceptance. The Hessian approximation is updated as:

with:

The iteration is stopped when:

or when the maximum number of iterations is reached.

Algorithm 1. Constrained BFGS optimization coupled with the HSDT dynamic model

Step 1. Define the initial design vector ${\mathbf{x}}_0={\left[n_0,k_0,{\alpha }_0,{\left(b/h\right)}_0\right]}^T$.

Step 2. Set ${\mathbf{B}}_0=\mathbf{I}$, choose ${ϵ}_g$, ${ϵ}_J$, and $N_{\mathrm{m}\mathrm{a}\mathrm{x}}$.

Step 3. Compute $E\left(y,z;{\mathbf{x}}_j\right)$ and $\rho \left(y,z;{\mathbf{x}}_j\right)$.

Step 4. Evaluate the stiffness and mass coefficients by cross-sectional integration.

Step 5. Assemble $\mathbf{K}\left({\mathbf{x}}_j\right)$ and $\mathbf{M}\left({\mathbf{x}}_j\right)$.

Step 6. Solve $\mathrm{d}\mathrm{e}\mathrm{t}\left[\mathbf{K}\left({\mathbf{x}}_j\right)-{\omega }^2\mathbf{M}\left({\mathbf{x}}_j\right)\right]=0$.

Step 7. Compute $f_1\left({\mathbf{x}}_j\right)$, $M_s\left({\mathbf{x}}_j\right)$, $J\left({\mathbf{x}}_j\right)$, and $g_{\sigma }\left({\mathbf{x}}_j\right)$.

Step 8. Approximate $\nabla J\left({\mathbf{x}}_j\right)$ using central finite differences, Eq. (32).

Step 9. Solve ${\mathbf{B}}_j{\mathbf{p}}_j=-\nabla J\left({\mathbf{x}}_j\right)$.

Step 10. Perform a Wolfe-condition line search ($c_1=$10^-4, $c_2=$ 0.9) and project the design vector onto the bound domain.

Step 11. Accept ${\mathbf{x}}_{j+1}$ only if the bound constraints and stress admissibility condition are satisfied.

Step 12. Update ${\mathbf{B}}_{j+1}$ using Eq. (34).

Step 13. Set $j\leftarrow j+1$ and repeat Steps 3-12 until the stopping criteria in Eq. (36) are satisfied or $N_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is reached; then stop.

6.1. Effect of porosity coefficient on natural frequencies

The porosity coefficient directly affects the stiffness and inertia distributions through Eq. (4-5). Increasing $\alpha$ generally reduces the effective stiffness of the beam, which leads to lower natural frequencies when the stiffness reduction dominates the mass reduction. This trend agrees with previous porous-beam studies, while its magnitude remains sensitive to the assumed porosity pattern, elastic foundation, boundary conditions, and shear-deformation model [12-16, 18-20, 35]. Fig. 3 illustrates the variation of the natural frequency with the porosity coefficient for selected values of $n$ and $k$.

Fig. 3Variation of the natural frequency with the porosity coefficient for selected bidirectional power-law indices

6.2. Effect of wave number and porosity distribution

The wave number influences the modal dispersion characteristics of the beam. For the same design vector, even and uneven porosity distributions lead to different stiffness and inertia distributions over the cross-section. Comparable dispersion sensitivity has been reported for porous FG beams and imperfect FG plates under higher-order kinematics [27-30, 36]. Fig. 4 compares the wave-number-dependent frequency response for even and uneven porosity models.

Fig. 4Natural frequency variation with wave number for even and uneven porosity distributions

6.3. Effect of power-law indices

The power-law indices $n$ and $k$ control the ceramic-rich and metal-rich regions along the thickness and width directions. Their influence on the first two natural frequencies is shown in Figs. 5(a) and 5(b). Earlier one-directional and multidirectional studies likewise demonstrate that power-law indices interact strongly with porosity and thermal or geometric parameters [11, 19, 30, 41-44, 51, 53, 79, 80]. These results should therefore be interpreted together with the porosity coefficient because the optimized stiffness-to-weight response depends on the combined effect of gradation and porosity.

The frequency decrease with increasing power index should be interpreted together with the associated mass change. Fig. 6 shows the corresponding weight-reduction trend, while Fig. 7 compares the normalized modal trends. These complementary plots emphasize that the optimization is governed by a coupled stiffness-inertia balance rather than by frequency alone.

Fig. 5Natural frequencies versus the power-law index: a) first natural frequency, b) second natural frequency

a)

b)

Fig. 6Weight-reduction trend versus the power-law index for the investigated parameter sweep

Fig. 7Normalized modal-frequency trends versus the power-law index

6.4. Phase velocity response

The phase velocity is computed from Eq. (23). Its variation with the power-law index reflects the combined influence of material gradation, porosity, wave number, scale effects, and, in broader service conditions, thermal fields [11, 27, 30, 36, 83]. Fig. 8 presents the phase velocity response for selected porosity coefficients.

Fig. 8Phase velocity variation with the power-law index for selected porosity coefficients

6.5. Stress admissibility

The stress constraint is evaluated using the dimensionless ratio in Eq. (30) and the mass-normalized first eigenvector defined in Section 4.1. A trial design passes the modal screening step when $g_{\sigma }\le 0$. Including such a criterion is consistent with failure-aware optimization strategies in composite mechanics, where frequency or stiffness improvement alone is insufficient to establish design admissibility [69]. Fig. 9 illustrates the dimensionless modal-stress ratio used in this screening. Because the quantity is modal and mass-normalized, it is interpreted as a reproducible comparative indicator rather than as a service-load strength prediction.

Fig. 9Dimensionless normalized modal-stress distribution σ^=σx/σallow used for the admissibility check

6.6. Optimization convergence

Fig. 10 reproduces the archived ten-iteration BFGS descent trace. The gradient norm decreases monotonically from 0.80 to 0.05, but the surviving log does not demonstrate the strict 10^-4 first-order tolerance stated in the generic stopping rule. The revision therefore no longer uses this trace as proof of strict convergence or of a unique mathematical optimum. It is retained only as an algorithmic descent diagnostic, while Eq. (36) states the stopping rule required for a fully documented new run. The recorded gradient-norm values are listed in Table 2, and the archived performance cases retained from the original numerical campaign are summarized in Table 3.

Fig. 10Archived BFGS descent trace. The curve is reported as a monotonic descent diagnostic and is not used to claim strict first-order convergence

6.7. Audit of prescribed design bounds

The rationale for the bounds is discussed in Section 4.1; here their influence is checked against the traceable archived response cases. For a scalar variable $x\in \left[x_{min},x_{max}\right]$, the normalized distance to the nearest bound is defined as:

A value of zero indicates an active box constraint, while larger values indicate an interior point. Table 4 shows that the traceable $n$, $k$, and $\alpha$ values are not located on their bounds. The historical $b/h$ entries are deliberately excluded for the provenance reason stated in Table 3. The archived one-dimensional $n=k$ sweep is also re-evaluated over $0\le n=k\le 8$ and $0\le n=k\le 10$ in Section 6.8; that comparison shows explicitly that a mass-dominant preference can become bound-active whereas a frequency-dominant preference remains at the high-frequency end of the explored response set.

6.8. A-posteriori sensitivity to objective weights

To address the dependence of the scalar objective on $w_f$ and $w_m$ without assigning unverifiable historical weights, the archived parameter sweep in Figs. 6-7 was reprocessed. The first normalized frequency branch ${\bar{f}}_1$ and the normalized mass ratio $\bar{M}=1-R_W/100$ were combined as:

where $R_W$ is the archived weight-reduction percentage. Table 5 summarizes representative preferences. The preferred point changes from the high-index, low-mass end of the archived sweep for mass-dominant weighting to the high-frequency end for frequency-dominant weighting. Thus, the preferred trade-off changes materially with the scalarization weights, and the revision does not claim a weight-independent optimum. This directly demonstrates why the weighting factors must be reported whenever a specific scalarized optimum is claimed. The same preference sensitivity is visualized in Fig. 11.

Fig. 11Preference sensitivity of the archived n=k frequency-mass sweep. The plot is an a-posteriori scalarization analysis of archived response data, not a retrospective reconstruction of missing BFGS weight metadata

The same data also provide a direct design-domain check. For mass-dominant weights the preferred point moves toward the largest available power index, so reducing the domain from $0\le n=k\le 10$ to $0\le n=k\le 8$ changes the selected endpoint. For frequency-dominant weights the preferred point remains at $n=k=0$. This confirms that bound sensitivity depends on the engineering preference and is explicitly reported rather than hidden inside the optimization setup.

6.9. Finite-difference perturbation consistency

A perturbation-size consistency check was performed at three locations of the archived normalized first-mode response, $n=2$, 5, and 8. A shape-preserving piecewise-cubic interpolant was constructed from the archived points in Fig. 7, and the derivative was evaluated with the same central-difference expression as Eq. (32). Table 6 reports the maximum relative difference, over the three audit points, between the central-difference estimate and the derivative of the interpolant. The error decreases systematically as the perturbation is reduced. At $\mathrm{\Delta }n=$10^-4 the maximum relative difference is approximately 1.33×10^-3%, providing direct numerical support for the adopted perturbation on the smooth archived response. The reduction of the derivative discrepancy with decreasing perturbation is shown in Fig. 12.

Fig. 12Finite-difference consistency at three archived response points. The monotonic reduction of the derivative discrepancy supports the use of Δxi= 10-4 for the smooth response audit

6.10. Robustness across porosity and boundary conditions

The direct verification in Table 1 is intentionally restricted to case-matched numerical comparisons. To broaden the robustness assessment without mixing non-equivalent models, Table 7 adds an independent trend-level comparison with studies that explicitly vary porosity pattern and boundary condition. The present trends are consistent with this literature: porosity distribution changes stiffness and inertia differently, boundary conditions materially shift the spectrum, and two-directional gradation interacts with porosity in determining the fundamental frequencies. This expanded comparison addresses robustness across modeling scenarios while preserving the distinction between quantitative verification and qualitative corroboration.

6.11. Physical interpretation of the coupled dynamic changes

The simultaneous changes in frequency, mass, and phase velocity originate from the fact that stiffness and inertia do not weight the spatial material distribution in the same way. Structural mass is governed primarily by the cross-sectional integral of density, whereas bending-related stiffness coefficients contain spatial weighting terms such as $z^2$. Consequently, changing $n$ and $k$ can redistribute the ceramic-rich phase toward mechanically influential regions and modify $D_{11}$ and the higher-order stiffness coefficients more strongly than the total mass.

In reduced form, the first natural frequency scales approximately as:

so a mass reduction does not by itself guarantee a frequency increase. Porosity introduces a second coupled mechanism because $\alpha$ modifies both $E\left(y,z\right)$ and $\rho \left(y,z\right)$. Even and uneven porosity distributions also place stiffness and inertia reductions in different regions of the cross-section. The optimizer therefore acts on the combined spatial redistribution of stiffness and inertia rather than on a single monotonic material parameter. The same redistribution enters the dispersion relation and explains the sensitivity of phase velocity to gradation and porosity.

This interpretation is consistent with the broader structural-dynamics literature in which stiffness-inertia tailoring is used to control vibration transmission and wave propagation. Recent lattice- and TPMS-based metastructure studies further demonstrate how spatial structural tailoring can affect acoustic transmission and vibration control [73-78]. The present beam formulation does not model fluid–structure interaction or acoustic radiation, but it provides a mechanics-based starting point for such extensions.

6.12. Comparison with existing studies

The present candidate response is placed in context with representative reference studies in Table 8. Because the listed models do not share all assumptions, this table is explicitly contextual and is not used as a validation table. Quantitative verification remains confined to Table 1.

The value 150.2 Hz for performance case 1 is approximately 24.6 % higher than the 120.5 Hz value listed for Benadouda et al. [27] and 15.5 % higher than the 130.0 Hz value listed for Batou et al. [28]. These percentages are reported as contextual performance differences, not as validation errors, because the reference formulations and structural models are not identical. Formal verification is restricted to case-matched comparisons such as those reported in Table 1.

6.13. Practical limitations and implications

The present framework assumes rule-of-mixtures effective properties, prescribed analytical porosity distributions, linear-elastic Touratier-HSDT beam kinematics, deterministic material parameters, and no explicit temperature-, moisture-, or damage-dependent constitutive evolution. These assumptions have direct practical implications. Temperature and moisture can modify constituent properties and therefore shift both the dynamic response and the design optimum [11, 17, 84]. Manufacturing scatter in porosity can move a deterministic optimum away from its nominal performance, motivating uncertainty quantification and robust optimization. The beam-level idealization also cannot resolve local three-dimensional stress concentrations, joints, detailed pore morphology, or manufacturing defects.

Finite-difference gradients are another numerical limitation: their cost grows with the number of design variables and their accuracy depends on the perturbation size. In the present four-variable problem this cost remains moderate, and the a-posteriori consistency audit in Section 6.9 shows stable differentiation of the archived smooth response at the nominal perturbation. Likewise, the scalar objective depends on the selected weights $w_f$ and $w_m$; the explicit sweep in Section 6.8 demonstrates that the preferred trade-off can change with the engineering preference. A Pareto formulation remains a natural future extension. Classical, multilayer, and $h$-$p$ finite-element formulations provide established routes for extending the beam model to larger structural discretizations [52, 54, 85-88].

Accordingly, the archived performance cases should be interpreted as design responses produced within the stated deterministic beam framework, not as proof of a unique global optimum or as universally optimal manufactured components. Future work should address thermo-hygro-mechanical properties, porosity uncertainty, Pareto/robust optimization, damage sensitivity, high-fidelity finite-element analysis, and experimental validation.