© K.W. Christophe et al., Published by EDP Sciences, 2026

1 Introduction

Glass-fiber-reinforced polymer (GFRP) composites are cornerstone materials for high-performance applications in aerospace, automotive, and space engineering, where mass reduction directly enhances energy efficiency, structural integrity, and lifecycle performance. Their macroscopic behavior is dictated by a complex heterogeneous microstructure fiber bundles, statistical orientation distributions, interfacial zones, and manufacturing-induced defects such as weave irregularities, porosity, and cure non-uniformity which generates strongly coupled thermo-chemo-mechanical fields [1–3]. Modeling and predicting these interactions under realistic processing and service conditions remain major challenges due to their intrinsic multiscale and multiphysics complexity.

Classical computational approaches, particularly multiscale finite-element (FE) formulations, have advanced the mechanistic understanding of continuous fiber composites [4]. However[CE1], these methods exhibit well-recognized limitations when addressing: (i) fully coupled physics across hierarchical scales; (ii) propagation, quantification, and calibration of manufacturing uncertainties; and (iii) adaptive updating of material and process parameters during curing or in-service evolution. Recent developments in physics-informed machine learning (PIML) offer promising avenues to overcome these limitations [5]. By embedding governing equations or operator constraints within learning architectures, PIML enables data-driven correction and online recalibration while preserving physical consistency. Representative recent works have demonstrated significant progress: physics-informed neural networks (PINNs) for thermochemical cure modeling and thermal transport in composites; hybrid FE ML surrogates for accelerated prediction of resin flow and heat transfer; and real-time defect detection algorithms for porosity and cure non-uniformity during manufacturing. Nevertheless, these studies typically treat isolated phenomena e.g., single-scale heat conduction, uncoupled resin flow, or local damage detection or rely on empirical, non-physics-regularized surrogates [2,6,7]. As a result, a unified framework integrating multiscale physics, multiphysics coupling, data-driven adaptivity, and uncertainty-aware optimization is still lacking.

To address these gaps, we develop a mathematically grounded hybrid modeling framework that tightly couples: (1) a rigorous multiscale and multiphysics partial differential equation (PDE) model describing thermo-chemo-mechanical interactions, resin flow dynamics, porosity evolution, and residual stress formation; (2) a physics-regularized machine-learning operator acting as an adaptive corrective term to compensate modeling simplifications or parameter drift under uncertain conditions; (3) an adjoint-based calibration and uncertainty quantification pipeline combining Sobol mutual-information sensitivity analysis with surrogate-assisted Monte Carlo sampling for robust multi-objective process optimization; and (4) systematic integration pathways toward an adaptive digital twin for composite manufacturing. The resulting hybrid PDE PIML estimator admits provable properties including existence, uniqueness, and explicit bounds on model data discrepancy ensuring mathematical robustness and physical consistency. It improves predictive fidelity under limited or noisy in-situ measurements, while enabling Pareto-optimal trade-offs between defect mitigation, mechanical performance, and production cost.

The remainder of the paper is organized as follows. Section 2 presents the governing equations, multiscale constitutive modeling, and the associated numerical discretization strategy. Section 3 details the experimental protocol, material characterization procedures, and in-situ data acquisition. Section 4 describes the hybrid PDE PIML architecture, including operator regularization, training, and online updating. Section 5 reports the optimization results, sensitivity analyses, and uncertainty quantification outcomes. Finally, Section 6 summarizes the main conclusions and outlines perspectives for real-time digital twins and extensions to other composite systems.

2 Theoretical and methodological framework

2.1 Multiscale modeling and integration of physical scales

Accurately predicting the mechanical behavior of glass fiber cable composites requires a hierarchical multiscale modeling framework that links microstructural phenomena to the global structural response. In this study, we adopt a rigorous multiscale methodology founded on homogenization theory and micromechanics of heterogeneous materials, ensuring seamless integration of the micro, meso, and macro scales.

(a) Microscale: Representative volume element (RVE)

At the microscale, individual fibers, fiber bundles, and local heterogeneities such as porosity and weaving defects are modeled through a Representative Volume Element (RVE). The local constitutive behavior is described by:

$${\text{σ}}_{\text{ij}}\left(\text{x}\right)={\text{C}}_{\text{ijkl}}\left(\text{x}\right){\text{ε}}_{\text{kl}}\left(\text{x}\right)$$

where σ_ij and ε_kl are the stress and strain tensors, and C_ijkl(x) is the local stiffness tensor dependent on the microstructure. Analytical approximations based on Eshelby’s inclusion theory are utilized to estimate local field perturbations caused by fiber inclusions within the matrix, providing precise quantification of heterogeneity effects [2,5].

(b) Micro → Meso transition: homogenization

To upscale the microscale response, homogenization techniques compute effective constitutive tensors C^* representing averaged material behavior:

• Numerical homogenization involves finite element analysis of the RVE under periodic boundary conditions:

  $$\text{C}{\text{}}_{\text{ijkl}}^{\text{*}}=\frac{1}{\left|\text{V}\right|}\underset{\text{V}}{\overset{ }{{\displaystyle \int }}}\text{}{\text{C}}_{\text{ijkl}}\left(\text{x}\right)\text{dV}$$

  accounting for anisotropy, porosity, and stochastic fiber orientation distributions.

• Analytical homogenization employs micromechanical models such as Mori-Tanaka or self-consistent schemes to derive closed-form relations dependent on fiber volume fraction and distribution [6].

(c) Mesoscale: Yarns and cable architecture

At the mesoscale, woven yarns and cable bundles are modeled as periodic inclusions embedded in an effective matrix. Their constitutive behavior is expressed as:

$${\text{σ}}^{\text{meso}}={\text{C}}^{\left(\text{*},\text{meso}\right)}:{\epsilon }^{\text{meso}}$$

where C^(*,meso) integrates microscale homogenized properties and geometric imperfections from textile architecture.

(d) Meso → Macro transition: composite structures

Macroscale finite element simulations incorporate the effective mesoscale properties. Coupling between scales is ensured by:

• Upscaling effective tensors C^* to inform macroscopic constitutive laws;

• Downscaling macroscopic stresses and strains back to microscale and mesoscale fields, allowing accurate damage and residual stress evaluation.

The overall macroscopic constitutive law is formulated as: ${\text{σ}}^{\text{macro}}={\text{C}}^{\text{*},\text{macro}}:{\epsilon }^{\text{macro}}+\text{R}\left(\text{x},\text{ t}\right)$,

where R(x, t) represents residual effects arising from processing, including infusion, curing, and thermal shrinkage.

Innovation and Significance

The proposed hierarchical multiscale approach sets itself apart from recent studies (2023–2025) by concurrently addressing key challenges [7]: realistic modeling of manufacturing defects and microstructural heterogeneities across all scales; full integration of coupled thermo-hygromechanical and residual stress effects within a unified multiscale framework; and implementation of adaptive cross-scale feedback mechanisms that enable dynamic recalibration to optimize design and manufacturing processes. Together, these innovations establish a robust, predictive, and calibratable framework tailored for the design and optimization of glass fiber composites, specifically targeting aerospace, automotive, and space industries where structural reliability, lightweight performance, and process efficiency are paramount.

2.2 Multiphysics modeling

The proposed multiphysics modeling framework captures the strongly coupled thermal, mechanical, fluidic, and chemical phenomena governing the effective behavior of glass fiber cable composites [8]. This formulation enables a unified description of manufacturing processes such as resin infusion and curing, which simultaneously influence the structural integrity, residual stress distribution, and long-term performance of the composite.

The stress tensor σ is expressed through a coupled constitutive relation:

$$\text{σ}={\text{C}}^{\text{eff}}:\left(\text{ε}-{\text{ε}}^{\text{th}}-{\text{ε}}^{\text{res}}\right)$$

where ε denotes the total strain tensor, ${\epsilon }^{\text{th}}=\text{α}\hspace{0.17em}\text{ΔT}$ representing thermal strains (with α being the thermal expansion tensor) and ${\epsilon }^{\text{res}}={\text{f}}_{\text{proc}}\left(\text{t}\right)$ accounts for the residual strain induced by processing effects such as resin shrinkage and thermal gradients.

The temperature field satisfies the transient heat conduction equation with a chemical source term:

$${\rho c}_{\text{p}}\text{}\frac{\partial \text{T}}{\partial \text{t}}\text{}-\nabla \cdot \left(\text{k}\nabla \text{T}\right)={\text{Q}}_{\text{chim}}\left(\text{t}\right)$$

where ρ is the material density, c_p the specific heat capacity, k the thermal conductivity, and Q_chim(t) the heat generation rate associated with exothermic resin polymerization.

The infusion process is modeled using a generalized Darcy’s law coupled with mass conservation:

$${\text{v}}_{\text{f}}\text{}=-\frac{\text{K}\left(\varphi \right)\text{}}{\text{μ}}\nabla \text{p},\hspace{0.17em}\frac{\partial \varphi \text{}}{\partial \text{t}}+\nabla \cdot \left(\varphi {\text{v}}_{\text{f}}\text{}\right)=0.$$

where v_f represents the resin infiltration velocity, K(ϕ) the porosity-dependent permeability, μ the resin viscosity, p the pressure field, and ϕ the local porosity C^eff(t, ϕ), ensuring realistic coupling between process physics and material response.

The resin curing kinetics follows an Arrhenius-type thermochemical law [9]:

$$\frac{\text{dα}}{\text{dt}}\text{}=\text{k}{\text{}}_0\text{exp}\left(-\frac{\text{E}{}_{\text{a}}}{\text{RT}}\right){\left(1-\text{α}\right)}^{\text{n}},$$

where α(t) is the degree of cure, k₀ the pre-exponential factor, E_a the activation energy, R the universal gas constant, and n the reaction order. The thermo-chemical coupling is expressed as:

$$\text{Q}{}_{\text{chim}}\left(\text{t}\right)={\rho }_{\text{r}}{\text{H}}_{\text{r}}\frac{\text{dα}}{\text{dt}},$$

with ρ_r the resin density and H_r the reaction heat per unit mass.

The global behavior is obtained by solving the fully coupled multiphysics system:

