Co-continuous polymer systems: A numerical investigation

A finite volume based implementation of the binary Cahn–Hilliard equation was implemented using an open source library, OpenFOAM. This was used to investigate the development of droplet and co-continuous binary polymer microstructures. It was shown that the initial concentrations of each phase define the final form of the resultant microstructure, either droplet, transition or co-continuous. Furthermore, the mechanical deformation response of the representative microstructures were investigated under both uniaxial and triaxial loading conditions. The elastic response of these microstructures were then compared to a classic representative microstructure based on a face centred cubic arrangement of spheres with similar volume fractions of each phase. It was found that the numerically predicted composite Young’s modulus closely followed the upper Hashin–Shtrikman bound for both co-continuous and classical structures, while significant deviations from analytical composite theory were noted for the calculated values of Poisson’s ratio. The yield behaviour of the composite microstructures was also found to vary between the co-continuous microstructures and the representative microstructure, with a more gradual onset of plastic deformation noted for the co-continuous structures. The modelling approach presented allows for the future investigation of binary composite systems with tuneable material properties. 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/3.0/).


Introduction
Structural adhesives are increasingly used for bonding components within critical load bearing engineering structures such as aerospace and automotives. Typically these adhesives are based on epoxy polymers. Epoxies are inherently brittle due to their homogeneous microstructure and highly cross linked nature. Thus, there has been much research focused on improving the fracture toughness of epoxy polymers by incorporating a second minority phase at the nano-scale. These modifiers fall into one of two main categories: inorganic additives, e.g. silica [1,2], glass [3], alumina [4], nano-clays [5] and carbon nanotubes [6,7] or organic, usually rubber particles. Rubbery additives can be either core-shell rubber particles [8][9][10] or can form during curing via reaction induced phase separation mechanisms [11,12]. The primary energy dissipation mechanisms for rubber toughened epoxies are known to be both plastic void growth and shear band development [13]. It has also been shown that a combination of the above additives to create a hybrid material can provide synergistic toughening effects, e.g. carbon nanotubes and silica nanoparticles [14] or rubber with silica nanoparticles [15][16][17].
Kinloch et al. [12] and Brooker et al. [18] have studied the morphology and fracture toughness of thermoplastic toughened epoxy polymers. They noted a change in morphology from a spherical-particulate morphology in an epoxy rich continuous phase at low concentrations of thermoplastic phase. The morphology changed to a co-continuous structure as the percentage of thermoplastic toughener was increased. Finally at very high concentrations of thermoplastic phase, a phase-inverted structure was observed where the epoxy was dispersed as spherical particulates in a continuous thermoplastic phase. The fracture toughness was found to increase with increasing thermoplastic content, although this was not attributed to the change in morphology but rather as a function of the percentage of the thermoplastic phase in the polymer blend.
More recently, significant advances [19] in polymer synthesis has allowed researchers and manufacturers the ability to create specific copolymers. There are many possible combinations of block copolymers (BCP), and amphiphilic BCPs are of particular interest for toughening epoxy polymers. The BCP modified epoxies have been shown to have a lower viscosity than conventional rubber modified epoxies [20] which make them attractive for resin infusion processes. Several researchers have already shown that these BCPs can form complex nanostructures via self-assembly [21,22] or reaction-induced phase separation [23] to significantly increase the mechanical and fracture performance of epoxy polymers. Some of these hierarchical substructures include spherical http micelles, vesicles and worm-like micelles [24]. Of particular interest in this work is the formation of a co-continuous, or bicontinuous gyroid, structure [25]. Co-continuous structures are structures which contain two or more interpenetrating and self-supporting phases. These types of structures have been observed for thermoplastic modified epoxies but the energy dissipation mechanisms for such a complex structure were not fully understood [12,26,27]. Most of the early work on BCP modified epoxies has been primarily focused on experimentally observing the morphology and resulting mechanical properties. The study of the toughening mechanisms and mechanics of deformation of the microstructure is required in order to fully understand the structure/property relations.
Furthermore, such complex designed binary ordered structures can provide outstanding combinations of material properties including stiffness, toughness and strength that would otherwise not be possible in a homogeneous material and have been the subject of intensive research efforts over the last number of years. Wang et al. [28] studied the energy dissipation in binary composite systems with the phase interface closely approximating a triply periodic minimal surface. Similarly, work by Lee et al. [29] has shown that the energy dissipated per unit volume in a voided epoxy nanostructure can be significantly enhanced for a system of lower relative density by close control of the nano-frame geometry. Torquato [30] has studied the developed of optimised co-continuous structures with targeted properties such as minimal thermal expansion. It is therefore clear that understanding the underlying physics of phase separation and the subsequent influence of the morphology on the mechanical properties of the composite blend, is key to guiding the future design of improved materials for structural applications.
In this work, we use the well-known Cahn-Hilliard equation to describe the temporal evolution of a binary phase field of two mixed polymers and subsequently investigate the mechanical response of the evolved microstructures under different loading conditions.

