On the Influence of Inhomogeneous Interphase Layers on Instabilities in Hyperelastic Composites

Polymer-based three-dimensional (3D) printing—such as the UV-assisted layer-by-layer polymerization technique—enables fabrication of deformable microstructured materials with pre-designed properties. However, the properties of such materials require careful characterization. Thus, for example, in the polymerization process, a new interphase zone is formed at the boundary between two constituents. This article presents a study of the interphasial transition zone effect on the elastic instability phenomenon in hyperelastic layered composites. In this study, three different types of the shear modulus distribution through the thickness of the interphasial layer were considered. Numerical Bloch-Floquet analysis was employed, superimposed on finite deformations to detect the onset of instabilities and the associated critical wavelength. Significant changes in the buckling behavior of the composites were observed because of the existence of the interphasial inhomogeneous layers. Interphase properties influence the onset of instabilities and the buckling patterns. Numerical simulations showed that interlayer inhomogeneity may result in higher stability of composites with respect to classical layup constructions of identical shear stiffness. Moreover, we found that the critical wavelength of the buckling mode can be regulated by the inhomogeneous interphase properties. Finally, a qualitative illustration of the effect is presented for 3D-printed deformable composites with varying thickness of the stiff phase.


Introduction
Composite materials are deeply integrated in the modern life, due to their excellent mechanical and functional properties, which until recently were unachievable. Composites containing two or more constituents can be tailored to meet specific requirements by design of their geometry and smart selection of constituent materials. One of the most challenging problems in composite science is associated with the prediction of their failure. While delamination was historically considered to be the main failure mode, the phenomenon of local buckling or loss of stability attracted significant attention recently. In contrast with delamination, when composites lose the integrity catastrophically, the loss of stability can be considered as reversible microstructure transformation mode; once the external loads are removed, the initial undeformed state is restored thanks to the stored elastic energy. This phenomenon of instability-induced microstructure transformation has been employed to design materials with the switchable properties and functions [1][2][3][4][5][6][7][8].
The pioneering works of Rosen [9] and Hill and Hutchinson [10] laid the basis for the theoretical understanding of the elastic instability phenomenon. The development of analytical