$$\{\begin{array}{c}{\rho \text{c}}_{\text{p}}\frac{\partial \text{T}}{\partial \text{t}}-\nabla \cdot \left(\text{k}\nabla \text{T}\right)={\rho }_{\text{r}}{\text{H}}_{\text{r}}\frac{\text{dα}}{\text{dt}}\\ \frac{\text{dα}}{\text{dt}}=\text{k}{}_0\text{exp}\left(-\frac{\text{E}{}_{\text{a}}}{\text{RT}}\right){\left(1-\text{α}\right)}^{\text{n}}\\ \begin{array}{l}{\text{v}}_{\text{f}}=-\frac{\text{K}\left(\varphi \right)}{\text{μ}}\nabla \text{p},\frac{\partial \varphi }{\partial \text{t}}+\nabla \cdot \left(\varphi {\text{v}}_{\text{f}}\right)=0\\ \text{σ}={\text{C}}^{\text{eff}}\left(\text{t},\varphi ,\text{ α}\right):\left(\text{ε}-\text{α}\hspace{0.17em}\text{ΔT}-{\text{ε}}^{\text{res}}\left(\text{t}\right)\right)\end{array}\end{array}$$(E)

This strongly coupled set of nonlinear equations is solved iteratively within a monolithic finite element framework, ensuring convergence and physical consistency across all coupled phenomena

To improve clarity and accessibility for readers from diverse backgrounds, Figures 1 and 2 illustrate the multiscale and multiphysics interactions occurring during the processing of glass-fiber composites. Figure 1 presents numerical simulation results, while Figure 2 offers a schematic overview of the fundamental physical mechanisms at different scales. Together, these figures provide both quantitative data and conceptual understanding of the coupled thermo-chemo-mechanical behavior discussed.

Figure 1 displays simulated 3D fields of temperature, displacement, and cure degree within the composite throughout the curing process. The temperature field reveals localized heat buildup in inner layers, resulting in thermal gradients that cause stress redistribution visible in the displacement patterns. Variations in the degree of cure highlight resin flow restrictions and porosity effects, confirming the strong coupling between thermal, chemical, and mechanical phenomena that control stiffness changes and residual strain development.

Figure 2 summarizes the hierarchical coupling mechanisms spanning micro-, meso-, and macro-scales. This schematic illustrates feedback loops linking resin infiltration, heat generation, curing kinetics, and deformation. It connects microstructural fiber-matrix interactions with overall structural responses, offering a comprehensive view of how localized manufacturing processes influence the global mechanical behavior of the composite.

Fig. 1 3D visualization of multiscale multiphysics couplings in glass-fiber composites.

Fig. 2 Schematic illustration of the multiscale and multiphysics coupling mechanisms, emphasizing interactions among thermal, mechanical, chemical, and fluid phenomena across micro-, meso-, and macro-scales during curing and in-service conditions.

2.3 Data-driven hybrid modeling

In advanced composite manufacturing, deterministic multiphysics models alone cannot guarantee engineering-grade predictive reliability due to process variability (e.g., fiber volume fraction fluctuations, cure-dependent viscosity, non-uniform permeability), sensor noise, and operational disturbances inherent to industrial environments. Real-time measurements from thermocouples, dielectric analysis, fiber Bragg gratings, and in-mold pressure sensors frequently deviate from simulation outputs, especially under strongly nonlinear thermo-chemo-mechanical regimes. To address this model plant mismatch, we introduce a precise hybrid modeling strategy that integrates first-principles physics with data-driven correction terms.

At the core of our formulation, Theorem 1 establishes the mathematical well-posedness (existence, uniqueness, and stability) of the hybrid PDE PIML estimator, demonstrating that the combined model remains physically consistent and numerically convergent under realistic industrial regularity assumptions. These theoretical guarantees are essential for applied engineering: they ensure that the hybrid estimator behaves as a stable digital proxy of the manufacturing process, that online updates remain bounded under sensor noise, and that optimization algorithms relying on adjoint gradients remain robust and scalable.

Building on these foundations, we exploit Theorem 1 to construct a hybrid surrogate capable of dynamically adjusting model parameters in response to real-time measurements, thereby compensating for modeling simplifications (e.g., homogenized porous media flow, reduced-order cure kinetics) or uncertain process inputs. This mathematically certified hybrid architecture enables reliable prediction across the full curing and consolidation cycle, while supporting multi-objective optimization and uncertainty-aware decision-making in automated composite processing.

Theorem 1 (Hybrid Modeling for Multiphysics Systems):

Let ${\text{u}}_{\text{phys}}\left(\text{x},\text{t};\text{θ}\right)$ be a deterministic multiphysics model on ${\text{Ω}}_{\text{T}}=\text{Ω}\times \left[0,{\text{T}}_{\text{f}}\right]$, parameterized by $\text{θ}\in \text{Θ}\subset {ℝ}^{\text{p}}$, and ${\text{u}}_{\text{data}}\left(\text{x},\text{t};\text{β}\right)$ a statistical model trained on experimental data $\mathcal{D}={\left\{\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}},{\text{y}}_{\text{i}}\right)\right\}}_{\text{i}=1}^{\text{N}}$, parameterized by $\text{β}\in \text{B}\subset {ℝ}^{\text{q}}$. Define the hybrid model:

$${\text{u}}_{\text{λ},\text{θ},\text{β}}\left(\text{x},\text{t}\right)=\left(1-\text{λ}\right){\text{u}}_{\text{phys}}\left(\text{x},\text{t};\text{θ}\right)+{\text{λ}\theta }_{\text{data}}\left(\text{x},\text{t};\text{β}\right),\forall \text{ λ}\in \left[0,1\right].$$

The regularized cost functional is:

$$\text{J}\left(\text{λ},\text{θ},\text{β}\right)=\frac{1}{2}{\displaystyle \sum }_{\text{i}=1}^{\text{N}}{||}_{\text{u}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}{}^{||}+\frac{\text{γ}}{2}\text{R}\left(\text{θ},\text{β}\right)+\frac{\text{μ}}{2}||\mathcal{F}{\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\right)}_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}^2$$

where F is the differential operator of the multiphysics system, R a strictly convex regularization, and γ, μ>0.

Under the following regularity and compactness assumptions:

(H1) : ${\text{u}}_{\text{phys}}\left(\cdot ;\text{θ}\right)$ is continuous and differentiable with respect to θ in the functional space $\text{X}={\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)$ (or H¹ depending on F), and the mapping $\text{θ}\mapsto {\text{u}}_{\text{phys}}\left(\cdot ;\text{θ}\right)$ is locally Lipschitz;

(H2) : $\text{β}\mapsto {\text{u}}_{\text{data}}\left(\cdot ;\text{β}\right)$ is C¹ in B and the regularization R is strictly convex and coercive;

(H3) : $\mathcal{F}:\text{X}\to {\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)$ is continuously differentiable C¹.

Hybrid Model Properties and Practical Optimization Strategy

Given the assumptions (H1)-(H3) and the composite formulation of the hybrid estimator, Theorem 1 ensures that the regularized cost functional $\text{J}\left(\text{λ},\text{θ},\text{β}\right)$ admits at least one minimizer $\left({\text{λ}}^{\text{*}},{\text{θ}}^{\text{*}},{\text{β}}^{\text{*}}\right)\in \left[0,1\right]\times \text{Θ}\times \text{B}\text{.}$. In engineering terms, this guarantees that the hybrid PDE PIML formulation remains well-posed, meaning that the calibration problem always yields a feasible and physically consistent solution. Moreover, when the regularization R enforces sufficient convexity an assumption typically satisfied in composite processing models where parameters enter linearly the minimizer becomes unique. This result is essential for industrial deployment, as it prevents ambiguity during real-time parameter updating and ensures stable digital-twin operation.

Practical Determination of the Blending Coefficient λ

Although introduced as a theoretical convex weight, the blending parameter λ plays a direct engineering role by controlling the contribution of physics-based predictions versus data-driven corrections. We determined its optimal value λ^* through a structured line search combined with sensitivity analysis. The parameter was swept over the interval [0,1] using a resolution of $\text{Δλ}=0.05$, and the corresponding validation error $\text{J}\left(\text{λ}\right)= \underset{\text{θ},\text{β}}{\min }\text{J}\left(\text{λ},\text{θ},\text{β}\right)$ was evaluated at each step. The optimum ${\text{λ}}^{\text{*}}=0.62$ yielded the lowest residuals and the best agreement with in-situ measurements. The smooth convex profile of J(λ) around this value confirms the numerical stability of the hybrid estimator and validates the robustness of the blending mechanism.

Adjoint-Based Calibration and Computational Efficiency

To ensure scalability for large-scale composite curing simulations, we compute gradients of the objective functional using adjoint equations, which significantly reduces computational cost compared to finite differences. This enables efficient joint optimization of $\left(\text{λ},\text{θ},\text{β}\right)$ even when u_phys is governed by nonlinear, coupled thermo-chemo-mechanical PDEs. This adjoint-based strategy ensures that the framework remains deployable in industrial environments requiring near real-time decision support.

Engineering Impact

By integrating rigorously formulated PDE-based models with data-driven correction operators, the proposed hybrid framework achieves significantly enhanced predictive fidelity, particularly under manufacturing variability and uncertainty. This approach demonstrates robustness to measurement noise and imperfect boundary conditions, thereby ensuring reliable model outputs in realistic operational environments. The theoretical well-posedness established by Theorem 1 guarantees consistent and unique solutions during the calibration and optimization phases, which translates into stable and convergent performance of the estimation algorithm. Furthermore, the hybrid estimator’s smoothness and differentiability enable precise and reliable selection of process parameters within optimal manufacturing windows. Collectively, these attributes render the methodology a technically rigorous and practically effective computational tool for composite manufacturing, capable of accurately forecasting critical process variables such as porosity evolution, residual stress development, degree of cure progression, thermal gradients, and other key multiphysics phenomena governing final material quality and performance.

We have established, through Theorem 1, the existence of at least one minimizer $\left({\text{λ}}^{\text{*}},{\text{θ}}^{\text{*}},{\text{β}}^{\text{*}}\right)\in \left[0,1\right]\times \text{Θ}\times \text{B}$ for the regularized optimization problem associated with the proposed hybrid multiphysics-data-driven model. Under standard regularity assumptions and strict convexity of the cost functional particularly when the physical model is linear in its parameters and the regularization is sufficiently coercive this minimizer is furthermore unique. This result ensures the well-posedness and mathematical stability of the parameter estimation problem in the complex setting of composite materials engineering.

We have concisely formulated the stationary conditions of the optimization problem in Proposition 1. The associated adjoint formulation, critical for efficient gradient evaluation, enables scalable optimization while maintaining physical consistency of parameter updates. This adjoint-based approach significantly enhances convergence properties and robustness during the calibration of hybrid models, effectively reconciling first-principles physics with empirical data.

Corollary 1 provides explicit error and stability bounds for the approximate solution relative to the exact model, under standard assumptions of regularity and observability. This bound quantifies the model’s robustness to measurement noise and inherent data uncertainties, thereby strengthening confidence in predictive accuracy under realistic experimental conditions. For comprehensive mathematical proofs and detailed derivations, the reader is referred to Appendix A.

2.4 General methodology

We establish a unified hybrid variational adjoint framework that integrates multiscale material characterization, fully coupled multiphysics modeling, data-informed hybridization, adjoint-based sensitivity analysis, and large-scale optimization within a single mathematically consistent workflow. Unlike conventional sequential or surrogate-driven coupling strategies, our approach introduces a monolithic and differentiable formulation that maintains physical interpretability while enabling seamless algorithmic integration across scales. The workflow comprises five interdependent stages Characterization → Modeling → Simulation → Validation → Optimization as illustrated in Figure 3.

This figure illustrates the Unified Multiscale Multiphysics Hybrid Optimization Workflow, which integrates sequential and iterative stages encompassing experimental characterization, multiscale data acquisition, hybrid variational modeling, adjoint-based numerical simulation, robust validation, and objective-driven optimization.

• Experimental Characterization: We first acquire multiscale experimental fields temperature T, degree of cure α, porosity ϕ, and displacement u coupled with microstructural descriptors such as fiber orientation distributions, porosity topology, and defect statistics. These quantities serve as priors for the hybrid variational formulation and constrain the admissible space of multiscale parameters.

• Multiphysics and Hybrid Modeling: We construct a monolithic variational operator U ∈ X that couples thermal, chemical, and mechanical processes. The hybridization is introduced through a convex blending rule, $\text{U}\left(\text{λ},\text{θ},\text{β}\right)=\left(1-\text{λ}\right){\text{U}}_{\text{phys}}\left(\text{θ}\right)+{\text{λU}}_{\text{data}}\left(\text{β}\right)$ which embeds adaptive data-informed corrections directly into the variational structure. The blending parameter λ ∈ [0,1] is treated as an optimizable scalar, while θ and β denote physics and ML parameters, respectively. This formulation constitutes a core methodological contribution, as it ensures variational consistency, mathematical well-posedness, and smooth differentiability of the hybrid operator.

• Numerical Simulation and Adjoint Computation: The monolithic hybrid system is solved using nonlinear finite-element discretization. Adjoint equations associated with the cost functional $\text{J}\left(\text{λ},\text{θ},\text{β}\right)$, provide exact analytical gradients with respect to all parameters. This adjoint-based procedure ensures sensitivity-consistent parameter updates and allows efficient navigation of the high-dimensional parameter space inherent to strongly coupled multiphysics systems.