Cahn-Hilliard model
The Cahn-Hilliard equation was originally proposed in 1958 [31] to model phase separation phenomena in binary alloys. Since then, it has been applied to a wide number of other diverse uses including the evolution of arbitrary microstructures [32,33], tumour growth [34] and image in-painting [35]. The Cahn-Hilliard equation can be written as: @c @t where c is the concentration of one phase in the mixture, t is time, M is defined as the mobility, and l is the local chemical potential, more commonly defined as: In Eq. (2), FðcÞ represents the free energy density of a homogeneous material of concentration c, and c is a constant related to the thickness of the interfacial regions between phases. The free energy, FðcÞ, is based on a Ginzburg-Landau functional of the form: The first term on the right hand side of Eq. (3) represents the homogeneous free energy while the second term penalises large gradients. The free energy is generally deduced via the Flory-Huggins model [36]: where h is a variable representing the Flory-Huggins interaction parameter, and R and T are the universal gas constant and temperature respectively. From Landau theory, an alternative quartic approximation can be written as [36]: where a; b and c can take suitable values. In order to derive an easily implementable numerical solver for Eq. (1), we first define an order parameter / which represents the concentration difference between phase A and B at any point in the domain of interest. Therefore, we have / ¼ c A À c B with c A þ c B ¼ 1. It is easy to observe that / ¼ 2c À 1. Eq. (1) can then be rewritten as: It can be seen in Fig. 1 that the new definition of Fð/Þ ¼ 1 4 ð/ 2 À 1Þ 2 is a double well potential with local minima at / ¼ fÀ1; 1g corresponding to the pure phases. The initial concentration of each phase in the mixture determines the subsequent pattern evolution. Specifically, for near equal concentrations of phase, a co-continuous pattern emerges. However, in the case where one phase is present in the mixture in a much greater concentration than the other, a droplet pattern emerges. In this case, one phase can be considered as the matrix and the other as dispersed spherical inclusions.

