PSI - Issue 79

A. Della Rocca et al. / Procedia Structural Integrity 79 (2026) 475–484

477

The work is structured as follows: beside a short introduction to depict the motivation of the study, the methodological approach explains the spinodal structure generation (section 2.1), the morphological and topological analysis (section 2.2), the finite element modeling (section 2.3) and the statistical analysis implemented to predict the effective stiffness. Preliminary results in terms of 2D and 3D morphological analysis are reported in section 3.1 and 3.2 and used for the implementation of the FEA, the evaluation of the effective stiffness and the predictive capability via linear regression modeling (section 3.3). The results obtained are then discussed in the light of the limitations encountered in Section 4 with the main conclusions and future work highlighted as well. 2. Materials and Methods 2.1. Spinodal structure generation and parametrization Spinodal structures were generated by numerically solving the Cahn–Hilliard equation, which models phases separation in a binary system as it evolves from an initially homogeneous mixture into two continuous and interpenetrating phases Cahn and Hilliard (1958). This process—known as spinodal decomposition —produces bicontinuous and isotropic morphologies that resemble the trabecular, lattice-like complexity of cancellous bone. In this formulation, the temporal evolution of the concentration field c(x,t), representing the local composition of the two phases, is governed by: డ డ ௖ ௧ ൌ ߘ ή ൬ ߘܦ ቀ డ௙ డ ሺ ௖ ௖ሻ െ ߘߛ ଶܿ ቁ൰ (1) where D is the diffusivity , controlling the dominant wavelength (or separation scale) of the emerging porous features, and γ is the gradient energy coefficient , which modulates the sharpness of the solid–void interfaces. The derivative of the free energy density, ∂ f(c)/ ∂ c, determines the thermodynamic driving force of phase separation. A MATLAB implementation of the Cahn–Hilliard equation was used, employing spectral methods to minimize numerical instabilities and efficiently simulate the dynamics of phase separation. The initial condition consisted of a random perturbation around an average concentration field (c=0.5±0.1), ensuring that the structures developed from stochastic but statistically controlled fluctuations rather than random noise. A parametric study was performed by systematically varying D and γ between 1 and 19 (step size = 2). Each parameter combination generated a unique spinodal composition, resulting in a total of 100 distinct structures for both 2D and 3D configurations (visible in Fig. 1). This approach provided a wide morphological spectrum—from fine to coarse porosities—while maintaining reproducibility and physical coherence. Values above 19 were avoided, as excessively high D and γ values caused two main issues: (i) structural degeneration and loss of isotropy due to overdeveloped features relative to the sample size, and (ii) numerical instability leading to computational crashes. 2.2. 2D and 3D Structure Generation Two-dimensional spinodal structures were first generated to explore the influence of the parameters D and γ on morphology and topology. Each simulation produced isotropic, bicontinuous patterns exhibiting interconnected pore and solid (respectively white and black/red areas in Fig. 2) phases. These 2D models were subsequently used for morphological and topological characterization, serving as a computationally efficient pre-screening stage before full 3D analysis. Three-dimensional spinodal structures were generated within a cubic domain of 64 × 64 × 64 voxels (pixels), corresponding to a physical size of 6.4 mm³ (voxel size = 0.1 mm). This domain size was chosen to balance structural representativeness and computational efficiency—large enough to capture characteristic spinodal features without excessive computational cost. All structures were oriented with the Y-axis aligned vertically, corresponding to the loading direction used in subsequent finite element analyses. Following phase separation, the evolved concentration field was binarized at a threshold c th =0.5, distinguishing solid and void regions. The porosity was then adjusted using iterative erosion–dilation operations to reach a target bone volume fraction (BV/TV) of