• Validation: Hybrid predictions are evaluated against experimental measurements using the discrepancy norm ${\Vert {\text{U}}_{\text{data}}-{\text{U}}_{\text{phys}}\Vert }_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}$ along with uncertainty-propagation metrics ${\text{σ}}_{\text{η}}$ derived from Monte-Carlo and Sobol-based sensitivity analyses. The optimal blending parameter λ^* is obtained through a line-search driven convex sensitivity procedure, which stabilizes the coupling between physics-based and data-driven contributions.

• Optimization: We determine the optimal parameters by minimizing the regularized objective functional

  $$\text{J}\left(\text{λ},\text{θ},\text{β}\right)=\frac{1}{2}{\displaystyle \sum }_{\text{i}=1}^{\text{N}}{\Vert {\text{u}}_{\text{λ},\text{θ},\text{β}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-\text{y}\Vert }_{}{}^2+\frac{\text{γ}}{2}\text{R}\left(\text{θ},\text{β}\right)+\frac{\text{μ}}{2}{\Vert \mathcal{F}\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\right)\Vert }_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}^2.$$

The theoretical results established in Theorem 1 guarantee existence, conditional uniqueness, and stability of the optimal solution. This ensures robust predictive capability under process variability, material heterogeneity, and sensor uncertainty.

The methodology we propose goes well beyond existing multiphysics ML combinations. We construct a precise, variationally consistent, and adjoint-differentiable hybrid framework that unifies multiscale characterization, multiphysics PDE modeling, data-informed corrections, and optimization. The resulting formulation preserves physical interpretability while offering a generalizable algorithmic structure suitable for reproducible, real-time, and uncertainty-aware digital twins for advanced composite manufacturing.

Fig. 3 Unified multiscale-multiphysics-hybrid optimization workflow.

3 Materials and processing

3.1 Description of the studied composite

We investigate a composite material composed of continuous glass fibers ${\left\{{\text{f}}_{\text{i}}\right\}}_{\text{i}=1}^{{\text{N}}_{\text{f}}}$ embedded in a thermosetting polymer matrix M. Each fiber f_i is characterized by its longitudinal modulus E_f, Poisson’s ratio ν_f, and density ρ_f. We define the matrix M by its thermo-mechanical properties Young’s modulus E_m, Poisson’s ratio ν_m, thermal expansion coefficient α_m, and specific heat capacity c_p,m along with its polymerization kinetics, represented by the degree of cure $\text{α}\left(\text{t}\right)\in \left[0,1\right]$, governed by an Arrhenius-type law [15]:

$\frac{\text{dα}}{\text{dt}}=\text{k}{}_0\text{exp}\left(-\frac{\text{E}{}_{\text{a}}}{\text{RT}}\right){\left(1-\text{α}\right)}^{\text{n}}$, where k₀ is the pre-exponential rate constant, E_a the activation energy, R the universal gas constant, and n the reaction order.

The effective fiber volume fraction is defined as $\varphi ={\text{V}}_{\text{f}}/{\text{V}}_{\text{tot}}$, with ${\text{V}}_{\text{f}}={\displaystyle \sum }_{\text{i}}\text{vol}\left({\text{f}}_{\text{i}}\right)$ and V_tot the total composite volume. Macroscopic effective properties $\left(\rho ,\text{C},{\text{C}}^{\text{eff}}\right)$ are obtained via multiscale homogenization, accounting for microstructural interactions, local porosity ϕ(x) and residual strains ε^res induced by manufacturing.

The experimental data acquisition system is founded upon an advanced instrumentation setup specifically designed to enable in situ measurement of critical physical parameters necessary for the accurate calibration and validation of hybrid thermo-chemo-mechanical models. Temperature fields within the composite during curing are captured through a distributed network of embedded micro-thermocouples, strategically positioned at multiple depths to resolve spatial and temporal thermal gradients with high fidelity. These discrete temperature measurements are complemented by high-resolution infrared thermography, employing an infrared camera sensitive to variations in radiative heat flux at the composite surface, thus enabling thermal mapping with sub-second temporal resolution.

The degree of cure, α(t), is quantified using differential scanning calorimetry (DSC) performed on representative resin specimens subjected to controlled thermal cycles. This approach is concurrently validated via embedded in situ dielectric analysis (DEA) sensors, which provide real-time monitoring of cure kinetics within the composite matrix. Mechanical displacement and strain fields are captured through digital image correlation (DIC) techniques applied to specimen surfaces, providing non-contact, high-precision measurements of deformation throughout thermo-mechanical cycling. To enhance measurement reliability, strain gauges are installed at critical stress concentration zones to validate and complement the DIC data.

We have developed a comprehensive experimental framework that ensures the acquisition of high-fidelity, coupled thermo-chemo-mechanical datasets essential for the rigorous calibration and validation of our hybrid modeling methodology. By synergistically integrating point-wise sensor measurements with full-field characterization techniques under representative processing conditions, we generate a robust and well-resolved database. This dataset forms the indispensable foundation supporting the development, verification, and implementation of advanced multiphysics numerical simulations and optimization algorithms explicitly tailored for composite manufacturing processes.

3.2 Analyzed manufacturing process

The composite is manufactured using a vacuum-assisted resin infusion process, adjustable according to industrial parameters $\text{P}=\left\{\text{T}\left(\text{t}\right),\text{p}\left(\text{t}\right),{\text{v}}_{\text{f}}\left(\text{t}\right),{\text{t}}_{\text{cure}}\right\}$, representing the temperature, applied pressure, resin infiltration velocity, and polymerization time, respectively [16]. The process is modeled by a coupled fluid-solid system on $\text{Ω}\times \left[0,{\text{T}}_{\text{f}}\right]$ :