Numerical Simulations
To analyze the mechanics of the layered composite with interphasial layers bonding two purely homogenous hyperelastic constituent materials, we utilized the finite element method with the help of COMSOL Multiphysics (v. 5.2, COMSOL AB, Stockholm, Sweden). Representative volume element (RVE) with dimensions t = a, h = 0.05a was used in the calculations, as shown in Figure 1. A mesh sensitivity analysis was performed, and the RVE containing 2498 quadrilateral elements with quadratic shape functions was used. Figure 1 shows the unit cell of an "ideal" composite without any interphasial layers (a), and the realistic composite, which contains the transition zone between two main constituents (b). In the ideal composite, the thickness of the stiff layer is t l = v l t, where v l is the volume fraction of the stiffer homogenous layer, and t is the period of the unit cell. We assumed that the interphasial layers, generated from the 3D-printing process (see Experimental section), contain equal amounts of the stiffer layer and soft matrix materials; the thicknesses of the pure stiff and the interphasial layers in this case are equal to t l = (1 − f )v l t and t i = f v l t, respectively. Here, 0 ≤ f ≤ 1 is the relative thickness of the interphasial layer. A mesh sensitivity analysis was performed, and the RVE containing 2498 quadrilateral elements with quadratic shape functions was used. Figure 1 shows the unit cell of an "ideal" composite without any interphasial layers (a), and the realistic composite, which contains the transition zone between two main constituents (b). In the ideal composite, the thickness of the stiff layer is = , where is the volume fraction of the stiffer homogenous layer, and is the period of the unit cell. We assumed that the interphasial layers, generated from the 3D-printing process (see Experimental section), contain equal amounts of the stiffer layer and soft matrix materials; the thicknesses of the pure stiff and the interphasial layers in this case are equal to = (1 − ) and = , respectively. Here, 0 ≤ ≤ 1 is the relative thickness of the interphasial layer. Since the interphasial layer is inhomogeneous, we need to define the variation of the local elastic modulus through the thickness of the interface. In this study, we considered three different types of shear moduli distributions through thickness, denoted as A, B and C on Figure 2. Shear modulus in the interphasial zone at position (see Figure 2) is defined as: where and denote the initial shear modulus of soft matrix and stiff layer, respectively. In expression (1), varies from zero to one (see Figure 2), ℎ = ( − ) where = 0.25, 0.5 and 0.75 for distributions A, B and C, respectively. Note, that for distribution B, the elastic modulus linearly increases between the values of the initial shear modulus of the soft matrix and stiffer layer . The distribution A corresponds to the case, when the average elastic modulus of the interphasial layer was lower than ( + )/2, as opposed to the distribution C, when the average shear modulus of the interphasial layer exceeded the value of ( + )/2. All constituents, including the interphasial layers, were considered as nearly incompressible hyperelastic materials with the neo-Hookean strain-energy function, integrated in COMSOL as: where Λ is the first Lame constant, = is the right Cauchy-Green tensor and = det ( ) is the determinant of the deformation gradient . To maintain the nearly incompressible behavior of the constituents, we set a high ratio between the shear modulus and the first Lame constant; in particular Λ = 1000 was used in our simulations. Since the interphasial layer is inhomogeneous, we need to define the variation of the local elastic modulus through the thickness of the interface. In this study, we considered three different types of shear moduli distributions through thickness, denoted as A, B and C on Figure 2. Shear modulus in the interphasial zone µ i at position x (see Figure 2) is defined as: where µ m and µ l denote the initial shear modulus of soft matrix and stiff layer, respectively. In expression (1), x varies from zero to one (see Figure 2), µ h = c(µ l − µ m ) where c = 0.25, 0.5 and 0.75 for distributions A, B and C, respectively. Note, that for distribution B, the elastic modulus linearly increases between the values of the initial shear modulus of the soft matrix µ m and stiffer layer µ l . The distribution A corresponds to the case, when the average elastic modulus of the interphasial layer was lower than (µ l + µ m )/2, as opposed to the distribution C, when the average shear modulus of the interphasial layer exceeded the value of (µ l + µ m )/2. All constituents, including the interphasial layers, were considered as nearly incompressible hyperelastic materials with the neo-Hookean strain-energy function, integrated in COMSOL as: where Λ is the first Lame constant, C = F T F is the right Cauchy-Green tensor and J = det(F) is the determinant of the deformation gradient F. To maintain the nearly incompressible behavior of the constituents, we set a high ratio between the shear modulus and the first Lame constant; in particular Λ = 1000µ was used in our simulations.  In order to analyze the stability of the considered layered composites, we employed the Bloch-Floquet analysis [21] superimposed on finite deformations. The procedure for identifying the onset of instabilities and associated wavenumbers was performed in two steps. First, the unit cell underwent static finite deformation, defined by means of the displacement periodic boundary conditions: Here, the indexes left, right, top and bottom denote the edges of the unit cell (see Figure 1); A and B correspond to the nodes located at the corners. The stiffness matrix, obtained during the solution of the problem (2), is stored for further solution of the eigenvalue problem on the next step, when the Bloch-wave conditions are superimposed on the deformed state by using the boundary conditions: Here, ̃ is normalized wavenumber corresponding to the Bloch wave vector. We swept the values of ̃ in the range from 0 to 10 with a step of 0.05, solving the corresponding eigenvalue problem, until the lowest eigenvalue became zero. If for the considered range of ̃ values only positive eigenvalues appear, the instability was not detected, and the composite remained stable for given . In this case, we increased the compressive strain and repeated the procedure, described above, for an increased , until the zero eigenvalue for non-zero ̃ was found. The strain and wavenumber ̃, for which the first zero eigenvalue was observed, corresponded to the buckling strain and critical wavenumber ̃, defining the buckling shape. Recall that the special case of ̃→ 0 corresponds to the macroscopic instability (long-wave mode); otherwise, the composite undergoes a microscopic loss of stability, developing finite size wavy shapes upon achieving the critical level of deformation [19,20]. Figure 3 illustrates the described numerical procedure, showing a typical evolution of the dispersion curves with applied deformation for dilute composite with = 0.025 (a) and non-dilute composite with = 0.2 (b). Recall that there is an intrinsic connection between the shear wave propagation and elastic instabilities in periodic composites [24]. The continuous black curve in Figure 3a,b describes the dispersion relation in the undeformed state ( = 0); the curves In order to analyze the stability of the considered layered composites, we employed the Bloch-Floquet analysis [21] superimposed on finite deformations. The procedure for identifying the onset of instabilities and associated wavenumbers was performed in two steps. First, the unit cell underwent static finite deformation, defined by means of the displacement periodic boundary conditions: Here, the indexes left, right, top and bottom denote the edges of the unit cell (see Figure 1); A and B correspond to the nodes located at the corners. The stiffness matrix, obtained during the solution of the problem (2), is stored for further solution of the eigenvalue problem on the next step, when the Bloch-wave conditions are superimposed on the deformed state by using the boundary conditions: Here, k is normalized wavenumber corresponding to the Bloch wave vector. We swept the values of k in the range from 0 to 10 with a step of 0.05, solving the corresponding eigenvalue problem, until the lowest eigenvalue became zero. If for the considered range of k values only positive eigenvalues appear, the instability was not detected, and the composite remained stable for given ε. In this case, we increased the compressive strain and repeated the procedure, described above, for an increased ε, until the zero eigenvalue for non-zero k was found. The strain ε and wavenumber k, for which the first zero eigenvalue was observed, corresponded to the buckling strain ε cr and critical wavenumber k cr , defining the buckling shape. Recall that the special case of k cr → 0 corresponds to the macroscopic instability (long-wave mode); otherwise, the composite undergoes a microscopic loss of stability, developing finite size wavy shapes upon achieving the critical level of deformation [19,20]. Figure 3 illustrates the described numerical procedure, showing a typical evolution of the dispersion curves with applied deformation for dilute composite with v l = 0.025 (a) and non-dilute composite with v l = 0.2 (b). Recall that there is an intrinsic connection between the shear wave propagation and elastic instabilities in periodic composites [24]. The continuous black curve in Figure 3a,b describes the dispersion relation in the undeformed state (ε = 0); the curves intersect with the x-axis only at k = 0, which corresponds to trivial rigid body motion. An increase in the applied compressive strain leads to a gradual change of the dispersion curves. Figure 3 illustrates these changes for the applied compressive strains ε = 0.032 (a) and ε = 0.01 (b) (dotted blue curves). Finally, for the strains ε = 0.0373 (a) and ε = 0.016 (b), we observed that zero eigenvalue appeared at the finite wavenumber (dashed red curve). Therefore, the corresponding composites, being deformed up to these critical strains lost their stability. Note, that for the dilute composite (v l = 0.025) k cr t = 2.3, which corresponds to the microscopic loss of the stability, while for the non-dilute composite (v l = 0.2) k cr ≈ 0, which corresponds to the macroscopic instability. Buckling modes of the composites with the corresponding critical parameters are shown in Figure 3 for the microscopic (a) and macroscopic (b) cases. Here and thereafter, we use the term "dilute" to refer to the composite undergoing microscopic loss of stability, as opposite to a "non-dilute" composite, which experiences macroscopic instability with formation of long-wave buckling shape. intersect with the x-axis only at = 0, which corresponds to trivial rigid body motion. An increase in the applied compressive strain leads to a gradual change of the dispersion curves. Figure 3 illustrates these changes for the applied compressive strains = 0.032 (a) and = 0.01 (b) (dotted blue curves). Finally, for the strains = 0.0373 (a) and = 0.016 (b), we observed that zero eigenvalue appeared at the finite wavenumber (dashed red curve). Therefore, the corresponding composites, being deformed up to these critical strains lost their stability. Note, that for the dilute composite ( = 0.025) = 2.3, which corresponds to the microscopic loss of the stability, while for the non-dilute composite ( = 0.2) ≈ 0, which corresponds to the macroscopic instability. Buckling modes of the composites with the corresponding critical parameters are shown in Figure 3 for the microscopic (a) and macroscopic (b) cases. Here and thereafter, we use the term "dilute" to refer to the composite undergoing microscopic loss of stability, as opposite to a "non-dilute" composite, which experiences macroscopic instability with formation of long-wave buckling shape.