Two-dimensional cases
Using the model described, we can examine the evolution of the concentrations of two phases with a random perturbation in a two-dimensional system. First, we prescribe the concentration difference between each phase as an initial condition (the concentration difference as defined by Eq. (6) can take any value between À1 and 1). To initialise the random perturbation, we chose a Gaussian normal distribution about the mean concentration with a standard deviation of 0.1, i.e. 5% of the total range.
A unit square domain with equal sides of 100 lm with a cell size of side length 1 lm was used. A time step of 0.1 ms was used in all cases while the mobility, M, was set to 0.01 m 2 /s, with c (which controls the width of the interface region) equal to 1 Â 10 À4 m 2 /s. A number of different simulations were performed at various volume fractions to examine the role of volume fraction on the resultant morphology. The boundary conditions were set to be periodic. Fig. 2 presents the temporal evolution of the phase morphology with an initial mean concentration of 30%. It can be seen that the mixture decomposes into discrete phases, comprising of a matrix (blue 1 phase) and dispersed droplets (red phase). The green phase demonstrates that in the Cahn-Hilliard implementation the phase boundaries are not discrete, but rather are diffuse over a given length. For very small particles the gradient term in Eq. (3) becomes dominant and the homogeneous mixture becomes the favourable energy state. This results in the preferential growth of larger particles over time as can be observed at t ¼ 10 and t ¼ 40. Fig. 3 presents the temporal evolution of a binary mixture with an initial concentration of 50%. A co-continuous morphology is observed, rather than the discrete decomposition observed in Fig. 2. Again, coarsening of the phases over time is observed. The characteristic length of each phase therefore changes over time.
The phase structures modelled can be observed experimentally. Fig. 4 shows the phase images of a block copolymer modified epoxy at 2 concentrations obtained using atomic force microscopy (AFM). The epoxy polymer is a polyetheramine (Jeffamine D-230 from Huntsman, UK) cured diglycidyl ether of bisphenol A (DGEBA) resin (Araldite LY556 from Huntsman, UK), and the modifier is a symmetric triblock copolymer (Nanostrength M52N from Arkema, France). The block copolymer was received as a dry powder and dissolved into the epoxy precursor by mechanical stirring at 120°C. The curing agent was then added and the mixture degassed in a vacuum oven. The epoxy is cured at 75°C for 3 h, followed by a ramp of 1°C/min to the post-cure temperature of 110°C for 12 h and cooled to room temperature in the oven.
A direct comparison between experimental and numerical results is not possible as the numerical procedure presented is a model structure for simplicity. The model assumes that the viscosity and hence mobility of each phase is both constant and equal to each other. Experimentally, this is not the case, allowing for the emergence of co-continuous structures at a much lower concentration than in the numerical cases presented here. Comparing Figs. 4(a) and 2(a), there are some clear similarities for the low initial concentration formulations. Spherical particles with a mean radius of approximately 117 ± 6 nm were measured from the AFM micrographs. At 10 wt% (see Fig. 4(b)), the co-continuous structure is comparable to that in Fig. 3(a). The high conversion rate of the epoxy results in a relatively quick curing time, thus the phases do not have the time to grow in size (see Figs. 2(c) and 3(c)). Moreover, the block copolymer in Fig. 4 is present in much lower concentrations than that investigated numerically.
It is important to be able to measure the characteristic length of each phase as the polymer morphology can significantly influence the resultant bulk mechanical behaviour. To this end, two point correlation functions were used to characterise the phase morphology as a function of time. These functions have been used previously by a number of researchers to quantify microstructures [50][51][52]. First, the diffuse phase boundaries are converted to a binary representation via an image thresholding operation. A set of functions, P ij , is then defined which denote the probability that a randomly placed vector of a given length, r, begins in phase i and ends in phase j (i; j = 0 or 1). The lighter phase in Fig. 5 is defined as phase 0, while the darker phase is designated as phase 1. We define f as the volume fraction of phase 1 and note that at the limits r ! 0 and r ! 1, the correlation functions are completely defined by f [53,54], i.e.: lim r!0 lim r!1 Additionally, it can be shown that the minimum value of P ii , or maximum value of P ij gives a reasonable estimate of the characteristic lengths of the microstructure. In the case of a dispersed structure this corresponds to the diameter of the dispersed phase or the distance between dispersed particles, while for a co-continuous structure, it represents the average thickness of each phase. Fig. 5 plots the two point correlations for a co-continuous structure for a range of values of r (r=r max is the ratio of the random vector length to the side length of the overall domain). The coarsening process was first formalised by Lifshitz and Slyosov [55] and Wagner [56] which predicts that the average particle size is linearly proportional to the cube root of time. Fig. 6 plots t 1=3 versus normalised structural length scale. It is clearly seen that the coarsening process in our numerical simulations obeys the Lifshitz-Slyosov-Wagner model. It can also be noted that the growth rate of the characteristic length is greater for c ¼ 0:5 than for c ¼ 0:3. The deviation from the straight line growth relationship can be attributed to the fact that the correlation function does not take into account the periodic nature of the system and that only 1000 random vectors were sampled at each value of r. In order to check that the numerical model has been implemented correctly and is valid, it is necessary to ensure that it is mass conservative, i.e. there must be no spontaneous creation or destruction of one phase over the other, and that the free energy in the system should be strictly non-increasing. Fig. 7 plots the normalised discrete energy calculated via Eq. (3) and the mass fraction of the minority phase versus time. This demonstrates that the developed numerical solver is both mass conservative and that the total energy in the system is always decreasing.