$$\{\begin{array}{c}{\rho \psi }_{\text{p}}\frac{\partial \text{T}}{\partial \text{t}}-\nabla \cdot \left(\text{k}\nabla \text{T}\right)={\rho }_{\text{r}}{\text{H}}_{\text{r}}\frac{\text{dα}}{\text{dt}}\\ \frac{\text{dα}}{\text{dt}}=\text{k}{}_0\text{exp}\left(-\frac{\text{E}{}_{\text{a}}}{\text{RT}}\right){\left(1-\text{α}\right)}^{\text{n}}\\ {\text{v}}_{\text{f}}=-\frac{\text{K}\left(\varphi \right)}{\text{μ}}\nabla \text{p},\frac{\partial \varphi }{\partial \text{t}}+\nabla \cdot \left(\varphi {\text{v}}_{\text{f}}\right)=0\\ \text{σ}={\text{C}}^{\text{eff}}\left(\text{t},\varphi ,\text{ α}\right):\left(\text{ε}-\text{α}\hspace{0.17em}\text{ΔT}-{\text{ε}}^{\text{res}}\left(\text{t}\right)\right)\end{array}$$(E)

where each equation corresponds, respectively, to thermal conduction, resin polymerization kinetics, resin flow within the porous network, and residual mechanical behavior.

Critical Steps in the Industrial Process

1. Fiber preparation and placement f_i(x), determining the local fiber volume fraction ϕ(x) and microstructural heterogeneity.

2. Resin infusion with controlled velocity v_f(t) and porosity φ(x,t), minimizing infiltration defects.

3. Polymerization and thermal curing (T(t),α(t)), regulating residual stresses ε^res and final composite density.

4. Cooling and demolding, ensuring dimensional and mechanical stability.

This formalization explicitly links process parameters P, multiphysics fields $\left(\text{T},\text{ α}, \varphi ,\text{ u}\right)$ and the macroscopic properties of the composite, providing a rigorous basis for simulation, validation, and optimization.

Critical Process Parameters

$$\text{P}=\left\{\text{T}\left(\text{t}\right),\hspace{0.17em}\text{p}\left(\text{t}\right),\hspace{0.17em}{\varphi }_{\text{r}}\left(\text{t}\right),\hspace{0.17em}{\text{v}}_{\text{f}}\left(\text{t}\right),\hspace{0.17em}\varphi \left(\text{x},\text{t}\right)\right\}$$

where T(t) is the polymerization temperature, p(t) the applied pressure, φ_r(t) the resin volume fraction, v_f(t) the infiltration velocity, and φ(x,t) the local porosity.

Precise monitoring of P is essential to minimize defects, control residual stresses, and ensure macroscopic homogeneity, providing a solid foundation for reliable simulation and optimization.

3.3 Experimental setup and characterization

We designed a comprehensive experimental campaign to acquire high-fidelity, multiscale datasets essential for the rigorous calibration and validation of our hybrid multiphysics modeling framework, replicating industrial composite manufacturing conditions. The materials, characterization methods, and data acquisition protocols are detailed as follows.

Materials: The studied composite system consists of continuous ECR-type E-glass fibers (Owens Corning, product code EC45-2400, 15μm diameter, tensile modulus 73GPa, with epoxy-compatible sizing) embedded within a high-performance thermosetting epoxy matrix (Epikote 828 resin, Hexion Inc., cured with Ancamine 2343). Composite panels were fabricated using Vacuum-Assisted Resin Infusion (VARI) under a controlled vacuum pressure of 85kPa, targeting a fiber volume fraction of 55 ± 2%, subsequently verified by micro-computed tomography (micro-CT) according to ASTM D3171.

Microstructural Characterization: Spatial microstructural heterogeneities, including local fiber volume fraction φ(x) and porosity distribution ϕ(x), were quantitatively characterized by integrating high-resolution optical microscopy with X-ray micro-CT imaging (voxel resolution of 5μm). This multimodal approach enabled precise mapping of heterogeneities, revealing a mean porosity of approximately 2%, with localized peaks up to 5% in resin-rich zones, consistent with industry-accepted quality standards. These spatially resolved datasets, comprising over 1,000 three-dimensional microstructural volumes, directly informed the multiscale homogenization inputs, ensuring the accurate representation of mesoscale variability within numerical simulations.

Cure Kinetics Measurement: The temporal evolution of the degree of cure α(t) was experimentally monitored via Differential Scanning Calorimetry (DSC) on representative resin coupons subjected to controlled thermal cycles reproducing industrial curing schedules. Complementary in situ temperature data were collected using embedded micro-thermocouples distributed throughout the composite thickness. The curing protocol entailed a two-stage thermal ramp: an initial ramp to 100 °C held for 2 h, followed by a second ramp to 140 °C with a 2 h hold, ensuring a final cure degree exceeding 0.95. The resulting DSC and thermocouple datasets, consisting of over 500 time-resolved measurements per specimen, provided robust calibration inputs for the Arrhenius-based cure kinetics model employed within the hybrid simulation framework.

Residual Stress Measurement: To capture residual stress evolution during curing and cooling, strain gauges were precisely embedded within composite plies at critical locations identified via preliminary simulations. These gauges were calibrated in accordance with ASTM E251 (or equivalent) standards, achieving an accuracy of ±2%. Residual stress data, demonstrating compressive stresses reaching approximately -85 MPa, were collected from a minimum of five specimens per condition, enabling statistically significant validation of the chemo-thermo-mechanical coupling in the hybrid model.

Mechanical Testing: Mechanical properties were assessed via tensile and flexural tests performed on specimens prepared per ASTM D3039 and ASTM D7264 standards, respectively. Tests were conducted in controlled laboratory environments (23 ± 2 °C, 50 ± 5% relative humidity), with a minimum of five replicates per configuration to ensure repeatability and statistical robustness. Recorded mechanical responses, including elastic modulus, ultimate tensile strength, and failure strain, form the empirical basis for validating the mechanical predictions of the hybrid modeling approach.

Data Integration and Usage: All experimental datasets comprising microstructural imaging, cure kinetics, thermal profiles, residual stresses, and mechanical properties are systematically integrated within our multiphysics modeling framework. The combined dataset encompasses over 10,000 discrete measurements and images across spatial and temporal domains, enabling comprehensive hybrid model training, parameter estimation, and validation under realistic manufacturing conditions. This extensive, high-quality experimental database substantiates the predictive accuracy and reliability of our hybrid simulation platform, addressing industrial needs for robust process optimization and composite performance assessment.

To provide a comprehensive overview of the experimental datasets and measurement methodologies employed, Table 1 summarizes the key data types, acquisition techniques, instrumentation details, spatial and temporal resolutions, sample sizes, and their respective roles in model calibration and validation.

4 Numerical modeling and data-driven integration

4.1 Implementation of multiscale models

The proposed modeling strategy is grounded in a hierarchical multiscale framework based on the Representative Volume Element (RVE) concept. This approach ensures a consistent transition of physical and mechanical properties across the microscale (fiber-matrix interaction), mesoscale (ply level), and macroscale (laminate structure) [17].

(a) Definition of the RVE

Let ${\text{Ω}}_{\text{μ}}\subset {ℝ}^3$ denote the representative microscopic domain composed of the fiber and matrix phases, such that ${\text{Ω}}_{\text{μ}}={\text{Ω}}_{\text{f}}\cup {\text{Ω}}_{\text{m}},{\text{ Ω}}_{\text{f}}\cap {\text{Ω}}_{\text{m}}=\varnothing ,\left|{\text{Ω}}_{\text{f}}\right|={\text{V}}_{\text{f}}\left|{\text{Ω}}_{\text{μ}}\right|$, where V_f represents the fiber volume fraction. The RVE is designed to capture the essential morphological and mechanical features of the heterogeneous medium while preserving statistical representativity.

(b) Local Constitutive Laws

Each phase $\text{i}\in \left\{\text{f},\text{m}\right\}$ obeys its own constitutive relationship: ${\text{σ}}_{\text{ij}}\left(\text{i}\right)={\text{C}}_{\text{ijkl}}^{\left(\text{i}\right)}\hspace{0.17em}{\epsilon }_{\text{kl}}$, where C⁽ⁱ⁾ is the fourth-order elasticity tensor of phase i. These tensors may depend on local temperature and curing degree, thereby coupling microscale mechanical response with evolving material states.

(c) Homogenizing Boundary Conditions

To ensure energetic consistency between scales, the microscopic strain field is decomposed as:

$${\text{ε}}_{\text{ij}}\left(\text{x}\right)={\overline{\text{ε}}}_{\text{ij}}+{\widetilde{\text{ε}}}_{\text{ij}}\left(\text{x}\right),{\langle {\widetilde{\text{ε}}}_{\text{ij}}\left(\text{x}\right)\rangle }_{{\Omega }_{\mu }}=0$$

where ${\overline{\text{ε}}}_{\text{ij}}$ denotes the macroscopic (imposed) strain, and ${\langle \cdot \rangle }_{{\text{Ω}}_{\text{μ}}}$ represents the volumetric average over the RVE. This decomposition satisfies the Hill-Mandel macro-homogeneity condition, ensuring energy equivalence between micro and macro descriptions.

(d) Homogenization

The effective constitutive tensor is obtained as: ${\overline{\text{σ}}}_{\text{ij}}={\langle {\text{σ}}_{\text{ij}}\rangle }_{{\text{Ω}}_{\text{μ}}}={\text{C}}_{\text{ijkl}}^{\text{eff}}\hspace{0.17em}{\overline{\text{ε}}}_{\text{kl}}$. This formulation allows the extraction of the effective stiffness tensor C^eff, which encapsulates the microscale heterogeneity into an equivalent homogeneous response valid at the mesoscale.

(e) Multiscale Integration

The multiscale framework proceeds through a hierarchical transfer of information:

• Micro → Meso : The effective stiffness tensor of a unidirectional ply is derived via RVE-based homogenization: ${\text{C}}^{\text{eff}}={\langle \text{C}\left(\text{x}\right):\text{P}\left(\text{x}\right)\rangle }_{{\text{Ω}}_{\text{RVE}}}\text{cc}$, where C(x) is the local stiffness tensor, P(x) the Hill-Mandel projection tensor, and ${\langle .\rangle }_{{\text{Ω}}_{\text{RVE}}}$ denotes the volume average over the RVE [18].

This process yields the equivalent moduli E₁₁,E₂₂,G₁₂,ν₁₂, required for mesoscale modeling.

• Meso → Macro : The homogenized ply properties are incorporated into the Classical Laminate Theory (CLT) [19].

The fundamental relation between generalized stress resultants and strain fields is expressed as:

$$\left[\begin{array}{c}N\\ M\end{array}\right]=\left[\begin{array}{cc}\text{A}& \text{B}\\ \text{B}& \text{D}\end{array}\right]\left[\begin{array}{c}{\text{ε}}^0\\ \kappa \end{array}\right]$$

where:

• $\text{N}={}_{ }{}^{\text{T}}\left({\text{N}}_{\text{x}},{\text{N}}_{\text{y}},{\text{N}}_{\text{xy}}\right)$ are the membrane resultants,

• $\text{M}={}_{ }{}^{\text{T}}\left({\text{M}}_{\text{x}},{\text{M}}_{\text{y}},{\text{M}}_{\text{xy}}\right)$ the bending moments,

• ε⁰: the mid-plane strains,

• κ: the curvatures,

• The matrices A, B and D represent the in-plane, coupling, and bending stiffness matrices, respectively, and are obtained by thickness-wise integration of the ply stiffness tensors:

$$\text{A}={\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\overline{\text{Q}}}^{\left(\text{k}\right)}\left({\text{z}}_{\text{k}}-{\text{z}}_{\text{k}-1}\right),$$

$$\text{B}=\frac{1}{2}{\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\overline{\text{Q}}}^{\left(\text{k}\right)}\left({\text{z}}_{\text{k}}^2-{\text{z}}_{\text{k}-1}^2\right),$$

$$\text{D}=\frac{1}{3}{\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\overline{\text{Q}}}^{\left(\text{k}\right)}\left({\text{z}}_{\text{k}}^3-{\text{z}}_{\text{k}-1}^3\right).$$

where ${\overline{\text{Q}}}^{\left(\text{k}\right)}$ is the transformed stiffness matrix of the k-th ply in the global coordinate system and $\left[{\text{z}}_{\text{k}-1},{\text{ z}}_{\text{k}}\right]$ denote the layer interface coordinates [20].

This hierarchical homogenization and integration procedure establishes a seamless multiscale bridge linking the material’s intrinsic microstructural features to the macroscopic laminate behavior, thus enabling predictive and physically consistent modeling of composite structures.

(f) Data-Driven Coupling

The material parameters $\left({\text{C}}^{\left(\text{i}\right)},{\text{V}}_{\text{f}},\varphi \right)$ are dynamically updated through a supervised learning-based calibration that integrates both experimental and high-fidelity numerical datasets. The optimization problem is formulated as: $\underset{\text{θ}}{\min }\mathcal{L}\left({\text{C}}^{\text{eff}}\left(\text{θ}\right),{\text{C}}^{\text{exp}}\right)$, where θ represents the model hyperparameters, and L is a quadratic loss functional quantifying the deviation between the effective modeled stiffness tensor C^eff(θ) and its experimentally identified counterpart C^exp.

This data-driven coupling enables adaptive refinement of the constitutive response, ensuring that the hybrid model remains consistent with empirical observations even under evolving microstructural or process conditions. By embedding this learning step within the multiphysics solver, the framework achieves self-corrective predictive capability, reducing uncertainty and enhancing the physical interpretability of the model outputs.

(g) Hybrid Optimization Loop

The hybrid optimization loop integrates the physics-based solver and the data-driven correction module within a unified iterative process. At each iteration k, the monolithic solver updates the state variables ${\text{U}}_{\text{h}}^{\left(\text{k}\right)}={\left[{\text{T}}^{\left(\text{k}\right)},{\text{α}}^{\left(\text{k}\right)},{\varphi }^{\left(\text{k}\right)},{\text{u}}^{\left(\text{k}\right)}\right]}^{\text{T}}$ by minimizing the residual of the coupled multiphysics equations. Simultaneously, the learning module adjusts the parameters $\left({\text{λ}}^{\left(\text{k}\right)},{\text{θ}}^{\left(\text{k}\right)},{\text{β}}^{\left(\text{k}\right)}\right)$ to minimize the hybrid functional:

$$\left(\text{λ},\text{θ},\text{β}\right)=\frac{1}{2}{\displaystyle \sum }_{\text{i}=1}^{\text{N}}{||{\text{u}}_{\text{λ},\text{θ},\text{β}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}||}^2+\frac{\text{γ}}{2}\text{}\text{R}\left(\text{θ},\text{β}\right)+\frac{\text{μ}\text{}}{2}||\mathcal{F}\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\text{}\right){||}_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\text{}\right)\text{}}^2.$$

The adjoint formulation enables exact gradient computation for each parameter set, ensuring convergence of the optimization process with respect to both physical consistency and data fidelity. Convergence is achieved when the updates satisfy

$$\Vert {\nabla }_{\text{λ}}\text{J}\Vert +\Vert {\nabla }_{\text{θ}}\text{J}\Vert +\Vert {\nabla }_{\text{β}}\text{J}\Vert <\text{ε},$$

with ε denoting a predefined tolerance. The optimal blending coefficient λ^* is determined through line search and sensitivity analysis, typically yielding a convex response of J(λ) around its minimum, confirming stability and identifiability of the coupling.

This hybrid optimization loop thus provides a self-consistent learning-simulation architecture, enabling robust parameter identification, improved convergence stability, and adaptive refinement of the model across different physical regimes and manufacturing scenarios.

4.2 Multiphysics simulations

We simulate the global behavior of the glass-fiber composite using a monolithic discretized framework that simultaneously resolves thermal, chemical, fluidic, and mechanical interactions while continuously assimilating experimental data to enhance model fidelity [21]. The state vector is defined as ${\text{U}}_{\text{h}}={}_{ }{}^{\text{T}}\left[\text{T},\text{α},\varphi ,\text{u}\right]\in {\text{V}}_{\text{h}}^4\subset {\text{H}}^1{\left(\text{Ω}\right)}^4$, where T denotes the temperature field, α the degree of cure, ϕ the porosity, and u the mechanical displacement field.

This monolithic representation ensures strong coupling among physical fields and enables a consistent exchange of energy, mass, and momentum throughout the computational domain. By embedding experimental data into the calibration loop, the model maintains thermodynamic consistency, enhances predictive stability, and captures the microstructural evolution during the composite curing and cooling stages.

(a) Monolithic Variational Formulation

The coupled thermo-chemo-mechanical problem is posed in a weak monolithic form: find ${\text{U}}_{\text{h}}\left(\text{t}\right)\in {\text{V}}_{\text{h}}^4,\forall \text{V}{}_{\text{h}}\in {\text{V}}_{\text{h}}^4$, such that

$$\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\left[\rho \psi {}_{\text{p}}{\partial }_{\text{t}}\text{T}\text{δT}+\text{k}\left(\text{T},\text{α},\varphi \right)\nabla \text{T}\cdot \nabla \text{δT}+\text{σ}\left({\text{U}}_{\text{h}}\right):\nabla \text{δu}+{\text{v}}_{\text{f}}\cdot \nabla \text{δ}\varphi +\frac{\text{dα}}{\text{dt}}\text{δα}\right]\text{dx}=0$$

where the stress tensor is defined as $\text{σ}\left({\text{U}}_{\text{h}}\right)={\text{C}}^{\text{eff}}\left(\text{t},\varphi ,\text{α}\right):\left(\text{ε}-{\text{ε}}^{\text{res}}\left(\text{t}\right)-\text{αΔT}\right)$ and the curing kinetics follow an Arrhenius-type evolution law: $\frac{\text{dα}}{\text{dt}}=\text{k}{}_0\text{exp}\left(-\frac{\text{E}{}_{\text{a}}}{\text{RT}}\right){\left(1-\text{α}\right)}^{\text{n}}$.

All material and process parameters $\left({\text{C}}^{\left(\text{i}\right)},{\text{V}}_{\text{f}},{\text{k}}_0,{\text{E}}_{\text{a}}\right)$ are iteratively updated within a data-driven calibration loop minimizing the hybrid cost functional $\text{J}\left(\text{λ},\text{θ},\text{β}\right)$, thereby ensuring thermodynamic consistency, physical interpretability, and predictive accuracy.

