This paper presents advanced computational techniques for predicting deformation and fracture mechanisms at the microscale in aluminium-silicon composite materials. Emphasising software engineering and computational modelling, finite element simulations have been carried out to analyse how inclusion size (100, 50, and 20 μm) and deformation velocity (1, 2, and 4 mm/s) influence micro-scale fracture phenomena. The aluminium matrix behaviour is described using a piecewise linear plasticity model, while the brittle silicon inclusions are characterised via the Johnson-Holmquist damage model. Simulation results indicate that larger inclusions fracture early and fragment significantly, causing pronounced damage propagation into the matrix at lower strain levels (~7 %). Medium-sized inclusions exhibit delayed fracture (~19-20 % strain) but still critically affect matrix integrity. The smallest inclusions remain intact under high strains but induce matrix failure through stress concentration and rotational effects. Slower deformation velocities significantly delay matrix failure initiation. Computational outcomes are validated by experimental data, highlighting the critical role of inclusion geometry and the advantages of spheroidisation in enhancing composite ductility. The simulation system incorporates modular software engineering methodologies, facilitating effective integration and reuse of computational assets, enhancing maintainability and adaptability. Future developments envisage incorporating ontological frameworks and artificial intelligence techniques to improve semantic interoperability and data-driven analysis capabilities.
Computational modelling has become a crucial component in modern engineering and materials
science, significantly reducing the need for costly and time-consuming experimental investigations.
Among various numerical methods, finite element analysis (FEA) stands out due to its flexibility,
accuracy, and wide applicability in solving complex mechanical problems. Most of the existing
computational studies have primarily addressed macroscopic behaviours of structures and
materials; however, micro-scale simulations, essential for heterogeneous or composite materials,
still pose considerable methodological and computational challenges [
In particular, metal matrix composites (MMCs), consisting of ductile matrices reinforced by
brittle inclusions, present complex microstructural interactions under mechanical loading.
Accurate prediction of deformation and fracture phenomena in these composites necessitates
advanced computational models and robust numerical algorithms. Prior research often simplified
the inclusion-matrix interaction, overlooking critical aspects such as inclusion fracture mechanics,
shape irregularities, and strain-rate dependence [
This study addresses these gaps by developing and implementing advanced numerical models
and computational algorithms for micro-scale finite element simulations of aluminium matrix
composites reinforced by silicon inclusions. Emphasis is placed on the software engineering
aspects of simulation, including algorithmic efficiency, stability of numerical integration schemes,
and computational validation procedures. Specifically, the Johnson-Holmquist damage model for
brittle materials and a piecewise linear plasticity model for ductile matrices have been integrated
within a comprehensive finite element framework [
The primary objectives of this research are twofold: firstly, to improve the predictive accuracy of micro-scale deformation and fracture behaviours in MMCs, and secondly, to refine the computational methodology, ensuring its suitability for high-performance computing environments. The outcomes of this work will provide valuable insights for further algorithmic improvements and practical applications in materials design, contributing significantly to the intersection of computational mechanics and software engineering.
The finite element method (FEM) serves as the core computational approach employed in this study to simulate the deformation and failure mechanisms within metal matrix composites (MMCs). FEM allows the spatial discretization of complex geometries and the solution of governing equations under prescribed boundary and initial conditions. It has been extensively applied to problems in structural mechanics, including impact dynamics, fracture propagation, and plastic deformation. A general scheme of the modelling technique with application to the considered problem of deformation and fracture simulation is presented in Figure 1.
In this work, 8-point hexahedral solid elements with reduced integration are used to discretize both the ductile aluminium matrix and the brittle silicon inclusions. This choice ensures computational efficiency while maintaining numerical accuracy in capturing stress gradients and localized damage. Stability of the FEM simulation is enhanced using a second-order objective stress update algorithm and an hourglass control mechanism based on the Flanagan-Belytschko formulation.
To accurately capture the differing mechanical responses of the two phases within the composite, material-specific constitutive models are embedded into the FEM formulation. The aluminium matrix is modelled using a piecewise linear plasticity approach, which represents the nonlinear stress-strain response by segmenting it into linear intervals. This method is computationally efficient and sufficiently flexible to represent yield plateau, strain hardening, and failure thresholds. For the brittle silicon inclusions, the Johnson-Holmquist (JH-2) damage model is implemented. This model is particularly suitable for high-strength brittle materials and accounts for both elastic and inelastic responses, as well as progressive damage evolution under compression and tension. It utilizes a scalar damage variable that accumulates over time, resulting in degradation of material stiffness and strength.
An essential aspect of the FEM setup is the precise definition and calibration of the material parameters within these models. The accuracy and stability of the numerical simulation heavily depend on these parameters, which include elastic moduli, plastic yield stress, fracture strain, and damage evolution constants. These values are derived from experimental data and prior literature, and their proper selection ensures the physical relevance of the simulation results [20].
The general structure of the computational system aimed for the considered simulation is presented in Figure 2 as a C4 diagram.
The simulation system developed for finite element analysis of deformation and fracture processes is architected as a modular and component-oriented application. Its structure is designed to ensure clarity of responsibilities across participants, modular reusability, and adaptability for multiscale physical modelling. The overall workflow begins with the Model Geometry Designer, a domain expert responsible for creating 2D or 3D geometrical representations of the object or system under investigation. These geometrical models are typically constructed using external CAD or modelling platforms, which export geometry data in formats such as .igs, .stp, or .stl. This exported geometry serves as an input to the simulation environment.
A computer modelling specialist, often a materials scientist for the considered subject area, defines the physical conditions, boundary constraints, and simulation parameters. These settings are specified via a Pre-processor component, which provides a graphical interface for assigning initial and boundary conditions, associating constituent materials, and managing interaction models among composite phases. The pre-processor produces a task descriptor (typically a structured keyword-based text file) that encapsulates all simulation metadata.
At the heart of the simulation workflow is the Solver, implemented as a compiled application. The solver performs numerical integration and deformation/fracture simulations using finite element methods. It interacts directly with the Materials Models Base, a container that holds calibrated constitutive models – such as the Johnson-Holmquist model for brittle phases and the piecewise linear plasticity model for ductile components. The material definitions are provided in a special internal format and loaded dynamically by the pre-processor and solver subsystems. The solver also communicates with the Computational Models module, a container housing compiled programmatic components that define numerical algorithms, storage schemes, and computation logic. These models ensure the modular execution of solvers and enable scalability and highperformance computation.
Upon task completion, simulation results are returned in the form of structured text output files, which are consumed by the Postprocessor. The postprocessor is a visualisation-oriented software tool that provides users with real-time access to stress-strain fields, damage distributions, crack initiation zones, and other key physical indicators. It enables detailed analysis of the simulation output and supports decision-making in materials evaluation and design.
In summary, the system is a robust, extensible platform tailored for high-fidelity simulation of composite material behaviour, supporting the research needs of computational mechanics, materials engineering, and fracture analysis.
where: M is the matrix domain (e.g., aluminium, ductile); I is the inclusion domain (e.g., silicon, brittle); G is the geometry of the composite, a spatial configuration G ⊂ R 3 ; P ={PM , PI } is he set of material parameters for M and I ; F is the space of applied forces and boundary conditions; N is the finite element discretisation (mesh); A is the set of numerical algorithms (e.g., integration scheme, damage evolution law).
Let the governing PDE system be compactly represented as a nonlinear operator: where: U is a sutable function space, e.g. U ⊂ H 1(Ω)d (Sobolev space for displacements); U * is its dual; u(x, t) is the known displacement field; f ∈ U * is the external force distribution.
The weak form (variational form) becomes:
L :: U → U * , such that L :(u) = f ∫ σ(u) : ∇δudx = ∫ f ⋅ δudx Ω Ω ∀δu ∈ U
Let the composite material system be defined as a structured tuple:
S =
M , I ,G, P , F , N , A (1) (2) (3) (4) (5) (6) (7) (8) (9) where: Φ includes: meshing, element-wise assembly, integration, and solution; R is the space of output fields:
Φ : S → R R = {u(x, t), σ(x, t), ε(x, t), D(x, t),...} Φ (S ) = Solve[L M (u) + L I (u) = f ]
N = {E j },{xi} ,
E j ⊂ Ω with L M and L I
denoting operators induced by their respective constitutive models on domains ΩM , ΩI ⊂ Ω.
Let the computational mesh be defined as:
Each finite element Ej is assigned local material parameters and evaluated using quadrature. Time evolution is performed using a discrete time integrator (explicit or implicit): where the stress field σ depends on the constitutive laws of M and I .
The materials models could be formally presented as constitutive functionals. Let σM = CM (ε) – for matrix, with piecewise linear plasticity:
CM : ε a σM , via PLP(ε) = U(Ei ⋅ ε) i σI = CI (ε, D) – for inclusions, via Johnson–Holmquist model:
CI : (ε, D) a σ ,
I
with D = D (ε&, p) Each Ck is a nonlinear mapping depending on strain ε, strain rate ε&, pressure p, and damage D. Let us define the full FEM simulation as a mapping: un+1 = un + Δt ⋅ u&n + ... (10)
The presented formalization enables modular implementation (e.g., FEM solver with interchangeable models), facilitates rigorous comparison of constitutive models via functionals, allows for mathematical optimisation, sensitivity analysis, and adjoint methods, lays foundation for uncertainty quantification and multi-scale coupling.
The general used here simulation technique is the finite elements method. The 8-point hexahedron
solid element type was used both for the inclusion and matrix geometry mesh, which possesses a
specific geometry and comprises eight integration points [
Second-order objective stress update [
The hourglass energy was computed and included in the energy balance. Hourglass modes are
element distortions that have zero strain energy. Thus, no stresses are created within the element
[
“Stonewall”, Rayleigh and sliding energy dissipations were computed and also included in the energy balance. Rayleigh energy dissipation refers to the energy loss in a dynamic system due to Rayleigh damping, a common damping model used in finite element and structural dynamics analyses. In this model, the damping forces are represented as a combination of mass- and stiffness-proportional contributions. Initial reference geometry energy was computed and included in the energy balance as part of the internal energy.
The Newmark time integration constants have the values β = 0.5 and γ = 0.25. The Newmark
time integration method is a widely used numerical scheme for solving dynamic problems in
structural mechanics [
The minimal time step was set to be 1.0×10-6 s and initial time step 1.0×10-5 s.
For the aluminium matrix simulation a piecewise linear plasticity model [
Compression:
ρ µ = −1
ρ 0
P = K1 ⋅µ P = K1 ⋅µ + K2 ⋅µ 2 + K3 ⋅µ 3 + ΔPn−1 120 100 aP80 M ,s 60 s tre 40 S 20 0
0 Density: ρ = 2712.6 kg/m3; Young's modulus: E = 68.948 GPa; Poisson's ratio: ν = 0.33; Yield stress: σY = 15.21 MPa;
Plastic strain to failure: ε = 0.7496
The stress strain curve used for the modelling as a table of values from the yield stress to the failure strain and stress is shown in Figure 3.
To describe the material behaviour of the silicon inclusion, which unlike the matrix is brittle,
the Johnson-Holmquist damage model [
Deformation: (11) (12) (13) (14) 20 40
Plastic strain, % 60 80
In equation (12) P corresponds to the bulking pressure of the material and is determined by the amount of accumulated damage.
The model introduces a scalar damage parameter D that ranges from 0 (undamaged, intact material) to 1 (fully damaged material). As D increases, the effective strength of the material degrades, reflecting the accumulation of micro-cracks and other damage mechanisms. The effective stress of the material is expressed as an interpolation between the intact (undamaged) material strength and the fully damaged (residual) strength (14):
σ * = (1 − D) ⋅σ i* + D ⋅σ d* where: σ i* is the intact strength; σ d* is the residual (fully damaged) strength.
This formulation ensures that as damage accumulates, the material response transitions smoothly from the higher, intact strength to a lower, damaged state.
The intact material strength is defined as:
σ i* = A ⋅ (P * +T *)N ⋅ (1+ C ⋅ lnε&)
σ *f = B ⋅ (P*)M ⋅ (1+ C ⋅ lnε&)
The equation (4) is used with the radial return method to determine the current increment in plastic strain. The current increment in damage can be determined as follows: (15) (16) (17) (18) ΔD = Δε P
ε f ε f = D1 ⋅ (P * +T*)D2 The plastic strain to fracture under a constant pressure is defined by the equation (18): The bulking pressure is zero for undamaged material, and the bulking pressure at the next time increment is given by the equation (19): ΔPn+1 = −K1 ⋅µ + (K1 ⋅µ + ΔP)2 + 2 ⋅β ⋅ K1 ⋅ ΔU (19) pairs of features, and the one that resulted in the greatest accuracy improvement was added to the model. The search continued until one of the following conditions was met: no further accuracy improvement, desired accuracy was reached, the feature set reached a predefined size limit, or no new features were available (which could occur in cases with very few extra features).
In the final stage, different classification methods were applied to the selected feature set to optimize the results. This process led to the development of a complex classifier capable of distinguishing between different class groups, significantly improving the ability to identify multiple classes within the dataset.
The surface-to-surface contact model has been used to describe the impacts between the inclusion and matrix, the parts of the inclusion after its fracture and the pasts of the matrix in the case of their contact at later deformation stages. The friction coefficients for the considered contact pairs are presented in Table 1.
In the inclusion / matrix contact pair inclusion has been treated as a “master” and the matrix as a “slave” segment.
The presented here study focuses on the cold deformation simulation (being carried out at room temperature) not considering hot deformation processing that assumes considerably different materials’ properties.
A single inclusion has been considered during each of the computation experiments. So no intersection between the inclusions has been taken into account (except the split parts of the fractured initial inclusion).
The inclusions had the shape of flat parallelepipeds. The deformation was performed perpendicular to the inclusion plane. Such conditions are primarily aimed at simulating silicon crystals in the cast Al-Si alloys. The thickness of the considered inclusions was 10 μm and the lengths (equal to width) were 100, 50, and 20 μm. The deformation velocity was a constant value given as a prescribed motion of the top and the bottom nodes of the matrix block. For the study of the inclusion size influence it was set at 2 mm/s (1 mm/s from top and 1 mm/s from bottom – towards one another). For the inclusion of 100 μm also the influence of deformation velocity has been studied using the values 1, 2 and 4 mm/s. The matrix block inside which the inclusion model was implemented had dimensions: height 100 μm, length (equal to width) 200 μm. Figure 4 presents a schematic example of the initial model geometry.
The simulation results reveal distinct differences in fracture mechanisms associated with the
inclusion sizes studied. A 100 μm inclusion fractures early in the deformation process, initially
splitting into two primary fragments around its midpoint. With further deformation, particularly
beyond approximately 18% settling, these primary fragments further disintegrate into smaller,
clustered fragments. Eventually, fragmentation continues predominantly near the primary fracture
zone, and the separated fragments undergo rotational and translational motion within the matrix.
This behaviour may account for experimentally observed arrangements of fragmented silicon
inclusions, including linear aggregations and "bag"-like clusters described previously [
In contrast, the 50 μm inclusion exhibits significantly lower susceptibility to early fracture. The first cracks in this intermediate-sized inclusion occur much later, at approximately 19–20% deformation compared to the 7% deformation threshold observed for the 100 μm inclusion. Fracture fragments for the 50 μm inclusion predominantly range between 25–30 μm, aligning closely with the computationally identified critical size for relatively stable inclusions under the investigated deformation conditions. Notably, fragmentation at this scale still critically impacts matrix integrity, as observed tearing is consistently initiated by rotational and translational stresses imposed by these smaller fragments.
The smallest inclusion considered (20 μm) demonstrates the highest fracture resistance in simulations. Brittle fracture is not initiated within this inclusion; however, significant deformationinduced rotations and movements occur at high deformation levels (above 26–27% setting), eventually causing severe stress concentrations and tearing in the aluminium matrix. This behaviour suggests that inclusions smaller than the critical size identified (approximately 25–30 μm) can still significantly degrade matrix performance through mechanical interactions despite not fracturing themselves.
The numerical outcomes align well with experimental observations documented in previous
studies [
From the Figure 6 we can see that deformation velocity modulates the onset of matrix failure, with slower rates enhancing ductility by delaying crack initiation.
Additionally, the computational framework presented here benefits from advanced software engineering methodologies, such as the systematic class diagram design approach described by Chebanyuk [22], which enhances software modularity and maintainability. Furthermore, adopting mentioned there methods for analysing software requirements and artifact reuse through artificial intelligence technologies [23, 24] can significantly increase the robustness and adaptability of simulation software, allowing for more efficient updates and expansions of computational models.
Future developments of the simulation environment could also incorporate advanced text-tomodel transformations [25] and multilingual question-driven approaches [26], facilitating broader accessibility, more efficient user interactions, and enhanced semantic interoperability in computational materials science.
Sdeetgtrieneg,% 1 mm/s 2 mm/s 4 mm/s
The numerical simulations performed in this study highlight the effectiveness and robustness of advanced computational models and finite element methods for predicting deformation and fracture behaviours of aluminium-silicon composite materials. Results confirm that inclusion size and deformation velocity significantly impact fracture mechanisms, demonstrating the sensitivity of numerical models to microstructural parameters. Larger silicon inclusions exhibit early brittle fracture and contribute to the propagation of cracks into the aluminium matrix, whereas smaller inclusions predominantly lead to matrix deformation through stress concentrations and geometric interactions.
A critical threshold inclusion size (~25–30 μm) was computationally determined, below which brittle fracture in silicon inclusions becomes negligible, yet structural integrity remains compromised due to increased stress localisation. Computational results align closely with experimental observations in Al-Si composites, validating the simulations' accuracy and highlighting the detrimental influence of inclusion geometry, particularly sharp-edged inclusions. Furthermore, the beneficial effects of morphological optimisation, such as spheroidisation, are confirmed by simulations, suggesting practical avenues for enhancing composite material performance.
From a software engineering perspective, this study underscores the importance of precise numerical algorithm implementation and efficient computational methodologies for accurate micro-scale simulations. The explicit integration algorithms and damage evolution models, such as the Johnson-Holmquist model and the piecewise linear plasticity model, demonstrated computational stability and predictive accuracy, establishing the potential of these tools for guiding material design and processing decisions.
The outcomes of this research emphasise the role of computational modelling as a powerful predictive tool for exploring microstructural deformation and fracture phenomena in metal matrix composites. Future work will extend these computational approaches by developing an ontological framework [27, 28, 29] that systematically captures and structures the simulation data. Such an ontology aims to facilitate deeper integration between computational mechanics and materials science, enhancing semantic interoperability [27, 28, 29], promoting data reusability, and supporting advanced computational analyses to uncover novel insights into microstructural influences on mechanical performance.
The research was conducted under project number 0225U002527 (2024) on "Development of a technology for manufacturing powder composite materials of the Al-Si-Ni system for highly stable device elements" [27, 30], based at the I. M. Frantsevich Institute for Problems in Material Sciences of the National Academy of Sciences of Ukraine.
The author(s) have not employed any Generative AI tools.
Vladislav Kaverinsky: Investigation, Formal analysis, Writing – Original Draft, Visualization, Software, Methodology. Zoya Sukhenko: Conceptualization, Resources, Validation, Writing – Review & Editing, Supervision Anna Litvin: Writing – Review & Editing. Sergii Kotlyk: Writing – Review & Editing. Kyrylo Malakhov: Validation, Supervision, Writing – Review & Editing. [20] A. Ruggiero1, G. Iannitti, N. Bonora, M. Ferraro, Determination of Johnson-Holmquist constitutive model parameters for fused silica, DYMAT 2012 - 10th International Conference on the Mechanical and Physical Behaviour of Materials under Dynamic Loading. 26 (2012). doi: 10.1051/epjconf/20122604011. [21] V.V. Kaverinsky, G.A. Bagliuk, A.D. Rud et al., Al–Si–Ni Sintered Alloy: Structure and
Properties. I. Powder Production, Metallofiz. Noveishie Tekhnol. 45, 1039–1050 (2023). [22] Chebanyuk, O. An approach to class diagram design. in Proceedings of the 2nd International Conference on Model-Driven Engineering and Software Development (MODELSWARD 2014), Lisbon, Portugal, 2014. Scitepress. doi:10.5220/0004763504480453. [23] Chebanyuk, O. Requirement Analysis Approach to Estimate the Possibility of Software Development Artifacts Reusing Consulting with Artificial Intelligence Technologies. CEUR Workshop Proceedings. 2024. Vol. 3806. pp. 62-74. https://ceur-ws.org/Vol3806/S_43_Chebanyuk.pdf. [24] Chebanyuk, O. An Approach to Software Assets Reusing. Lecture Notes of the Institute for Computer Sciences, Social-Informatics and Telecommunications Engineering. 2022. LNICST 450. pp. 73-83. https://doi.org/10.1007/978-3-031-17292-2_6. [25] An approach of text to model transformation of software models. ENASE 2018: proceedings of the 13th International Conference on Evaluation of Novel Approaches to Software Engineering (Funchal, Madeira, March 23–24, 2018). Funchal, Madeira (Portugal), 2018. pp. 432-439. https://doi.org/10.5220/0006804504320439. [26] Chebanyuk, O. Multilingual Question-Driven Approach and Software System to Obtaining Information From Texts. CEUR Workshop Proceedings. 2022. Vol. 3501. pp. 256-265. https://ceur-ws.org/Vol-3501/s24.pdf. [27] M. G. Petrenko, E. Cohn, O. Shchurov, K. S. Malakhov, Ontology-driven computer systems: Elementary senses in domain knowledge processing, South African Computer Journal 35 (2023) 127–144. doi:10.18489/sacj.v35i2.17445. [28] O. Palagin, K. S. Malakhov, V. Y. Velychko, T. Semykopna, O. Shchurov, Hospital information smart system for hybrid e-rehabilitation, in: CEUR Workshop Proceedings, volume 3501, CEUR-WS, 2022, pp. 140–157. URL: https://ceur-ws.org/Vol-3501/s50.pdf. [29] O.V. Palagin, V.Yu. Velychko, K.S. Malakhov, O.S. Shchurov, Distributional Semantic Modeling: a Revised Technique to Train Term/Word Vector Space Models Applying the Ontology-related Approach, in: CEUR Workshop Proceedings, volume 2866, CEUR-WS, p. 342–353. URL: http://ceur-ws.org/Vol-2866/ceur_342-352palagin34.pdf. [30] V.V. Kaverinsky, G.A. Bagliuk, Z.P. Sukhenko, Numerical Simulation of In Situ Reaction Synthesis of TiC Reinforced Aluminum Matrix Composite from Elemental Al-Ti-C Powders, J. Mater. Eng. Perform. 33, (2024) 9976–9986. doi:10.1007/s11665-023-08650-6.