Three-dimensional cases
In order to examine the energy dissipation mechanisms in these polymer structures, three dimensional     c ¼ 0:5. For each system, both an early stage decomposition structure and a late stage decomposition structure were chosen for analysis as shown in Table 1. In each case the minority phase is shown in red with the majority phase in blue. It can be seen from Figs. 6 and 7 that in the early stage structures, the characteristic lengths of the microstructures are much smaller and this could influence the mechanical response of the composite. The choice of early or late stage structures is also influenced by the fact that in experimental studies the non-equilibrium co-continuous microstructures are trapped by gelation of the matrix phase. The resultant microstructures were created by taking the concentration field and thresholding c such that c ¼ v f , where v f is the volume fraction of the minority phase. Furthermore, representative volume elements comprising 1=8 of a face centred cubic structure, therefore consisting of spheres, corresponding to each value of v f were also analysed to compare against the mechanical response of the co-continuous structures as shown in Table 1. This unit cell arrangement is popular amongst researchers as it is the simplest structure which allows for investigation of phase interactions and each particle is equidistant from each other [57]. In the cases where c ¼ 0:4 and c ¼ 0:5 a co-continuous structure was observed to have formed, whereas for c ¼ 0:26 an intermediate structure that is between a co-continuous and discrete spherical particles, was observed. Fig. 8 details the difference between a co-continuous structure and a so-called transition structure with the second phase removed for clarity.

Stress analysis
Finite volume based stress analysis was carried out on each of the microstructures shown in Table 1 using foam-extend-3.0 [38]. Phase 1 was treated as a linear elastic-plastic solid, while phase 2 was treated as linear elastic, with the properties as given in Table 2. It should be noted that the elastic and cohesive properties chosen are not indicative of any particular material but can be thought of as representing a generic epoxy-rubber system. In general the minority phase (red phase) represents the rubber phase except when phase inverted structures are discussed.
The volume averaged stress, r c , and strain, e c , were found by averaging the local stress and strain in each cell using Eqs. (13) and (14) [58]. This is known as homogenisation.
where V X is the total volume of integration and r and e represent the local stress and strain at any point within the cell.