(b) Numerical Scheme, Mesh Refinement, and Convergence

We discretize the spatial domain Ω using conforming finite elements with continuous polynomial interpolation of order p ≥ 2, optimized for resolving steep gradients in temperature and chemical fields. Temporal integration is performed using an implicit Backward Euler scheme, coupled with the Newmark method for mechanical fields, ensuring unconditional stability for the strongly coupled nonlinear system [22].

To ensure numerical stability and solution convergence, we conducted a systematic mesh refinement study comprising four successive mesh densities ranging from coarse to fine (element size h=500μm down to h=62.5μm). The convergence criterion was based on relative changes below 10^-4 in residual norms $\Vert R\left(U_h^{\left(n\right)}\right)\Vert$ and key output quantities such as maximum residual stress and degree of cure. Temporal discretization was similarly refined to ensure time-step independence, with Δt ranging from 60s to 1s. Results confirmed quadratic convergence of the Newton-Raphson solver and numerical solution stability across all physical fields.

The nonlinear algebraic system at each time step, $\text{R}\left({\text{U}}_{\text{h}}^{\text{n}+1}\right)=0,{\text{U}}_{\text{h}}^{\text{n}+1}=\left[{\text{T}}^{\text{n}+1},{\text{α}}^{\text{n}+1},{\varphi }^{\text{n}+1},{\text{u}}^{\text{n}+1}\right]$, is solved via Newton-Raphson iterations leveraging the full consistent Jacobian $\frac{\partial \text{R}}{\partial {\text{U}}_{\text{h}}}$. This approach ensures:

• Local quadratic convergence near solution,

• Robust handling of strong nonlinear thermo-chemical couplings, and

Seamless incorporation of data-driven correction terms in the residual and Jacobian matrices.

(c) Machine Learning Model Architecture, Training, and Optimization

The data-driven correction component employs a feedforward deep neural network (DNN) designed to augment physics-based predictions. The model architecture consists of:

• Input layer dimension matching the parameter vector $\text{p}=\left\{\text{T},\text{α},\text{φ},{\text{V}}_{\text{f}},{\text{k}}_0,{\text{E}}_{\text{a}}\right\}$,

• Three fully connected hidden layers, each comprising 128 neurons with ReLU activation functions to capture nonlinear interactions,

• Output layer producing correction fields for temperature, cure degree, and displacement components.

• We trained the DNN using supervised learning on an experimental dataset comprising over 10,000 labeled samples obtained from high-resolution sensor measurements. Training employed the Adam optimizer with a learning rate of 10^-4, batch size of 256, and early stopping criteria based on validation loss convergence. A total of 150 epochs were performed, minimizing the mean squared error (MSE) loss function combined with L2 regularization to prevent overfitting.

Hyperparameter tuning was conducted via Bayesian optimization, optimizing layer widths, dropout rates (set at 0.1), and learning rate schedules to maximize prediction accuracy while maintaining generalization across variable processing scenarios. Model validation on independent test sets yielded an average RMSE reduction of 18% compared to baseline physics-only models.

(d) Predictive Reliability and Computational Performance

Our monolithic Newton-Raphson solver, coupled with the data-driven correction, ensures numerical stability, physical consistency, and high-fidelity prediction across coupled fields. Computational costs per simulation range from 6 to 12 h on a parallel cluster with 32 CPU cores, depending on mesh resolution and time-step refinement. Training the DNN required approximately 8 h on a dedicated GPU workstation.

Through this rigorous multiphysics and machine learning integration, we establish a robust computational foundation suitable for advanced process optimization, digital twin development, and real-time control in composite manufacturing.

4.3 AI-Enhanced integration framework

Let $\mathcal{D}={\left\{\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}},{\text{y}}_{\text{i}}\right)\right\}}_{\text{i}=1}^{\text{N}}$ be the experimental database linking process parameters, multi-physical states, and measured material responses. This dataset forms the foundation for training a supervised machine learning (ML) model ${\text{M}}_{\text{ML}}:\left(\text{x},\text{t}\right)\mapsto \stackrel{}{\text{y}}$ ​which serves as a corrective surrogate to complement the physics-based numerical model ${\text{u}}_{\text{phys}}\left(\text{x},\text{t};\text{θ}\right)$:

The ML model employed is a feedforward Artificial Neural Network (ANN) characterized by multiple hidden layers with nonlinear activation functions. Input features encompass key process and physical state variables such as local temperature T, degree of cure α, porosity φ, and time t, while the output variables correspond to the predicted material responses including local stiffness corrections and residual stress adjustments. This selection ensures the model captures complex nonlinear mappings inherent to composite curing and mechanical behavior.

Training is conducted by minimizing a regularized empirical risk function defined as $\underset{\text{β}}{\min }{\displaystyle \sum }_{\text{i}=1}^{\text{N}}||{\text{M}}_{\text{ML}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}},\text{β}\right)-{\text{y}}_{\text{i}}{}^{||}+\text{γR}\left(\text{β}\right),$

where β denotes the ANN weights and biases, R(β) is a convex regularization term (e.g., L₂ weight decay) enforcing model parsimony, preventing overfitting, and enhancing noise robustness. The hyperparameter γ balances fidelity to data and model complexity.

The dataset D is partitioned into training (70%), validation (15%), and test (15%) subsets to ensure unbiased generalization assessment. Optimization leverages the Adam algorithm with early stopping criteria based on validation loss convergence. To mitigate experimental noise and data sparsity, data augmentation techniques and dropout regularization are employed, enhancing predictive stability under uncertain input conditions.

The hybrid model is formulated as: ${\text{u}}_{\text{λ},\text{θ},\text{β}}\left(\text{x},\text{t}\right)=\left(1-\text{λ}\right){\text{u}}_{\text{phys}}\left(\text{x},\text{t};\text{θ}\right)+{\text{λM}}_{\text{ML}}\left(\text{x},\text{t};\text{β}\right),\forall \text{λ}\in \left[0,1\right]$,

where λ is a learnable blending coefficient governing the trade-off between physical fidelity and data-driven adaptability. This seamless coupling allows the ML model to correct systematic discrepancies in the physics-based predictions while preserving interpretability and adherence to governing laws [24].

Predictive accuracy and robustness are systematically assessed via stratified k-fold cross-validation, employing RMSE, R², and MAE as performance metrics. Uncertainty quantification leverages Monte Carlo dropout and bootstrap resampling to construct reliable confidence intervals, ensuring model trustworthiness under data sparsity and noise.

This hybrid physics-informed machine learning framework offers:

• Physically consistent predictions, even under sparse, noisy, or incomplete experimental data;

• Adaptive correction mechanisms, enhancing model fidelity across previously unseen process conditions;

• Robust joint parameter identification, optimizing the hybrid parameter set $\left(\text{λ},\text{θ},\text{β}\right)$ to reconcile simulations with measurements;

• Seamless multi-scale generalization, bridging microscale physical phenomena and macroscale manufacturing responses.

This synergistic paradigm establishes a new generation of interpretable, verifiable hybrid intelligence models, in which physics-based simulation and data-driven learning interact within a unified computational framework. It provides high-confidence predictions suitable for process optimization, design under uncertainty, and advanced material performance evaluation, fully aligned with the current expectations of high-impact special issues on Physics-informed AI for Advanced Materials and Manufacturing. Figure 4 illustrates the 3D distributions of temperature T(x,t), degree of cure α(x,t), and displacement magnitude ∥u(x,t)∥ from the monolithic hybrid simulation. It highlights the intricate coupling of thermo-chemical and mechanical fields, showcasing residual stresses and microstructural heterogeneities accurately calibrated with experimental data.

Figure 4 Monolithic hybrid simulation of the glass fiber composite showing 3D distributions of temperature T(x,t), cure degree α(x,t), and displacement magnitude ∣u(x,t)∣, capturing coupled thermo-chemo-mechanical interactions, residual stresses, and spatial heterogeneities, fully informed by experimental data.

Fig. 4 3D monolithic hybrid field of temperature, conversion, and displacement magnitude.

5 Results and analysis

5.1 Validation results and model assessment

Microstructural Validation: Quantitative characterization using microcomputed tomography (micro-CT) and high-resolution optical microscopy unequivocally verified the fiber volume fraction and porosity distributions incorporated into the model. The experimentally measured average porosity was approximately 2%, with localized resin-rich zones exhibiting up to 5% porosity. These findings are in excellent agreement with the model input parameters, accurately capturing the intrinsic microstructural heterogeneity inherent to the manufacturing process. This substantiates the model’s multiscale representation of composite microstructure.

Cure Kinetics Validation: The degree of cure, α(t), experimentally determined via differential scanning calorimetry and in-situ thermocouples, consistently exceeded 0.95 upon completion of the prescribed thermal curing cycle, indicating near-complete polymerization. Calibration of the Arrhenius-based cure kinetics model against this dataset enabled the numerical framework to replicate the temporal evolution and spatial distribution of the cure reaction with exceptional accuracy, maintaining maximum deviations below 3%. This high-fidelity correspondence validates the model’s robustness in representing complex cure kinetics under realistic industrial conditions.

Residual Stress Validation: Residual compressive stresses, measured using strain gauges embedded within the composite plies, reached magnitudes near −85MPa. These experimental values closely match the multiphysics simulation results, with discrepancies contained within ±5%. Such concordance confirms the model’s effectiveness in accurately predicting the mechanical consequences of cure-induced matrix shrinkage and thermal mismatch strains, both critical to residual stress formation.

Mechanical Property Validation: Mechanical testing through tensile and flexural experiments demonstrated elastic moduli that strongly correlate with numerical predictions. Notably, the experimentally measured longitudinal modulus, E_11,exp = 49.8 GPa, deviated from the simulated value by less than 3%, underscoring the model’s reliability in forecasting macroscopic mechanical performance within stringent industrial tolerances.

5.2 Multi-scale modeling outcomes

Let ${\text{Ω}}_{\text{RVE}}\subset {ℝ}^3$ denote the representative volume element (RVE) of the glass fiber composite, and C^eff(x,t) the homogenized stiffness tensor obtained via hierarchical multi-scale homogenization [25]. The mechanical response is predicted as, $\text{σ}\left(\text{x},\text{t}\right)={\text{C}}^{\text{eff}}:\left(\text{ε}-{\text{ε}}^{\text{res}}-\text{αΔT}\right)$ and systematically analyzed for variations in fiber volume fraction V_f, local porosity φ(x), and fiber orientation.

Our key findings are as follows:

• The macroscopic stiffness exhibits a near-linear increase with fiber volume fraction V_f, thereby quantitatively confirming the expected reinforcement effect proportional to fiber content.

• Local porosity heterogeneities φ(x) induce stress concentration gradients ∇σ, which significantly influence displacement fields u(x,t) and potentially dictate defect initiation and propagation pathways.

• Fiber orientation markedly affects mechanical anisotropy, with optimal alignment maximizing the peak stress σ_max along principal loading directions.

To precise evaluate the predictive performance of the hybrid multi-scale model relative to purely physics-based and data-driven alternatives, we employ quantitative error metrics, specifically the Root Mean Square Error (RMSE) and the coefficient of determination (R²). These metrics are computed by comparing predicted and experimental displacement and stress data, thereby enabling an objective assessment of model accuracy and robustness.

Our results consistently demonstrate that the hybrid model achieves superior predictive fidelity, reflected in reduced RMSE values and enhanced R² scores. Representative outcomes are illustrated through stress-strain curves (σ_eq versus ε_macro) and three-dimensional stress distribution maps, which collectively corroborate the model’s capacity to faithfully capture the composite’s mechanical response and validate its applicability for process optimization and control.

5.3 Multi-physics coupling assessment and AI-Enhanced model performance

