The strength of flax fiber-reinforced polymer (FFRP) composites is a critical factor influencing their reliability in structural applications. However, the prediction of FFRP strength remains challenging due to the nonlinear stress–strain relationship of flax fibers. To address this limitation, this study proposes a three-dimensional progressive failure model (PFM) based on the fiber bundle chain method, with a specific focus on investigating the influence of size effects on the model's predictive performance. The PFM was constructed by generating a representative volume element (RVE) in which each fiber element is assigned an independent tensile strength, \({\sigma }_{\text{p,q}}^{\text{u}}\) , using Weibull statistical theory, enabling the model to simulate progressive damage in unidirectional composites. To enhance the tensile performance of FFRP, flax fibers were surface-modified via grafting with multi-walled carbon nanotubes (MWCNTs). The PFM was subsequently employed to numerically simulate the tensile behavior of unmodified FFRP, MWCNT-grafted flax fiber composites, carbon–flax hybrid composites, and basalt–flax hybrid composites. The results demonstrate that for the MWCNT-grafted flax fiber composites, the PFM coupled with Zhou's stress concentration factor (SCF) calculation model predicted a tensile failure strain of 1.26% and a tensile strength of 406.8 MPa. Compared with experimental data, the relative errors were 5.9% and 10.3%, respectively, indicating that the PFM can effectively predict the tensile stress and strain of the aforementioned composite systems. Furthermore, the model elucidates the influence of size effects: as the cross-sectional dimensions of the RVE increase, or as the total length of the RVE and the length of fiber elements increase, the tensile strain and stress predicted by the PFM exhibit a decreasing trend. Therefore, the proposed PFM is suitable for simulating the tensile properties of a wide range of composite materials.