Results
Before considering non-ideal composites with interphasial layers, let us firstly make some remarks on the instability in ideal hyperelastic composites without interphases. Figure 4 shows the dependence of the onset of instability on the shear modulus contrast in hyperelastic layered composites under the incompressibility assumption of all constituent materials (Poisson's ratios of layers and matrix are = = 0.5, respectively).

Results
Before considering non-ideal composites with interphasial layers, let us firstly make some remarks on the instability in ideal hyperelastic composites without interphases. Figure 4 shows the dependence of the onset of instability on the shear modulus contrast in hyperelastic layered composites under the incompressibility assumption of all constituent materials (Poisson's ratios of layers and matrix are p l = p m = 0.5, respectively).
Solid black and red lines correspond to the dilute (v l = 0.025) and non-dilute ideal composites (v l = 0.2), respectively. It is clear that in the logarithmic scale this dependence looks almost linear regardless of the instability type that the composite undergoes. Parnes and Chiskis [35] derived the estimation for the onset of instability, occurring in dilute layered composites with linear elastic constituents. Under the assumption of incompressibility of both phases in plane strain conditions, the critical strain can be estimated as: This function is shown in Figure 4 with the dotted blue curve. Remarkably, even for high critical strains ε cr > 20%, expression (5), derived for the elastic case, provided very accurate results with negligible variation from the exact value of the critical strain in the dilute layered composite developing microscopic instability mode. Another estimation for the onset of instability in the non-dilute composites, with elastic constituents originally obtained by Rosen [1], is provided in [35]. Under plane strain conditions, this estimation takes the form: As one may see from Figure 4, where expression (6) is shown by the dotted blue curve, the difference between the exact and estimated values of the critical strain is relatively small for composites with high contrast between elastic moduli; however, it increases with a decrease in the contrast. These observations allowed us to conclude that the critical buckling strain in hyperelastic neo-Hookean composites can be accurately estimated by these expressions, initially derived for the layered composites with linear elastic constituents for dilute as well as for non-dilute cases. However, we note that this good agreement may be due to the fact that the buckling develops at relatively small strains, where the linear model can approximate the nonlinear behavior.  Solid black and red lines correspond to the dilute ( = 0.025) and non-dilute ideal composites ( = 0.2), respectively. It is clear that in the logarithmic scale this dependence looks almost linear regardless of the instability type that the composite undergoes. Parnes and Chiskis [35] derived the estimation for the onset of instability, occurring in dilute layered composites with linear elastic constituents. Under the assumption of incompressibility of both phases in plane strain conditions, the critical strain can be estimated as: This function is shown in Figure 4 with the dotted blue curve. Remarkably, even for high critical strains > 20%, expression (5), derived for the elastic case, provided very accurate results with negligible variation from the exact value of the critical strain in the dilute layered composite developing microscopic instability mode. Another estimation for the onset of instability in the non-dilute composites, with elastic constituents originally obtained by Rosen [1], is provided in [35]. Under plane strain conditions, this estimation takes the form: While the continuous curves in Figure 4 represent the ideal composites, the dashed black and red curves correspond to the composites with the same µ l /µ m and v l , which contain the interphasial layers with f = 0.5 and linear variation of the shear modulus through their thickness (curve B on Figure 2). We can see that the critical strain in the composites with interphasial layers exceeds the critical strain in their "ideal" counterparts. It is worth mentioning that the dashed lines corresponding to non-ideal cases are almost parallel to the solid curves, representing ideal composites, and the difference between the critical strains (in log scale) remains virtually the same regardless of the elastic modulus contrast.
Thereby, the existence of the interphasial layers might make the layered composite more stable; however, thus far, we can state this only for linear variation of shear modulus through the thickness. At the same time, the variation of the shear modulus in the interphasial layers through their thickness Materials 2019, 12, 763 7 of 13 might be highly nonlinear, for instance, due to the complexity of the mixing and curing processes at the interphase between the different materials. Figure 5 shows the dependencies of the critical strain ε cr on the relative thickness of interphasial layers in the non-dilute composite (v l = 0.2) for different variations of the shear modulus inside the transition zone. It can be seen, that for distributions A and B, the composites with the interphasial layer were more stable than the ideal composite. Note that composites with linear distribution of shear modulus in interphasial layer (B) had the same effective macroscopic shear modulus as its ideal counterpart. The effective shear modulus is defined as an integral over the thickness of the unit cell, namely: Thus, under the uniaxial compression along the layer's directions, the incompressible hyperelastic layered composite showed the same relation between applied force and displacement as the homogeneous material with shear modulus equal to µ e f f . Since the response of the ideal composite and non-ideal composite with linear distribution B were the same for the uniaxial deformation, the observed increase in the stability of the non-ideal composite could not be explained by the change of the effective properties. Therefore, the observed increase in the stability in the non-ideal composite was directly caused by the existence of the smooth transition between soft matrix and stiff layer instead of "jump" or discontinuity of the shear modulus value on the boundary between main homogenous materials. Interestingly, for distribution C, the critical strain remainednearly the same for different values of f , and it may even be lower than the critical strain in ideal layered composite. Similar to case B, in the composite with distribution C, the smooth transition of shear modulus between soft matrix and stiff layer had a stabilizing effect on the uniaxial deformation. However, the effective shear modulus of distribution C was larger than that of the ideal layered composite, which reduced the stabilizing effect of the smooth transition zone. This is reflected by the variation of critical strain for distribution C (see Figure 5).  While the non-dilute composites lost their stability by the developing long-wave mode (wavenumber → 0), the dilute composites form wrinkling patterns with finite wavelengths = 1/ . Figure 6b shows the dependence of the critical wavenumber on the relative thickness of the interphasial layer . We can see that for the shear modulus distributions B and C, the wavenumber drastically decreased with an increase in the interphase layer thickness . Interesting behavior was observed in the composite with type A shear modulus variation through the interphasial layers. In this case, the dependence of the wavenumber on the relative thickness of the interphasial layer was not monotonic. In particular, the normalized critical wavenumber increased first, but then started decreasing with a further increase in the interphase thickness, f. Note that the effective shear modulus of the combined stiff phase (interphasial zone and stiff layers) for distribution A was lower than the one of the ideal composites; the total effective volume fraction The dependence of critical strain on relative interface thickness f was rather similar for the dilute composite, which underwent buckling by microscopic mechanisms, as shown in Figure 6a. Similar to the non-dilute case, the critical strain increased with an increase in f for shear modulus distributions A and B. At the same time, we did not observe the qualitative difference for dilute composites between the distributions A and C as in the case of non-dilute composite. Non-ideal dilute composites demonstrated more stable behavior (require higher compressive strains for onset of From Figure 7, it is clear that the interphasial layers had significant influence on the buckling behavior only if their dimensions were comparable with the dimensions of the stiff layer. In this case, non-ideal layered composite containing interphasial layer were more stable in comparison with ideal composites with the same effective shear modulus. The effect of the interphasial layer decreased with an increase in the volume fraction of the stiff layer, until the transition from microscopic to macroscopic type of instability occurred. At the same time, it should be noted that the observed convergence of the microscopic curves with an increase in volume fraction does not imply that interphasial layers in general have a marginal influence in the composites undergoing macroscopic instabilities. In fact, that Figure 7 shows the results for the composites with fixed and not and, as it was shown in Figure 5b, the interphasial layers played a significant role in the development of instabilities in the non-dilute composites as well. While the non-dilute composites lost their stability by the developing long-wave mode (wavenumber k cr → 0 ), the dilute composites form wrinkling patterns with finite wavelengths l = 1/k. Figure 6b shows the dependence of the critical wavenumber k cr on the relative thickness of the interphasial layer f . We can see that for the shear modulus distributions B and C, the wavenumber drastically decreased with an increase in the interphase layer thickness f . Interesting behavior was observed in the composite with type A shear modulus variation through the interphasial layers. In this case, the dependence of the wavenumber on the relative thickness of the interphasial layer f was not monotonic. In particular, the normalized critical wavenumber increased first, but then started decreasing with a further increase in the interphase thickness, f. Note that the effective shear modulus of the combined stiff phase (interphasial zone and stiff layers) for distribution A was lower than the one of the ideal composites; the total effective volume fraction for the combined stiffer zone was increased. However, the smooth change in the local shear modulus-in contrast to the sharp jump in the ideal composite-created a competing stabilizing effect. The contribution of these effects is reflected in the non-monotonic dependence of the normalized wave number for composites with distribution A. For the non-ideal composites with distribution A with higher values of relative thickness of interphasial layer ( f 0.6), buckling modes are characterized by lower normalized critical wavenumber as compared to ideal composites. Gao and Li [34] reported that the layered composites, in which the interphasial layers have a constant shear modulus, can buckle with a formation of patterns where the interphasial and stiff layer have different wavelengths in the buckled mode. These hierarchical patterns appear for specific combinations of geometrical parameters and materials constants [34]. However, in our study-when the transition between two main homogeneous (stiff and soft) constituents was smooth-such hierarchical buckling modes were not observed, and the interphasial and stiff layers demonstrated similar buckling shapes.
Thus far, we considered the ideal and non-ideal composites with the fixed volume fraction of the stiff layer; we examined the interplay between microscopic and macroscopic instabilities and focused on the dependence of critical strain on the volume fraction. From previous studies [3,4], it is known that layered composite undergoes the microscopic type of instabilities if the volume fraction of the stiff layers does not exceed some critical value, depending on the shear modulus contrast. Otherwise, the macroscopic loss of the stability is observed. In the case of the layered composites, containing interphasial layers, we fixed the thickness of the interphasial layer t i and found critical stretch ratio λ cr = 1 − ε cr , for which the composite lost its stability. Figure 7 shows the dependence of the critical stretch ratio on the volume fraction for the ideal composite (black continuous curve) and composites with interphasial layers (dashed and dashed dotted curves). We considered the composites with µ l /µ m = 15 and thicknesses of interphasial layers t i = 0, 0.01t, 0.025t. The interphasial layer in these composites had the distribution of the shear modulus, corresponding to the mode B ( Figure 2); therefore, all considered composites have the equal effective shear modulus µ e f f . The black dotted line represents the onset of the macroscopic instability in the ideal composite, calculated according to the explicit formula [3]: where: From Figure 7, it is clear that the interphasial layers had significant influence on the buckling behavior only if their dimensions were comparable with the dimensions of the stiff layer. In this case, non-ideal layered composite containing interphasial layer were more stable in comparison with ideal composites with the same effective shear modulus. The effect of the interphasial layer decreased with an increase in the volume fraction of the stiff layer, until the transition from microscopic to macroscopic type of instability occurred. At the same time, it should be noted that the observed convergence of the microscopic curves with an increase in volume fraction does not imply that interphasial layers in general have a marginal influence in the composites undergoing macroscopic instabilities. In fact, that Figure 7 shows the results for the composites with fixed t i and not f and, as it was shown in Figure 5b, the interphasial layers played a significant role in the development of instabilities in the non-dilute composites as well.

Conclusions
The numerical study presented in this paper showed that the existence of inhomogeneous interphases between two main constituents in hyperelastic layered composites significantly changed their buckling behavior, affecting the onsets of the instabilities as well as the developing

Conclusions
The numerical study presented in this paper showed that the existence of inhomogeneous interphases between two main constituents in hyperelastic layered composites significantly changed their buckling behavior, affecting the onsets of the instabilities as well as the developing buckling patterns. In particular, we found that for non-dilute as well as dilute cases, the "non-ideal" composites-which contain the mixing zone between two materials-were usually more stable in comparison with their "ideal" counterparts with the same effective shear modulus. However, it appears that the buckling characteristics depended not only on the thickness and effective shear modulus of the interphasial layer, but also on the distribution of the shear modulus through its thickness. Buckling responses predicted by our analysis were based on certain scenarios for distributions of the properties in the interphase material created due to mixing of the two phases during layer-by-layer curing. A careful experimental characterization is needed to provide the information on these actual distributions, and their dependence on various material fabrication and curing conditions. In the future, it is planned to perform a systematic experimental study of the formed wavy patterns in 3D-printed laminates with controllable interphase mixing zones. In Appendix A, we included a qualitative comparison of the results with limited experimental observations to illustrate the possible effect of the interphase properties on the buckling and postbuckling behavior of the periodic laminates a further careful experimental study should be performed in the future. The dependence of the buckling mode on the relative thickness of interphasial layer, shown in Figure 6a, was used to make some qualitative predictions and estimations of the role of the interphases in 3D-printed layered materials. By using Objet Connex multi-material 3D-printer, we produced the samples, which contain stiff layer embedded in soft matrix. The matrix was printed using soft hyperelastic TangoPlus resin; the layer was printed using so-called digital material DM85 with shear modulus an order of magnitude higher than in TangoPlus [28]. The digital materials in the applied inkjet technology were produced by local mixing of two base materials: TangoPlus and VeroWhite with different mass ratios. During the printing process, the print-head moves above the already printed part of the specimen injecting TangoPlus and VeroWhite inks, and the UV-light cures the deposited layer after each pass. During the inkjet printing, local mixing on the boundary between soft matrix and stiff layer occurs. As a result, the printed composites contain the "pure" homogenous soft matrix (TangoPlus), "pure" stiff homogenous layer (DM85) and the inhomogeneous interphasial layer, where the mixing occurs. We assumed that the amount of TangoPlus in the interphasial layer decreased through the thickness of the interface in the direction towards the stiff homogenous layer. Since the interphasial layer locally is the mixture of TangoPlus and VeroWhite, we treated it as a new "digital" material; and under the assumption that the amount of TangoPlus in the mixed phase decreased towards the stiff layer, we further assumed that the shear modulus of the interphasial layer monotonically increased towards the stiff layer. Another assumption was that the thickness of the interphasial layer or phase mixing zone did not change with the thickness of the stiff layer. Thus, the ratio between thicknesses of interphasial and stiff layers f changed with the since the thickness of the stiff layer. In accordance to the numerical results, shown on Figure 6a, we performed compression experiments on the 3D-printed composites with different thicknesses of stiff layer in order to measure the wavelengths of the buckling patterns. We used only one layer in our experiments; however, the critical strain and wavelength in this case matched the critical strain and wavelength of the corresponding dilute composite almost perfectly. Table A1 shows the experimental snapshots for the composites with different thicknesses of the embedded stiff layer. Regardless of the layer thickness, we observed that buckling occurred in the microscopic scenario with the formation of wrinkling patterns. Under the assumptions above about the interphase and the measurements of the wavelengths, estimates for the dependence of critical wavenumber k cr on the relative interface thickness can be summarized as shown in Figure A1. We did not know exactly the thickness of the interphasial layer, and we could only qualitative follow the change in k with an increase in the relative thickness of the interphasial layer f . amount of TangoPlus in the mixed phase decreased towards the stiff layer, we further assumed that the shear modulus of the interphasial layer monotonically increased towards the stiff layer. Another assumption was that the thickness of the interphasial layer or phase mixing zone did not change with the thickness of the stiff layer. Thus, the ratio between thicknesses of interphasial and stiff layers changed with the since the thickness of the stiff layer. In accordance to the numerical results, shown on Figure 6a, we performed compression experiments on the 3D-printed composites with different thicknesses of stiff layer in order to measure the wavelengths of the buckling patterns. We used only one layer in our experiments; however, the critical strain and wavelength in this case matched the critical strain and wavelength of the corresponding dilute composite almost perfectly. Figure A1 shows the experimental snapshots for the composites with different thicknesses of the embedded stiff layer. Regardless of the layer thickness, we observed that buckling occurred in the microscopic scenario with the formation of wrinkling patterns. Under the assumptions above about the interphase and the measurements of the wavelengths, estimates for the dependence of critical wavenumber on the relative interface thickness can be summarized as shown in Figure  A2. We did not know exactly the thickness of the interphasial layer, and we could only qualitative follow the change in with an increase in the relative thickness of the interphasial layer . assumption was that the thickness of the interphasial layer or phase mixing zone did not change with the thickness of the stiff layer. Thus, the ratio between thicknesses of interphasial and stiff layers changed with the since the thickness of the stiff layer. In accordance to the numerical results, shown on Figure 6a, we performed compression experiments on the 3D-printed composites with different thicknesses of stiff layer in order to measure the wavelengths of the buckling patterns. We used only one layer in our experiments; however, the critical strain and wavelength in this case matched the critical strain and wavelength of the corresponding dilute composite almost perfectly. Figure A1 shows the experimental snapshots for the composites with different thicknesses of the embedded stiff layer. Regardless of the layer thickness, we observed that buckling occurred in the microscopic scenario with the formation of wrinkling patterns. Under the assumptions above about the interphase and the measurements of the wavelengths, estimates for the dependence of critical wavenumber on the relative interface thickness can be summarized as shown in Figure  A2. We did not know exactly the thickness of the interphasial layer, and we could only qualitative follow the change in with an increase in the relative thickness of the interphasial layer .  As one may see from Figure A2, wavenumber decreased with an increase in the relative thickness . A qualitative comparison of the experimental observations with the numerical results  As one may see from Figure A2, wavenumber decreased with an increase in the relative thickness . A qualitative comparison of the experimental observations with the numerical results (see Figure 6b) may hint that the interphasial layers in the studied 3D-printed composites has the distribution of the shear modulus through the thickness similar to the cases B or C (Figure 2). The mixing mass ratio between TangoPlus and VeroWhite in the interphasial layer might change As one may see from Figure A2, wavenumber decreased with an increase in the relative thickness . A qualitative comparison of the experimental observations with the numerical results (see Figure 6b) may hint that the interphasial layers in the studied 3D-printed composites has the distribution of the shear modulus through the thickness similar to the cases B or C (Figure 2). The mixing mass ratio between TangoPlus and VeroWhite in the interphasial layer might change linearly through the thickness of the interphasial layer; however, the change in the shear modulus through the thickness of the mixing zone might be highly non-linear. This may explain the Figure A1. Dependence of critical wavenumber on relative thickness of interphasial layer.
As one may see from Figure A1, wavenumber k cr decreased with an increase in the relative thickness f . A qualitative comparison of the experimental observations with the numerical results (see Figure 6b) may hint that the interphasial layers in the studied 3D-printed composites has the distribution of the shear modulus through the thickness similar to the cases B or C (Figure 2). The mixing mass ratio between TangoPlus and VeroWhite in the interphasial layer might change linearly through the thickness of the interphasial layer; however, the change in the shear modulus through the thickness of the mixing zone might be highly non-linear. This may explain the observed dependence of k cr on relative thickness f . Thus, the combination of the experimental observations and numerical simulations potentially may be used for indirect characterization of interphasial layers in 3D-printed soft composites, based on their buckling behavior.