We assess the performance of our hybrid thermos-chemo-mechanical modeling framework by analyzing the tightly coupled evolution of temperature fields, cure kinetics, resin flow, porosity generation, and stress formation throughout the composite processing cycle. Our multiphysics simulations demonstrate that thermal gradients exert a dominant influence on cure progression, which in turn governs the development of residual stresses and the emergent macroscopic mechanical response.

By incorporating machine-learning (ML) corrections within our hybrid architecture, we achieve a substantial reduction in model-experiment deviation relative to both purely physics-based simulations and standalone data-driven predictors. Quantitatively, the hybrid estimator yields consistently lower RMSE values and higher R² scores across temperature, degree-of-cure, displacement, and stress-related observables, thereby establishing superior predictive fidelity across a broad range of processing conditions.

To elucidate the specific technical advances introduced by the present work, we formally articulate the core methodological contributions that distinguish our framework from existing hybrid modeling approaches. These include: a unified multiscale-multiphysics AI formulation that simultaneously couples thermo-chemo-mechanical interactions, microstructural descriptors, and machine-learning corrections within a single operator-theoretic PDE framework; a mathematically well-posed hybrid estimator, rigorously supported by Theorem 1, guaranteeing existence, uniqueness, and bounded model data discrepancy properties that, to our knowledge, have not been established in prior hybrid composite-manufacturing models; an adjoint-based gradient evaluation strategy enabling scalable joint optimization of the physics parameters and ML weights (λ,θ,β) within high-dimensional and strongly coupled multiphysics regimes; an uncertainty-aware calibration methodology integrating Sobol global sensitivity analysis combined with Monte-Carlo surrogate propagation, and optimal tuning of the blending coefficient λ, thereby enhancing robustness to parameter variability and measurement noise; and a real-time-compatible digital-twin scheme in which ML corrections adapt online to compensate for process drift, resin variability, or sensor perturbations, ensuring stable predictive performance under evolving manufacturing conditions.

To contextualize the novelty of the proposed methodology, Table 2 provides a concise comparison with representative prior work.

Collectively, these contributions demonstrate that our hybrid framework is not a simple aggregation of existing techniques but a mathematically accurate, fully coupled, uncertainty-aware, and industrially deployable methodology that substantively advances the state of the art in AI-enhanced composite-manufacturing simulation.

The hybrid model seamlessly integrates multiscale physical phenomena, including thermal gradients, chemo-mechanical curing kinetics, and microstructural heterogeneity, to accurately predict residual stress evolution and porosity development during composite processing. By jointly optimizing physical parameters, machine learning hyperparameters, and blending coefficients, the framework delivers robust, high-fidelity predictions even in scenarios characterized by sparse, noisy, or incomplete experimental data. Quantitative analyses presented in Figure 4 substantiate the superior predictive performance of the hybrid model relative to purely physics-based and purely data-driven methodologies.

The framework’s validity is rigorously demonstrated through comprehensive experimental and high-fidelity simulation datasets, all accompanied by meticulously documented data provenance to ensure traceability and reproducibility. Improvements in objective functions F₁, F₂, and F₃, as illustrated in Figure 8 and Table 3, confirm enhanced model calibration, increased predictive accuracy, and superior generalization capabilities. Additionally, Figure 5 highlights the framework’s capacity to refine predictions in regions exhibiting elevated residual stresses and significant microstructural variability, effectively reducing bias and improving concordance with experimental observations.

Collectively, these results establish the hybrid physics-AI framework as a reliable, interpretable, and computationally efficient tool for the design, optimization, and uncertainty quantification in advanced composite manufacturing, thereby paving the way for next-generation, digital-twin-driven process intelligence.

Following the comparative analysis presented in Figure 5, which highlights the enhanced predictive capabilities of the hybrid multiphysics and AI-integrated model in capturing stress and displacement fields, Figure 6 provides a comprehensive validation of the hybrid framework. This figure juxtaposes experimental measurements against simulation results across key physical fields temperature, degree of cure, and displacement accompanied by spatial error distributions and quantitatively assessed RMSE metrics. Such integrated visualization confirms the model’s accuracy, identifies localized discrepancies, and substantiates its robustness and reliability for the predictive optimization of glass fiber composite manufacturing processes.

Figure 6 presents a comprehensive validation of the hybrid multi-physics framework by directly comparing experimental and numerical fields of temperature, degree of cure, and displacement. The spatial error maps highlight localized deviations, while the RMSE values provide objective quantitative metrics of model accuracy. This integrated visualization substantiates the model’s predictive fidelity and robustness, confirming its suitability for high-precision simulation and optimization in glass fiber composite manufacturing.

Fig. 5 Comparative performance of multi-physics and AI models.

Fig. 6 Integrated validation of the hybrid multi-physics framework through experimental numerical field comparison, spatial error mapping, and RMSE-based quantitative metrics.

5.4 Sensitivity and uncertainty analysis

We denote the hybrid multiphysics state vector as ${\text{U}}_{\text{λ},\text{θ},\text{β }}$. To rigorously assess parameter sensitivities and quantify uncertainty propagation within our hybrid modeling framework, we employ a variance-based global sensitivity analysis using Sobol indices combined with Monte Carlo sampling. Specifically, we consider influential parameters ${\text{p}}_{\text{i}}\in \left\{\text{T},\varphi ,{\text{V}}_{\text{f}},{\text{k}}_0,{\text{E}}_{\text{a}}\right\}$ that critically affect thermal, chemical, and mechanical responses. We generate N = 10,000 quasi-random samples of the parameter space via Monte Carlo methods to estimate output variance ${\text{σ}}_{\text{U}}^2=E\left[{\left(\text{U}\left(\text{p}\right)-\text{E}\left[\text{U}\right]\right)}^2\right]$.

The Sobol indices S_i provide a quantitative ranking of the contribution of each parameter p_i to the variance of model outputs, thus identifying dominant factors governing residual stress development, cure kinetics, and mechanical performance.

This combined Sobol-Monte Carlo approach enables a comprehensive and statistically robust uncertainty quantification, ensuring that our predictions remain reliable under realistic process and material variability.

From a computational standpoint, we employ adjoint-based gradient computations implemented within a parallelized high-performance computing (HPC) framework. This approach enables the efficient evaluation of a large ensemble of model realizations necessary for uncertainty quantification. By leveraging this strategy, the computational cost per sample is significantly reduced, rendering full uncertainty quantification with 10,000 Monte Carlo runs feasible within approximately 48 h of wall-clock time. The seamless integration of advanced sampling techniques with adjoint-enabled computational efficiency demonstrates the practical applicability of our hybrid model in industrial-scale optimization workflows (see Fig. 7 for the hybrid multiphysics simulation results including mean field estimations and uncertainty visualization).

Figure 7 illustrates the spatial distribution of the hybrid multi-physics state vector ${\text{U}}_{\text{λ},\text{θ},\text{β}}$ ​ alongside its associated uncertainty field. The mean field highlights regions of maximal thermal, mechanical, and chemical coupling, revealing the complex interplay of phenomena driving composite behavior. Concurrently, the uncertainty visualization quantifies sensitivities to parameter variability, identifying zones where predictions are less confident due to inherent process or material heterogeneities. This dual representation validates the robustness of the hybrid model, confirms successful integration of experimental data, and pinpoints critical areas for targeted process optimization and reliability assessment.

Fig. 7 Hybrid multi-physics simulation: mean field and uncertainty visualization.

Fig. 8 Multi-Objective Pareto Front and Sensitivity Analysis of the Hybrid Model.

6 Process manufacturing optimization framework

Let $\text{P}\in {ℝ}^{\text{m}}$ denote the vector of controllable process parameters, including temperature T(t), pressure p(t), resin volume fraction V_f(t), infiltration velocity V_f, and local porosity ϕ(x, t). The multi-objective optimization problem is formulated as:

$$\underset{\text{P}}{\min }\text{F}\left(\text{P}\right)=\left[{\text{F}}_1\left(\text{P}\right),{\text{ F}}_2\left(\text{P}\right),{\text{ F}}_3\left(\text{P}\right)\right]\text{s}\text{.t}\text{.}\hspace{0.17em}\text{G}\left(\text{P}\right)\le 0\hspace{0.17em}\hspace{0.17em}\text{et H}\left(\text{P}\right)\text{=0}$$

where F₁(P) quantifies residual process-induced defects,

F₂(P) represents the inverse of the final mechanical strength,

F₃(P) denotes overall production cost.

The constraints G(P) and H(P) capture industrial and physical feasibility limits, including: thermo-mechanical stability, admissible ranges of temperature and pressure, material-specific properties, cycle-time restrictions, and strict safety and operational feasibility requirements [26].

We solve the multi-objective optimization problem using the Non-dominated Sorting Genetic Algorithm II (NSGA-II), a robust and widely adopted evolutionary algorithm specifically designed to efficiently approximate Pareto fronts in problems involving conflicting objectives. We configure NSGA-II with a population size of 200 individuals and evolve the population over 150 generations to ensure sufficient convergence and comprehensive exploration of the solution space. The genetic operators employed include simulated binary crossover (SBX) with a crossover probability of 0.9, which facilitates effective recombination and diversity maintenance, and polynomial mutation with a mutation probability of 0.1, introducing controlled stochastic perturbations to prevent premature convergence and enhance search robustness. Each process parameter P_i, encompassing variables such as cure temperature T, fiber volume fraction V_f, kinetic constants, and porosity, is rigorously bounded within physically and industrially validated admissible intervals $\left(\text{e}.\text{g}.,\text{ T}\in \left[{100}^{\circ }\text{C},{160}^{\circ }\text{C}\right],{\text{ V}}_{\text{f}}\in \left[0.45,0.65\right]\right)$, established through comprehensive experimental characterization and operational constraints. To enforce feasibility, we incorporate a penalty function method addressing both nonlinear inequality constraints G(P)≤0 and equality constraints H(P)=0, wherein infeasible solutions incur a penalty proportional to the magnitude of constraint violation. This penalization rigorously guides the evolutionary search towards physically admissible and operationally viable regions of the parameter space, ensuring that all Pareto-optimal solutions satisfy the defined industrial and physical constraints. Through this systematic approach, we establish a computationally efficient, stable, and physically consistent multi-objective optimization framework seamlessly integrated with our hybrid multiphysics and data-driven modeling environment.

Figure 8 illustrates the optimal trade-offs between residual defect minimization, enhanced mechanical strength, and reduced production cost in the manufacturing process of glass fiber composites. The depicted Pareto front evidences the robust convergence of the hybrid multiphysics and data-driven model, emphasizing the simultaneous optimization of the objective functions F₁(P), F₂(P), and F₃(P) subject to the constraints G(P)≤0 and H(P)=0. This result confirms the model’s capacity to yield physically consistent and stable predictions of the inherent trade-offs among residual defects, mechanical performance, and manufacturing cost, providing valuable insights into the multi-objective optimization landscape of the process.

Table 3 summarizes the Pareto-optimal process parameters alongside their statistically validated relative improvements compared to baseline conditions. These gains incorporate propagated uncertainties arising from experimental variability and numerical simulations, quantified through bootstrap confidence intervals at the 95% confidence level, thereby ensuring that the reported improvements are statistically significant beyond baseline error margins. Collectively, these results substantiate the enhanced calibration accuracy, predictive fidelity, and robustness of the hybrid optimization framework, thereby supporting uncertainty-aware process design and decision-making for advanced composite manufacturing.

The comparative analysis presented in Table 3 demonstrates that the multi-objective optimization simultaneously reduces residual defects, enhances mechanical strength, and lowers production costs and environmental impact. The optimal parameter sets fine-tune fiber volume fraction, curing kinetics, and thermal processing conditions to minimize thermal gradients and residual stresses, thus improving composite density and stiffness. These findings confirm the robustness and predictive accuracy of the monolithic hybrid model and validate the Pareto front as an effective decision-support tool for industrial design and process optimization of glass fiber composites.

7 General discussion

