The range of applications for fiber-reinforced polymer matrix composites continues to expand. These heterogeneous, direction-dependent materials pose significant challenges for modeling, as their mechanical behavior depends on parameters such as layer thickness, fiber orientation, and stacking sequence. Numerous laminate configurations – including various fillers – remain the subject of ongoing research [1,2]. As a result, stress and strain analyses are carried out at both micro- and macro-scales.
The increasing industrial demand for models capable of reliably estimating the strength and mechanical properties of polymer matrix composites – while ensuring repeatability of manufactured components [3] – means that most existing approaches rely on phenomenological hypotheses of damage accumulation. This trend is driven by advances in deformation and failure theories grounded in structural considerations [4], particularly the onset of delamination and the development of microcracks in the reinforcement and matrix [5–8].
Early work by Halpin et al. [9], building on the modified Reifsnider models [10,11], addressed the role of critical and subcritical stress states within layered composites and attempted to correlate the ultimate strength of the composite with specific microdamage mechanisms. These include the models proposed by Song and Otani, and Gamstedt and Sjögren [12].
Despite the variety of available methods for determining or defining material strength, many of them do not adequately characterize polymer matrix composites – particularly honeycomb sandwich structures [13,14]. Optimal mechanical properties for a minimum-mass sandwich panel with an appropriate honeycomb-core-to-panel-mass ratio were presented in [15]. Further work has documented the typical damage mechanisms and locations in sandwich composites (Fig. 1), including matrix cracking, fiber fracture, and core delamination [13, 16].
A well-known example illustrating the severity of such damage mechanisms is the May 2005 crash of the Grumman G-73T Turbo Mallard seaplane operated by Chalk’s Ocean Airways (N2969, built in 1947; Fig. 2) [17]. One of the identified failure modes was the separation of the wing caused by a fatigue crack in the lower spar, previously concealed beneath a layer of sealant.
Grumman Seaplane (G-73/G-73T Wing Box Component) [17].
Other researchers, such as Zhang et al. [18], identify hail impact as a key initiator of delamination in hybrid sandwich structures. According to their findings, delamination is governed primarily by the interlaminar cracking force rather than by fiber distribution. The authors analyze cladding–core debonding under impacts at various velocities (including the influence of corrosion) using a nonlinear FEM model. Damage mechanisms in composite structures, such as surface scratches, dents, and moisture-driven degradation in modern composite aircraft wings, are discussed in [19] and [20], respectively.
In the present study, we investigate the feasibility of processing fatigue-test results to estimate the residual undamaged area of a composite material using a Poisson process model for the initiation of weak microvolumes (WMVs). Detecting subsurface defects is particularly challenging because the true shape and depth of real delamination regions cannot be directly determined. This uncertainty complicates the development of appropriate repair procedures for composite structures. The proposed mathematical model enables the construction of a probabilistic description of fatigue-damage propagation at an arbitrary positive stress ratio and facilitates the prediction of both residual strength and the extent of damage.
We assume that failure of a composite specimen under cyclic loading occurs after the destruction of n weak microvolumes (WMVs). Each WMV consists of two components: perfectly elastic (brittle) longitudinal fibers or fiber bundles, and a plastic component (the matrix), where plastic strains accumulate during cyclic loading (Fig. 3). It should be emphasized that, in addition to the longitudinal fibers, the matrix contains all other constituents of the composite – namely the polymer matrix itself and all plies with orientations other than the longitudinal direction.
Within the model, we assume that cyclic loading reduces the number of elastic load-bearing elements within each microvolume by a quantity rR. Once these elastic elements can no longer sustain the applied load, their elastic component fails – and consequently, the entire microvolume (and ultimately the specimen) fails. In Figure 3, the inclined lines schematically indicate the possibility of accumulating irreversible strains up to the level εy. When this level is reached, the microvolume – and therefore the specimen – is considered destroyed.
Note that this graphical interpretation is only symbolic when applied to composite materials. It originates from descriptions more typical of metals, in which plastic strain accumulation corresponds to a distinct “flow act” (e.g., displacement of slip planes). In classical formulations, each change in the internal state is treated as a transition to a new state of the system [21,22].
However, how justified is this assumption for a modern composite structure? If we treat an alloy of different metals as a composite, the analogy may indeed hold (Fig. 4). In the present study, however, we examine composite systems in which plastic-deformation zones cannot develop. The stiffness of such materials prevents plastic flow, and failure is therefore governed by brittle mechanisms. As in ceramics and similar brittle materials, damage does not manifest as a gradual increase in crack length; instead, it appears as sudden area-type damage involving sliding or separation across the material volume.
Example of delamination and approximate linear size in relation to the initial damage width.
When considering damage in a honeycomb sandwich composite, delamination typically initiates at the edges of the structure and propagates inward. This propagation does not occur within the fibers themselves, but rather along the bond line between the honeycomb core and the skin. For this reason, assessing the final damage extent based solely on visible skin-edge indications is highly unreliable. Delamination may also originate within the interior of the honeycomb core, but this requires specific conditions – such as the presence of Euler-type compressive stresses and localized weak adhesion between the core and skin. Under such conditions, the honeycomb–skin interface becomes a preferential path for delamination growth.
Figure 5 shows the onset and development of delamination in a flap component. At the early stage – before the damage reaches the honeycomb core – the deformation is absorbed partially by the matrix and partially by the composite fibers. While surface damage may appear minor, primary delamination is significantly less critical than the stage at which the process penetrates the honeycomb core. Crucially, once delamination extends into the honeycomb, the ratio of its linear propagation length to the visible width at the exposed skin edge may exceed 20:1. When compared with the skin thickness, this ratio may reach 200:1. This illustrates why subsurface damage in honeycomb composites is particularly dangerous: it can spread over large areas while remaining invisible. Such subsurface damage may take years to develop. Although it can be detected using non-destructive testing (NDT), this requires intensive and costly maintenance procedures. Repairing honeycomb structures is also technically challenging.
Damage to an aircraft flap element, before and after preparation for repair.
Figure 5 illustrates this issue: after the skin is removed during maintenance, the actual damage area – including completely destroyed honeycomb cells and partially degraded core regions – proves to be far larger than what is visible from the exterior. Prior to skin removal, the extent of damage cannot be visually assessed and must be inferred solely through NDT methods.
In honeycomb sandwich structures, one of the primary imperfections arises from the limited contact area between the honeycomb core and the skins. Core–skin bonding occurs mainly at the nodal attachment points, where individual honeycomb cells intersect. These nodes largely determine the residual peel strength and initiate the onset and progression of delamination, as clusters of weak nodes form preferential failure paths.
Imperfections in the honeycomb core – such as dislocations, vacancies, and grain-boundary irregularities – also influence damage initiation. Their occurrence depends on the manufacturing method, adhesive quality, and the precision of the curing process. Fig. 6 illustrates a representative point defect and its influence on the microvolume of the honeycomb. The presence of a vacancy distorts the local structure, making the affected microvolume more susceptible to further damage accumulation.
Point defects.
A point defect creates favorable conditions for the rapid propagation of damage and, importantly, determines the direction of delamination spread. Identifying the exact location of such defects is challenging; instead, one must determine probabilistic zones of weakness and assess them statistically. Because failure propagation depends on the distribution and strength of honeycomb nodes – which serve as the primary load-transfer sites – statistical models allow only for estimating the probable damage growth rate and the residual strength of the affected region.
From a materials-science perspective, dislocation motion involves the sequential displacement of atomic half-planes, with bonding across slip planes breaking as dislocations move. This process separates slipped material from unslipped material. Engineering materials rarely achieve their theoretical strength because microflaws create stress concentrations that lead to premature failure; similarly, sharp geometrical features significantly amplify local stresses.
In honeycomb composites, the fabric pattern of the skins and their local thickness variations affect stress distribution and reinforcement behavior. Because the panel as a whole is highly complex, it is impractical to analyze the entire structure directly. Instead, a representative sub-volume of the honeycomb core is extracted and examined as the basis for mechanical and probabilistic modeling. Figure 7 shows an example panel with reinforced sheathing and a corresponding honeycomb fragment selected for analysis.
Core alignment for the panel in sandwich.
In the selected representative volume, we identify and count the honeycomb nodal points that contribute to skin reinforcement. For illustration, a circular region is chosen, and all node points within this region are counted row by row. As long as at least one node point within the region remains intact, the element is considered undamaged. Damage is assumed to occur only when all node points included in the selected core region have failed (Fig. 8).
Numerical values of nodal points in each row of a honeycomb fragment
By performing a simple row-based count, we obtain the following number of nodal points within the analyzed fragment: 3 × 2 + 11 × 2 + 13 + 15 + 17 × 2 + 16 × 2 + 14 = 136 nodal points. (These represent nodes, not cells.)
These values will serve as reference parameters in subsequent calculations related to residual strength, fatigue behavior, and predicted failure time of the composite sample.
To model the initiation and progression of delamination, we use the inverse transform method together with the probabilistic formulation of a Poisson process. Since a computer can generate independent and identically distributed (i.i.d.) random variables that are uniformly distributed on the interval (0, 1), it is essential to transform these uniform random values into samples from any desired probability distribution (in our case, a Poisson distribution).
Let F(r), r ∈ ℝ, denote any cumulative distribution function (CDF) (continuous or not). Recall that F(x): ℝ → [0,1] is thus a non-negative and non-decreasing (monotone) function that is continuous from the right has left hand limits, with values in [0, 1]; moreover F(∞) = 1 and F(-∞) = 0. Our objective is to generate (simulate) a random variable X distributed as F, that is, we want to simulate a random variable X such that P(X ≤ x) = F(x), x ∈ ℝ.
Define the generalized inverse of F, F−1 : [0.1] → ℝ as F−1(y) = min{x:F(x) ≥ y}, y ∈ [0,1]. If F is continuous, then F is invertible (since it is thus continuous and strictly increasing), in which case F−1(y) = min{x : F(x) = y}. The ordinary inverse function is thus F(F−1(y)) = y and F−1(F(x)) = x. In general, it holds that F(F−1(y)) ≥ y and F−1(F(x)) ≤ x and F−1(y) is a non-decreasing (monotone) function in y. This fact yields a simple method for simulating a random value X distributed as F.
For a Poisson distribution with mean λ = 5, the probability mass function is
Note that Y = min{n ≥ 1: tn > 1} = min{n ≥ 1: X1 + … + Xn >1}, which is a stopping time. Using the inverse transform to generate the i.i.d. exponential interarrival times Xi, we can represent Xi = – (1/λ) In( Ui), where Ui are independent uniform random variables on (0,1). We can then re-write expressions using the identity In(xy) = In(x) + ln(y), which allows us to simulate Y by simply generating independent uniforms Ui.
To justify the inverse transform, recall that F(F−1(y)) = y and so by the monotonicity of F, if F−1(U) ≤ x, then U = F(F−1(U)) ≤ F(x), that is, U ≤ F(x). Similarly, F(F−1(x)) = x and so if U ≤ F (x), then F−1(U) ≤ x. Therefore, the two events are equivalent (as was to be shown). In the general case (continuous or not), it is easily shown that {U < F(x)} {F−1(U) ≤ x} {U≤F(x)}, which yields the same result after taking probabilities (since P(U = F(x)) = 0 when U is a continuous random value).
In our example, however, we use a different algorithm. Instead of applying inverse functions, we directly compute probability values using built-in Excel functions. This provides a tabulated dependence of the CDF on its argument. By comparing a generated uniform random value with the tabulated probabilities, we automatically obtain the corresponding inverse discrete output without any explicit transformation. The interval probabilities are defined as
The discrete selection step may be symbolically expressed as
In the numerical example, the time points are Ti = 1, …, 5 and the line Pt = with the values 15, 19, 45, 59, 74. Above this, the corresponding “increments by steps” are 15, 4, 26, 14, 15, whose sum is 74. This is precisely the cumulative sum evaluated at successive time steps.
Let T1 < T2 < … < Tm – a be discrete “stages” (m = 5 in this example), and let ΔP ≥ 0 be the number of new damaged nodes that appeared at step. The cumulative damage function is then
For our example, (ΔP1, ΔP2, ΔP3, ΔP4, ΔP5) = (15, 4, 26, 14, 15). Thus,
If the increments are modeled by independent Poisson variables ΔPi ~ Pois(μi), then at every stage
This behavior is fully consistent with the stochastic formulation adopted in this work, in which the initiation and growth of defects are represented by a Poisson process.
To obtain the cumulative distribution function (CDF) of the damage process, we use the Poisson distribution with mean λ = 5, interpreted as the average rate of cell failures. For illustration, we consider five discrete stages of a conditional time variable, T1, T2, T3, T4, T5. At each stage, damage accumulates in a given honeycomb fragment as a combination of failures across all cell rows. In the numerical implementation (Fig. 9), we construct two columns, X and Y: X contains the discrete damage levels (or time steps), Y contains the corresponding CDF values computed from the Poisson distribution.
Calculation of the Poisson distribution function and mapping to honeycomb nodes in the Excel model.
To calculate the distribution, we implement an inverse comparison procedure. After modelling in the “random match values” column with the Excel formula:
This is a complex array formula, and it is necessary to understand the entire recalculation process in order to use the data correctly. In this procedure, the values of the variables and the resulting data are determined by comparison. Each stage passes through a simulation of random numbers, which are then normalized through the Poisson distribution using the Excel formula:
In order to visualize the spread of delamination along the honeycomb element of the structure, we selected nearby nodes and shaded the adjacent cell areas (see Fig. 10). Although this was not the main task of the solution, it makes it easier to understand the process obtained in the simulation. Physical experiments are needed to accurately select the distribution parameters, and in this work we do not determine the residual strength or the exact time of damage development. However, they can be estimated by dividing the conditional strength of the entire sandwich element by the number of nodal points. As the time parameter, one should take the time of damage initiation and then multiply it by the number of stages obtained from the calculations.
Example visualization of simulated delamination spreading across honeycomb cells.
Modeling a large number of cells, or a larger number of damage events, becomes difficult to perform manually. However, such simulations are of interest when considering 1000 points or more. For this purpose, we used an R package (Fig. 11), which allowed us to perform extended evaluations and comparisons in future work. As for determining the residual strength of the product and the exact time of damage propagation, this task remains open and will be addressed later, possibly in future studies.
Regular initialization estimation program.
The listing of the program in the R package (Fig. 11) demonstrates the complexity of transforming and processing the simulated values. Using this approach, the size of the dataset can be increased substantially, enabling deeper analysis. However, in such large-scale simulations, visual representation becomes more difficult compared to simple examples. Fig. 12 shows the program working with 1000 points within a honeycomb structure; creating such a representation manually would be impractical and extremely time-consuming. Therefore, these stages need to be generated programmatically, which makes it possible to track and analyze the dynamics of damage development.
Universal propagation time.
This study has proposed that delamination in composite and honeycomb sandwich structures follows a damage-propagation mechanism fundamentally different from those observed in metals or in homogeneous polymeric materials. In the systems considered, the absence of significant plastic deformation and the dominant role of nodal attachment points result in a brittle, avalanche-type progression of damage. For this reason, a Poisson process was adopted to describe the initiation and cumulative growth of weak microvolumes.
The modelling framework developed here demonstrates that Poisson-driven defect nucleation can capture essential features of delamination spread and provides a basis for estimating cumulative damage and residual integrity. The numerical examples, implemented both in Excel and in R, illustrate how probabilistic methods can support the analysis of subsurface damage that cannot be directly observed in honeycomb structures.
However, several aspects of the model require further development. A more accurate selection of distribution parameters will necessitate controlled physical experiments. Additionally, the relationship between defect nucleation rate, structural geometry, and service loads should be examined in detail. Large-scale simulations, as shown in this study, are feasible but require fully programmatic approaches for meaningful visualization and interpretation.
Overall, this paper is intended as a preliminary framework that organizes existing concepts and demonstrates how stochastic modelling can be applied to composite delamination. The authors hope that this approach will support future research on damage evolution and residual strength assessment in advanced composite structures.