Elastic behaviour
Each microstructure was loaded using a fixed axial strain rate of _ e u ¼ _ e h ¼ 1 Â 10 À3 s À1 in both uniaxial extension and triaxial extension configurations. _ e u and _ e h represent the uniaxial strain rate and hydrostatic strain rate respectively. This yields the composite bulk and constrained moduli respectively via appropriate application of Eqs. (13) and (14) and Hooke's Law, allowing complete characterisation of each composite in the linear elastic regime. Hashin and Shtrikman [59] developed theoretical upper and lower bounds for the bulk modulus, k, and the shear modulus, l, of quasi-isotropic composites with an arbitrary phase geometry: where l and u represent the lower and upper bounds respectively, m is the volume fraction and the subscripts 1 and 2 represent the two phases in the composite. From this it is possible to calculate upper and lower bounds for the Young's modulus, E, and Poisson's ratio, m, using the usual relations (for i equal to u and l): In the case of a voided epoxy structure, it is appropriate to use the scaling laws for cellular structures. Gibson where q Ã is the relative density of the foam, q s is the density of the foam material and E s is the Young's modulus of the foam material. C is an empirical constant which is approximately equal to 1 for an equi-axed cellular structure. Furthermore, Gibson and Ashby suggested that the Poisson's ratio for an open cell foam can be adequately approximated as 0.33, although they concede that there can be a large variability in this property between structures, which they attribute to differences in cell geometry. A simple adaptation of this formula to take into account the presence of the second phase in the binary system, i.e. E Ã ¼ E r at q Ã ¼ 0 and E Ã ¼ E s at q Ã ¼ 1, can be written as: where E r represents the Young's modulus of the compliant phase and the other variables have the same meaning as in Eq. (19). Fig. 9 shows the effect of volume fraction on the Young's modulus and Poisson's ratio of the different microstructures discussed in Table 1. The higher volume fractions of 0.6 and 0.74 were obtained by simply inverting the material phases of the microstructures with v f ¼ 0:4 and 0.26, except in the case of the representative face centred cubic structures where the soft rubbery phase was confined to the discontinuous phase at the corners of each microstructure. This is to maintain the stiff central columns in the directions of loading [28]. It should be noted that it is not possible to create a representative face centred cubic structure for v f ¼ 0:8 as the maximum packing fraction for spheres of the same size is approximately 0.74048 [61]. The Hashin-Shtrikman bounds and the adapted empirical relationship (Eq. (20)) are also plotted in Fig. 9. It can be seen that the Young's modulus, E c , of each of the structures investigated lies close to the upper Hashin-Shtrikman bound. This is not surprising given that the elastic modulus response would be dominated by the much stiffer component in the composite. There is no significant difference between the predicted Young's modulus of the co-continuous structures and the representative sphere based structures. However there is a significant difference in the predicted Poisson's ratio, m c , and this is due to the mutual constraints of each phase within the complex co-continuous geometry.   The effective elastic properties of a voided structure were also computed. These structures were achieved by removing the rubbery phase prior to simulation. The numerical predictions for E c and m c are given in Fig. 10. The predictions for E c lie close to the upper Hashin-Shtrikman bounds, but in this case the empirical Eq. (19) suggested by Gibson and Ashby gives a better fit to the data. As in the case of the epoxy rubber microstructures, there is little difference in E c between the co-continuous structures and the representative face centred cubic structure.
However, it can be seen that the agreement of calculated Poisson's ratio between the co-continuous structures and representative structures diverge significantly as the volume fraction of the voids increase. In the case of the representative microstructures, the calculated Poisson's ratio agree with the suggestion by Gibson and Ashby that for any open cell structure, m % 0:33. By extrapolating the Poisson's ratio data for the co-continuous structures investigated in the current work, it can be expected that an auxetic structure exists at high void volume fractions. Indeed, there is significant research interest in auxetic materials (À1 < m < 0) with microstructures similar to those investigated in the current work [62].
In the case of the epoxy-rubber system, the high bulk modulus of the rubber phase provides an additional constraint to the epoxy preventing local deformation into the regions occupied by the rubber phase. In the case of the voided structure, where a network of voids replaces the rubber phase, no such constraints are present and the epoxy phase is free to locally deform into the voided network under a macroscopic load. This is the primary reason for the different mechanical responses between the epoxy-rubber microstructures and the voided microstructures.