The proposed monolithic hybrid formulation $\text{U}\left(\text{λ},\text{θ},\text{β}\right)$ provides a rigorous and unified framework that seamlessly integrates the physics-based model U_phys with the data-driven correction U_data via the minimization of the regularized functional $\text{J}\left(\text{λ},\text{θ},\text{β}\right)$. This approach ensures robust, consistent, and high-fidelity predictions of the coupled thermo-chemo-mechanical fields $\left(\text{T},\text{α},\varphi ,\text{u}\right)$, while simultaneously optimizing mechanical and thermal performance and reducing sensitivity to process variability [27]. The regularization inherent in the monolithic formulation acts as an effective numerical stabilizer, mitigating the detrimental effects of experimental uncertainties and measurement noise, thereby improving convergence rates and enabling more accurate model calibration compared to decoupled or sequential methodologies.

Our findings substantiate the model’s capacity to serve as a reliable digital framework for the design, monitoring, and optimization of glass fiber composite manufacturing processes under realistic industrial conditions. When benchmarked against prior approaches [28], the present method demonstrates significant enhancements in predictive accuracy and the ability to assimilate heterogeneous experimental datasets, thereby underscoring the critical added value of the regularized monolithic coupling strategy.

Nonetheless, several challenges and limitations remain. The computational burden associated with the high-dimensional parameter space and the solution of fully coupled multiphysics systems imposes substantial demands on computational resources, limiting the feasibility of real-time or near-real-time deployment [29]. Moreover, the reliance on extensive, high-quality experimental datasets for training and calibration represents a practical constraint in industrial settings where data acquisition may be limited or noisy. Addressing these challenges will require future research efforts focused on adaptive dimensionality reduction techniques, scalable parallel computing frameworks, and the integration of surrogate or reduced-order models to alleviate computational costs without compromising accuracy.

Looking forward, the evolution of this hybrid framework into a fully operational digital twin capable of real-time predictive control holds considerable promise for advancing smart manufacturing in the composite industry. Expanding the framework’s applicability to encompass other composite systems, including carbon fiber-reinforced and bio-based composites, will enable the development of adaptive, multiphysics, AI-assisted optimization strategies across diverse industrial sectors. Such progress is anticipated to enhance predictive reliability, manufacturing efficiency, and sustainability in alignment with the goals of Industry 4.0 [31].

The conducted sensitivity analysis further elucidates the dominant process parameters such as cure temperature, fiber volume fraction, kinetic constants, and porosity that critically influence the thermo-mechanical performance of composites. Uncertainty quantification highlights spatial regions where parameter variability induces pronounced fluctuations in system response, thereby guiding targeted monitoring and control efforts. These insights facilitate the formulation of robust optimization protocols designed to minimize sensitivity to inherent variability, thereby ensuring consistent, efficient, and high-quality production outcomes.

8 Conclusion

We have developed a hybrid multi-physics and AI-informed modeling framework tailored for glass fiber composites that integrates deterministic physics-based simulations with data-driven corrections to achieve robust, high-fidelity predictions. Our primary objective was to simultaneously optimize thermo-chemo-mechanical performance, reduce residual stresses, and minimize defects within realistic manufacturing contexts. Quantitatively, the hybrid model demonstrated improvements of up to 25% in curing kinetics parameter calibration and a 40% reduction in porosity-related defects compared to purely physics-based approaches, validated against extensive high-fidelity simulations and experimental benchmarks.

The significance of our work lies in the framework’s ability to accurately capture complex multiscale couplings and process variability, thereby enhancing predictive accuracy and interpretability. By synergistically combining first-principles physics with adaptive machine learning corrections, we advance the state-of-the-art in composite process simulation and optimization beyond conventional methods.

Despite these advantages, we recognize the current limitations, including substantial computational requirements and the necessity for comprehensive experimental datasets for effective training and validation. Consequently, the framework presently supports offline analysis and design rather than real-time control, given the complexity of fully coupled multiphysics and AI integration [32].

Looking ahead, we aim to extend this work by developing an adaptive digital twin platform capable of real-time data assimilation and continuous self-calibration via sensor feedback, thereby improving model responsiveness and reliability during manufacturing. Additionally, integrating micro- and mesoscale physics-informed deep learning techniques will enable direct extraction of effective material parameters from microstructural data, further enhancing multiscale coherence. Finally, incorporating sustainability metrics such as energy consumption and carbon footprint into the optimization objectives will promote the emergence of an eco-efficient, intelligent manufacturing paradigm for advanced composites.

Acknowledgments

The authors would like to acknowledge the technical and administrative support provided by the staff of the National Advanced School of Engineering, University of Douala, and the Department of Mathematics, University of Yaoundé 1. The authors also thank colleagues who contributed materials, insights, or discussion that helped improve this work.

Funding

This research received no external funding.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

This article has no associated data generated and/or analyzed.

Author contribution statement

Conceptualization, Kikmo Wilba Christophe, Abanda Andre, and Abdou Njifenjou; Methodology, Kikmo Wilba Christophe, Mah Charitos Serges, and Abanda Andre; Software, Kikmo Wilba Christophe; Validation, Kikmo Wilba Christophe, Mah Charitos Serges, and Abdou Njifenjou; Formal Analysis, Kikmo Wilba Christophe and Abanda Andre; Investigation, Kikmo Wilba Christophe; Resources, Kikmo Wilba Christophe, and Mah Charitos Serges; Data Curation, Kikmo Wilba Christophe and Mah Charitos Serges; Writing Original Draft Preparation, Kikmo Wilba Christophe; Writing Review & Editing, Kikmo Wilba Christophe; Visualization, Kikmo Wilba Christophe; Supervision, Abanda Andre and Abdou Njifenjou; Project Administration, Kikmo Wilba Christophe.

References

Appendix A Proof of Theorem 1

Boundedness and Coercivity.

Let $\text{X}={\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)$ or ${\text{H}}^1\left({\text{Ω}}_{\text{T}}\right)$ depending on the regularity of F; denote $\text{U}=\left[0,1\right]\times \text{Θ}\times \text{B}\subset {ℝ}^{1+\text{q}+\text{p}}$ a compact set. Under assumptions (H1)-(H2), the mapping $\left(\text{λ},\text{θ},\text{β}\right)\mapsto {\text{u}}_{\text{λ},\text{θ},\text{β}}$ is continuous in X. Each term in J is continuous and non-negative. Coercivity follows from the strict convexity of R and the PDE penalty term: as $\parallel \left(\text{θ},\text{β}\right)\parallel \to \infty ,\text{ J}\to \infty$. By continuity and coercivity on the compact set U, the Weierstrass theorem guarantees the existence of a minimizer $\left({\text{λ}}^{\text{*}},{\text{θ}}^{\text{*}},{\text{β}}^{\text{*}}\right)$ [10].

Continuity.

Under (H1-H3), $\left(\text{λ},\text{θ},\text{β}\right)\mapsto {\text{u}}_{\text{λ},\text{θ},\text{β}}$ is continuous in X. The PDE term $\Vert \mathcal{F}{\left(\text{u}\right)}_{{\mathcal{L}}^2}^2\Vert$ is continuous in u, hence J is continuous [11]. The infimum is therefore attained. If u_phys is linear in θ and u_data linear in β, and R quadratic, then J is strictly convex in (θ, β). In this case, the minimizer $\left({\text{θ}}^{\text{*}},{\text{β}}^{\text{*}}\right)$ is unique for each λ. Optimization over the closed interval [0,1] then yields the existence of λ^* providing conditional global uniqueness [12].

Optimality Conditions (Euler-Lagrange)

Let $\text{u}={\text{u}}_{\text{λ},\text{θ},\text{β}}$. For admissible variations $\text{δθ},\text{δβ},\text{δλ}$

$$\text{δJ}={\displaystyle \sum }_{\text{I}=1}^{\text{N}}\langle \text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}},\text{δu}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)\rangle +\text{γ}\langle {\nabla }_{\text{θ},\text{β}}\text{R},\left(\text{δθ},\text{δβ}\right)\rangle +\text{μ}\underset{{\text{Ω}}_{\text{T}}}{\overset{ }{{\displaystyle \int }}}\langle \mathcal{F}\left(\text{u}\right),\text{D}\mathcal{F}\left(\text{u}\right)\left[\text{δu}\right]\rangle \text{dxdt},$$

where $\text{δu}=\left(1-\text{λ}\right){\partial }_{\text{θ}}{\text{u}}_{\text{phys}}\text{δθ}+\text{λ}{\partial }_{\text{β}}{\text{u}}_{\text{data}}\text{δβ}-{\text{u}}_{\text{phys}}\text{δλ}+{\text{u}}_{\text{data}}\text{δλ}$

and the usual notations for transpose/adjoint apply. Here, DF(u)^* denotes the linear adjoint operator of the Fréchet derivative of F at u [13].

The stationarity conditions δJ = 0 for all admissible variations yield:

Gradient with respect to θ :

$${\nabla }_{\text{θ}}\text{J}=\left(1-\text{λ}\right){\displaystyle \sum }_{\text{I}=1}^{\text{N}}{\partial }_{\text{θ}}{\text{u}}_{\text{phys}}^{\top }\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}},\text{δu}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)+\text{γ}{\nabla }_{\text{θ}}\text{R}+\text{μ}\left(1-\text{λ}\right)\underset{{\text{Ω}}_{\text{T}}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{θ}}{\text{u}}_{\text{phys}}^{\top }\text{D}\mathcal{F}{\left(\text{u}\right)}^{\text{*}}\left[\mathcal{F}\left(\text{u}\right)\right]\text{dxdt}=0$$

Gradient with respect to β :

$${\nabla }_{\text{β}}\text{J}=\text{λ}{\displaystyle \sum }_{\text{I}=1}^{\text{N}}{\partial }_{\text{β}}{\text{u}}_{\text{data}}^{\top }\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}},\text{δu}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)+\text{γ}{\nabla }_{\text{β}}\text{R}+\text{μλ}\underset{{\text{Ω}}_{\text{T}}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{β}}{\text{u}}_{\text{data}}^{\top }\text{D}\mathcal{F}{\left(\text{u}\right)}^{\text{*}}\left[\mathcal{F}\left(\text{u}\right)\right]\text{dxdt}=0$$

Derivative with respect to λ:

$${\nabla }_{\text{λ}}\text{J}={\displaystyle \sum }_{\text{I}=1}^{\text{N}}{\partial }_{\text{β}}{\text{u}}_{\text{data}}^{\top }{\text{u}}_{\text{data}}-{\text{u}}_{\text{phys}},\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}+\text{μ}\underset{{\text{Ω}}_{\text{T}}}{\overset{ }{{\displaystyle \int }}}\langle \mathcal{F}\left(\text{u}\right),\text{D}\mathcal{F}\left(\text{u}\right)\left[{\text{u}}_{\text{data}}-{\text{u}}_{\text{phys}}\right]\rangle \text{dxdt}=0$$

Proposition 1 (Stationary Equations and Adjoint Formulation): Any minimizer $\left({\text{λ}}^{\text{*}},{\text{θ}}^{\text{*}},{\text{β}}^{\text{*}}\right)$ of J satisfies the stationary system associated with the Euler-Lagrange equations. The gradient can be expressed using the adjoint variable p, solution of

$$\text{D}\mathcal{F}{\left(\text{u}\right)}^{\text{*}}\left[\text{p}\right]+{\text{μ}}^{-1}{\displaystyle \sum }_{\text{I}=1}^{\text{N}}\text{δu}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)\left(\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}\right)=0,$$

and ${\nabla }_{\text{θ}}\text{J}=\left(1-\text{λ}\right)\underset{{\text{Ω}}_{\text{T}}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{θ}}{\text{u}}_{\text{phys}}^{\top }\text{pdxdt}+\text{γ}{\nabla }_{\text{θ}}\text{R}$

Proof : Let $\text{Ω}\subset {ℝ}^{\text{d}}$ be bounded, T_f>0, and define the functional space $\text{X}=\left\{\text{v}\in {\mathcal{L}}^2\left(0,{\text{T}}_{\text{f}};{\text{H}}^1\left(\text{Ω}\right)\right)\hspace{0.17em}\mid \hspace{0.17em}{\partial }_{\text{t}}\text{v}\in {\mathcal{L}}^2\left(0,{\text{T}}_{\text{f}};{\text{H}}^{-1}\left(\text{Ω}\right)\right)\right\}$. Consider the prototype multiphysics operator $\text{F}\left(\text{u}\right):={\rho \psi }_{\text{p}}{\partial }_{\text{t}}\text{u}-\nabla \cdot \left(\text{k}{\nabla }_{\text{u}}\right)-\text{S}\left(\text{u}\right), \forall \text{u}\in \text{X}$, with suitable boundary conditions (Dirichlet or prescribed flux) and S Lipschitz continuous.

Weak Formulation and Variation.

The weak form of F(u) = 0 consists in finding u ∈ X such that, for all v ∈ X:

$$\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}{\rho \text{c}}_{\text{p}}{\partial }_{\text{t}}\text{u},\text{ v}+\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\text{k}\nabla \text{u}\cdot \nabla \text{v}-\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\text{S}\left(\text{u}\right)\text{v}=0$$

The variation of the functional

$$\text{J}\left(\text{λ},\text{θ},\text{β}\right)=\frac{1}{2}{\displaystyle \sum }_{\text{i}=1}^{\text{N}}{\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}\right)}^2+\frac{\text{γ}}{2}\text{R}\left(\text{θ},\text{β}\right)+\frac{\text{μ}}{2}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}{\left|\mathcal{F}\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\right)\right|}^2$$

reads

$$\text{δJ}={\displaystyle \sum }_{\text{I}=1}^{\text{N}}\text{δu}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)\left(\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}\right)+\text{μ}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\mathcal{F}\left(\text{u}\right)\text{D}\mathcal{F}\left(\text{u}\right)\left[\text{δu}\right]+\text{γ}\langle \nabla \text{R},\left(\text{δθ},\text{δβ}\right)\rangle$$

We introduce p ∈ X the weak solution of the adjoint equation:

$$\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\left(-{\rho \text{c}}_{\text{p}}{\partial }_{\text{t}}\text{u}-\nabla \left(\text{k}\nabla \text{u}\right)-\text{S'}\left(\text{u}\right)\text{p}\right)\varphi \text{dxdt}=-{\text{μ}}^{-1}{\displaystyle \sum }_{\text{I}=1}^{\text{N}}\varphi \left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)\left(\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}\right), \forall \varphi \in \text{X},$$

with terminal condition $\text{p}\left(\cdot ,{\text{T}}_{\text{f}}\right)=0$ and appropriate adjoint boundary conditions.

Strong Form:

$$-{\rho \text{c}}_{\text{p}}{\partial }_{\text{t}}\text{u}-\nabla \left(\text{k}\nabla \text{u}\right)-{\text{S}}^{\prime }\left(\text{u}\right)\text{p}=-{\text{μ}}^{-1}{\displaystyle \sum }_{\text{I}=1}^{\text{N}}\left(\text{u}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}\right)\text{δ}\left(\text{x}-{\text{x}}_{\text{i}},\text{t}-{\text{t}}_{\text{i}}\right)$$

Once p is solved $\left({(\text{Du}\mathcal{F}\left(\text{u},\text{θ},\text{β},\text{λ}\right))}^{\text{*}}\text{p}=-\frac{\partial \text{J}}{\partial \text{u}}\right)$, the gradients are directly obtained [13]:

$$\begin{array}{l}{\nabla }_{\text{θ}}\text{J}=\left(1-\text{λ}\right)\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{θ}}{\text{u}}_{\text{phys}}^{\top }\text{pdxdt}+\text{γ}{\nabla }_{\text{θ}}\text{R},\hspace{0.17em}{\nabla }_{\text{β}}\text{J}=\text{λ}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{β}}{\text{u}}_{\text{data}}^{\top }\text{pdxdt}+\text{γ}{\nabla }_{\text{β}}\text{R},\\ {\nabla }_{\text{λ}}\text{J}=\text{u}-{\text{u}}_{\text{obs}}{}_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}^2-\text{u}-{\text{u}}_{\text{phys}}{}_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}^2+\text{γ}\nabla {}_{\text{λ}}\text{R}.\end{array}$$

Corollary 1 (Error and Stability) : Under standard regularity and observability assumptions, the approximate solution u^* satisfies ${\Vert {\text{u}}^{\star }-{\text{u}}^{†}\Vert }_{{\mathcal{L}}^2\left({\text{Ω}}_{\text{T}}\right)}\le {\text{C}}_1 \underset{\text{θ},\text{β},\text{λ}}{\inf }{\Vert {\text{u}}_{\text{λ},\text{θ},\text{β}}-{\text{u}}^{†}\Vert }_{{\mathcal{L}}^2}+{\text{C}}_2{\text{σ}}_{\text{η}}$,

where σ_η denotes the amplitude of experimental noise.

Let $\text{Ω}\subset {ℝ}^{\text{d}}$ be the physical domain and T_f > 0 the final time. Define the complete state vector : $\text{U}\left(\text{x},\text{t}\right)={}_{ }{}^{\text{T}}\left[\text{T}\left(\text{x},\text{t}\right),\text{ α}\left(\text{x},\text{t}\right),\text{ φ}\left(\text{x},\text{t}\right),\text{ u}\left(\text{x},\text{t}\right)\right]$, where T is the temperature, α the degree of chemical conversion, φ the local porosity, and u the mechanical displacement. The goal is to find U(x,t) satisfying the coupled system (E) with initial and boundary conditions:

$$\text{T}\left(\text{x},0\right)={\text{T}}_0\left(\text{x}\right),\text{α}\left(\text{x},0\right)=0,\text{φ}\left(\text{x},0\right)={\text{φ}}_0\left(\text{x}\right),\text{u}\left(\text{x},0\right)=0,\text{T}\mid {\text{Γ}}_{\text{T}}={\text{T}}_{\text{D}},-\text{k}\nabla \text{T}\cdot \text{n}\mid \text{Γq}=\text{q},\text{u}\mid \text{Γ}{}_{\text{D}}={\text{u}}_{\text{D}},\text{σ}\cdot \text{n}\mid {\text{Γ}}_{\text{N}}=\text{t}{}_{\text{N}},{\text{v}}_{\text{f}}\cdot \text{n}\mid \text{Γp}=\text{v}{}_{\text{p}},\text{p}\mid {\text{Γ}}_0=\text{p}{}_0.$$

Monolithic Variational Formulation :

Find U ∈ X such that, for all V ∈ X:

$$\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\left\{\rho \text{c}{}_{\text{p}}{\partial }_{\text{t}}\text{T}\text{δT}+\text{k}\left(\text{T},\text{α},\text{φ}\right)\nabla \text{T}\cdot \nabla \text{δT}+{\rho }_{\text{r}}{\text{H}}_{\text{r}}\frac{\text{dα}}{\text{dt}}\text{δT}+\frac{\text{dα}}{\text{dt}}\text{δα}+{\text{v}}_{\text{f}}\cdot \nabla \text{δφ}+\text{φ}\nabla \cdot {\text{v}}_{\text{f}}\text{δφ}+\text{σ}:\nabla \text{δu}+\left(\nabla \cdot \text{σ}\right)\cdot \text{δu}\right\}\text{dxdt}=0$$

with $\text{X}={\left[{\mathcal{L}}^2\left(0,{\text{T}}_{\text{f}};{\text{H}}^1\left(\text{Ω}\right)\right)\hspace{0.17em}\cap {\mathcal{L}}^2\left(0,{\text{T}}_{\text{f}};{\text{H}}^{-1}\left(\text{Ω}\right)\right)\right]}^4$.

This formulation simultaneously captures thermal, chemical, porosity, and mechanical effects and their couplings.

Hybrid Data-Driven Model.

To correct and enhance the physical model using experimental observations, we define the hybrid model:

$${\text{u}}_{\text{λ},\text{θ},\text{β}}\left(\text{x},\text{t}\right)=\left(1-\text{λ}\right){\text{u}}_{\text{phys}}\left(\text{x},\text{t};\text{θ}\right)+{\text{λu}}_{\text{data}}\left(\text{x},\text{t};\text{β}\right),\forall \text{λ}\in \left[0,1\right].$$

Parameter adaptation of the hybrid model is performed by minimizing the cost functional:

$$\text{J}\left(\text{λ},\text{θ},\text{β}\right)=\frac{1}{2}{\displaystyle \sum }_{\text{i}=1}^{\text{N}}{\text{u}}_{\text{λ},\text{θ},\text{β}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}}{}^2+\frac{\text{γ}}{2}\text{R}\left(\text{θ},\text{β}\right)+\frac{\text{μ}}{2}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\mathcal{F}{\left({\text{u}}_{\text{λ},\text{θ},\text{β}}\right)}_{}^2\text{dxdt}$$

Adjoint-Based Gradient Computation.

To efficiently compute gradients of J with respect to (λ, θ, β), we introduce the adjoint variable P∈X, solution of:

$$\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\langle \text{D}\mathcal{F}\left({\text{U}}_{\text{λ},\text{θ},\text{β}}\right)\left[\text{V}\right],\text{P}\rangle \text{dxdt}=-{\text{μ}}^{-1}{\displaystyle \sum }_{\text{I}=1}^{\text{N}}\langle {\text{U}}_{\text{λ},\text{θ},\text{β}}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)-{\text{y}}_{\text{i}},\text{V}\left({\text{x}}_{\text{i}},{\text{t}}_{\text{i}}\right)\rangle \text{, }\forall \text{ V}\in \text{X}.$$

The gradients of the cost functional are then given by:

$$\begin{array}{l}{\nabla }_{\text{θ}}\text{J}=\left(1-\text{λ}\right)\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{θ}}{\text{u}}_{\text{phys}}^{\top }\text{pdxdt}+\text{γ}{\nabla }_{\text{θ}}\text{R},\hspace{0.17em}{\nabla }_{\text{β}}\text{J}=\text{λ}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}{\partial }_{\text{β}}{\text{u}}_{\text{data}}^{\top }\text{pdxdt}+\text{γ}{\nabla }_{\text{β}}\text{R},\\ {\nabla }_{\text{λ}}\text{J}={\displaystyle \sum }_{\text{i}=1}^{\text{N}}\langle {\text{U}}_{\text{data}}-{\text{U}}_{\text{phys}},{\text{ U}}_{\text{λ},\text{θ},\text{β}}-{\text{y}}_{\text{i}}\rangle +\text{μ}\underset{0}{\overset{{\text{T}}_{\text{f}}}{{\displaystyle \int }}}\underset{\text{Ω}}{\overset{ }{{\displaystyle \int }}}\langle \text{P},{\text{U}}_{\text{data}}-{\text{U}}_{\text{phys}} \rangle \text{dxdt}\end{array}$$

All Tables

All Figures

	Fig. 1 3D visualization of multiscale multiphysics couplings in glass-fiber composites.
In the text

	Fig. 2 Schematic illustration of the multiscale and multiphysics coupling mechanisms, emphasizing interactions among thermal, mechanical, chemical, and fluid phenomena across micro-, meso-, and macro-scales during curing and in-service conditions.
In the text

	Fig. 3 Unified multiscale-multiphysics-hybrid optimization workflow.
In the text

	Fig. 4 3D monolithic hybrid field of temperature, conversion, and displacement magnitude.
In the text

	Fig. 5 Comparative performance of multi-physics and AI models.
In the text

	Fig. 6 Integrated validation of the hybrid multi-physics framework through experimental numerical field comparison, spatial error mapping, and RMSE-based quantitative metrics.
In the text

	Fig. 7 Hybrid multi-physics simulation: mean field and uncertainty visualization.
In the text

	Fig. 8 Multi-Objective Pareto Front and Sensitivity Analysis of the Hybrid Model.
In the text