Post-elastic behaviour
In order to simulate the yield behaviour of the composite microstructures, it is important to be able to understand the state of stress under a range of different types of loading. The constraint factor, H ¼ r H =r eq ¼ 1=3, for a bulk material under uniaxial loading conditions where r H is the hydrostatic stress and r eq is the von Mises stress. For other forms of loading, H is typically greater than 1/3. Where a crack exists in a material, the crack tip region is normally expected to be under a highly triaxial state of stress. For an adhesively bonded tapered double cantilever beam geometry, Hadavinia et al. reported H values of between 1.5 and 2.3 for a range of bond-gap thicknesses [63]. Cooper et al. attributed the bond-gap thickness effect on fracture toughness for a rubber toughened epoxy adhesive to the variation in constraint ahead of the crack tip rather than the variation in 'far-field' plasticity [64].
Therefore, in the current work, the yield and post-yield behaviour of each model microstructure is studied under triaxial extension. Fig. 11 plots the homogenised applied strain versus von Mises stress for triaxial loading for both the sphere based structures and the late stage co-continuous structures. The 0.2% offset proof stress is also shown for each case. The difference in behaviour between the face centred cubic structures and the physically based co-continuous is emphasised in Fig. 12.
It can be seen that as the volume fraction of the rubbery phase increases the 0.2% offset stress decreases, with a lower value for the co-continuous structure being recorded in all cases. A similar decrease with increasing volume fraction is noted for both the sphere based FCC structures and the co-continuous structures. No significant effect is noted in the transition from droplet type co-continuous structures to fully co-continuous structures. It is not surprising that the co-continuous structure has a lower proof stress than the corresponding FCC structure. During the evolution of the microstructures, the phase interface rarely has a constant curvature as in the case of the FCC structures (curvature, K ¼ 1=R 2 where R is the radius of the sphere). In some areas along the phase interface the local curvature is much greater, leading to an increased stress concentration factor at that point. This leads to the material yielding locally at a lower global stress and consequently a lower proof stress is measured. This implies that the yield stress is only dependent on the local stress concentration factors at the phase interfaces within the microstructure.
For the sphere based FCC structures, the strain at which the 0.2% offset is reached is steadily predicted to decrease marginally from 1.27% for v f = 0.2 to 1.16% for v f = 0.5. Whereas, in the co-continuous structures, the trend for the 0.2%, with an increasing value of v f , initially follows that for the offset stress before increasing again at high volume fractions. This is directly attributable to the change in microstructure from a droplet type microstructure to a fully-cocontinuous structure at high volume fractions. Therefore, the effect of a co-continuous structure appears to be to increase the proof strain. On the other hand, the proof stress is independent of morphology since it is only dependent on the local stress concentration factors at the phase interfaces. The variation in 0.2% yield stress agrees with the trends noted by Wang et al. [28] in their study of co-continuous structures defined by minimal periodic surfaces.
In the co-continuous structures localised plasticity occurs early, and at multiple locations in the epoxy phase, before gradually spreading throughout the entire structure. This is illustrated by the more gradual transition between the linear elastic regime and the fully plastic regime compared to the representative structures. The location of this early stage plasticity coincides with junctions in the geometry where multiple ligaments of the epoxy phase connect. In the case of the face centred cubic microstructures, plasticity does not occur until later in the loading history and then spreads rapidly between each rubber particle. This is clearly demonstrated by the fringe plots in Fig. 13 which detail the maximum principal strains in both the face centred cubic microstructure and co-continuous microstructure at 1% and 2% axial applied strain.
Finally, it is worth noting that once plasticity has diffused throughout the entire structure, i.e. when e ii ¼ 5%, the equivalent stress, r eq , in the co-continuous structure is higher than that in the corresponding representative structure.

Conclusions
A method of generating physically realistic co-continuous microstructures based on the Cahn-Hilliard equation has been developed and implemented within an open-source finite volume package. The method was shown to accurately capture the physics of phase separation with a minimum of computational effort. A number of representative microstructures generated using the Cahn-Hilliard equation have been investigated under different mechanical loading conditions and compared to conventional representative structures, namely 1/8 face centred cubic (spherical) structures. The results reveal that the mechanical properties of the co-continuous structures can differ significantly from the classic structures, and must be taken into account when analysing such complex nanostructures. The results presented in this paper provide a numerical roadmap for developing co-continuous polymer composites with tuneable mechanical behaviour and energy absorption properties. For example, Suo and Lu [65] demonstrate the possibility to create self assembling periodic structures by introducing discontinuities at a coarse level in the concentration field in a binary system such as the one studied here. Clearly, such structures have the potential to offer a unique combination of mechanical, thermal and electrical properties.