Multiscale Multiphysics Simulation in Composite Materials

1 Introduction

Nowadays, numerous materials used in engineering industrial field, as well as in biological and biomedical area, are characterized by two or more constituents, also called phases.

Carbon fiber composites, polycrystalline structures, concrete, bone or wood are only some examples of components where the constituents could be distinguished at certain length scale. These materials are called inhomogeneous. The length scale as well as the type of constituents are only some example of possible criterion to classify them. Indeed, microstructure take an important role during the analysis of the material because of the topology and the properties of each phase that determines its behavior.

It is also possible to say that defining a microscale structure means to study the behavior of each phase in order to determine the overall behavior of the global structure in terms of physical properties such as mechanical, thermal, electrical, etc.

Geometry, structure, properties and behavior of the inhomogeneous materials could be separated in two different fields depending on its length scale. On one hand, the microscale (or fine-scale) corresponds to the microstructure where each phase is distinguishable. On the other hand, the macroscale (or coarse-scale) is the largest scale (and length) of the model, in contrast to the microscale that represents the lowest scale.

The most important assumption that has to be made when studying the micromechanical models is the scale separation between micro and macro structures. That assumption allows to separate the contribution of the microscale, defined as fluctuation part, from the macroscale ones, the constant part.

As it is explained in Bohm [1] these two contributions can be distinguished as fast and slow variables, where the fast variables, microscale fluctuations, can influence the behavior of the macroscale only through their volume average. Instead, the slow variables, the macroscale contribution, are not significant at the microscale level and can be considered as locally constant.

Otero et al. [2] and Petracca et al. [3], also provides an efficient multiscale strategy that, in case of first-order homogenization (FE2), ensures the macro and micro mesh independence taking into account the conservation of the energy dissipation of through the scales.

2 Multiscale Methods

One of the most popular multiscale methods is the first-order homogenization. In this procedure, the strain obtained when analyzing the macroscopic structure is used to define the boundary conditions, applied on a Representative Volume Element RVE, to solve the Boundary Value Problem (BVP) at the microscale.

The basic principles of homogenization method were provided by Suquet [4] to obtain the constitutive equation for the homogenized properties of a heterogeneous material. The unit cell is defined as a microscopic subregion that is representative of the entire microstructure in an average sense. The RVE is employed to obtain the effective properties for the homogenized material because it is assumed that it must contain a sufficient number of heterogeneities [5] [6].

The solution of the problem at the microscale, under such conditions, acts as an equivalent constitutive law for the macroscale, and it provides material stiffness and stresses as the volume average of the microscopic ones. This equivalent constitutive law is used in all the integration points of the macroscopic model to obtain the global response of the structure. When dealing with nonlinear microstructures, it will lead to an iterative procedure in which the RVE must be solved for different boundary conditions until both scales reach equilibrium, ensuring consistency between the micro- and macroscale solutions.

The first-order homogenization technique developed assumes a scale separation between the macro and the microscale, that is, the characteristic length of the microscale should be much smaller than the length of the macroscale elements, L: l<<L [7].

The main steps of the classical FE2 technique can be resumed in the figure below. On each integration point of the discretized macroscale domain, the macroscopic strain tensor provides the input variables for the microscale domain. Then, the solution of the microscopic behavior of the RVE provides the macroscopic output and properties of the equivalent homogeneous medium.

Figure 1. Representative homogenization scheme.

In this paper, for the sake of simplicity, the procedure assumes small displacements and a quasi-static behavior in both macro and micro scales. However, the methodology proposed can be extended to other cases.

3 Discrete Multiscale Threshold Surface (DMTS) and Discrete Multiscale Constitutive Model (DMCM)

The purpose of this paper is to provide two multiscale optimization techniques, directly derived from the classical first order multiscale theory. Considering these two methods, the main hypothesis is to construct a stress-strain state collector based on multiple RVE simulations able to describe the microscale behavior. Using an Equivalent Damage Parameter ( ${\textstyle {d}_{eq}}$ ), index of the difference between the real and the elastic stress, the Discrete Multiscale Threshold Surface (DMTS) will identify with a finite number of points within the strain space the initiation of the non-linearity of the RVE. The DMTS is intended to reduce the computational cost of FE2, by identifying, during the analysis, where the RVE generation is needed and where the homogenized elastic properties are sufficient to obtain the macroscale response. As an extension of the DMTS, the Discrete Multiscale Constitutive Model (DMCM) aims to describe the complete behavior of the microstructure identifying multiple threshold surfaces related to increasing levels of ${\textstyle {d}_{eq}}$ decoupling the analysis from the FE2 methods.

Below, the authors will describe in detail all the aspects of both methods, providing examples and algorithms used for the numerical implementation.

3.1 Strain database definition

In general, due to the complexity of the RVE, we cannot know a-priori the real stress behavior of the microscale. To overcome this problem, we will study the strain space doing a preliminary analysis of the RVE for different loading cases storing the stress response of the structure. These strains, defined as strain histories, are in the multidimensional spaces ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\gamma }_{xy})}$ for 2D or ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{zz},{\gamma }_{xy},{\gamma }_{yz},{\gamma }_{xz})}$ for 3D mechanical problems.

To make accessible this information during the analysis we will store the strain histories in a Homogenized Strain Database. In this chapter, for simplicity, we will only describe the 2D case, but it is perfectly extendable to 3D.

To obtain a complete definition of the microscale behavior we will impose as strain histories the value of evenly spaced point projected on the sphere of unitary radius centered in the origin of the axes ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\gamma }_{xy})}$ . We can uniquely define these points in 3D spherical coordinates system as the combination of three parameters ( ${\textstyle \theta ,\varphi ,\lambda }$ ), as we can see in the Figure 2, where the angles ${\textstyle \theta }$ and ${\textstyle \varphi }$ represents the direction of the strain loads that we applied to the microscale and ${\textstyle \lambda }$ is the unitary strain intensity.

Then we can obtain the components of the strain vector in cartesian coordinates as:

(1)

(2)

From Eq. 1 we can observe that the strain is periodic with sin and cos functions of ${\textstyle \theta }$ and ${\textstyle \varphi }$ and these two angles can be varying between ${\textstyle [-}$ $\pi ,\pi ]$ . The subdivision of this interval determines the number of analyses to preform and the precision of the discretization.

Figure 2. Total number of analysis for 2D and 2D Reduced case and 3D Strain Space.

Introducing the parameter ${\textstyle \,m\in \,\left[1,inf\right[}$ as subdivision of the interval ${\textstyle [0,\pi ]}$ we can uniquely define the strain direction with a pair of integer parameter ${\textstyle \in [-}$ $m,m]$ , that we will call tag. With this method, each point was equally spaced from the other and the angles ${\textstyle \theta }$ and ${\textstyle \varphi }$ varying between ${\textstyle [{\frac {-\pi }{m}},{\frac {\pi }{m}}]}$ . In that way, we can observe that the total amount of analysis is ${\textstyle {\left(2\cdot m+1\right)}^{2}}$ . Considering the superposition of the same resulting strain direction the amount of strain histories will be reduced to ( ${\textstyle 2\cdot m}$ )( ${\textstyle m-}$ $1$ )+2. Despite this reduction, we will remark that using a high value of m implies an exponential increasing of the analysis, see Figure 2 above.

Once the strain loading directions are computed, we apply them to the RVE doing a classical First Order Multiscale Analysis [8]. During the analysis all the degrees of freedom are fixed and for each ${\textstyle n}$ time step we will solve the microscale problem, evaluating the homogenized stress and constitutive tensor. With this information we can determine the corresponding value of equivalent damage at each time step.

3.2 Equivalent damage definition

During a classical full multiscale analysis, the RVE is solved at each time step even if the linear elastic limit of the material is not achieved. In terms of optimization, it represents wasted time consumption. In this paper, we will introduce a key parameter, called equivalent damage, that provide this limit as the relation between the homogenized stress and the corresponding elastic one. This parameter is defined as:

(3)

where ${\textstyle {\sigma }_{real}}$ is the homogenized stress, ${\textstyle {\sigma }_{el}=}$ ${C}_{el}^{H}:\epsilon$ is the elastic stress and ${\textstyle {C}_{el}^{H}}$ is the Homogenized Constitutive Elastic Tensor.

The discretization of the RVE response with a unique parameter is needed for the interpolation of the DB information. In this way, we can obtain the stress flux surface at certain equivalent damage level and reconstruct the RVE behavior for any strain direction interpolating the stored homogenized stress data.

3.3 Discrete Multiscale Threshold Surface

Fixing the ${\textstyle {d}_{eq}}$ equal to 0.1, that the author defined as the linear elastic threshold of the microstructure, the Discrete Multiscale Threshold Surface is represented by the corresponding stress at each tag studied.

In Figure 3 and 4 we report a 3D visualization of a portion of the DMTS and the algorithm of the proposed method, where the point E represents the intersection of the direction of the stress and strain with the elastic limit surface. The points A, B, C, D correspond to the nearest tag surrounding this direction used to interpolate the results.

3.4 Discrete Multiscale Constitutive Model

Starting from the model presented above, the second proposed method is an extension of the DMTS. Indeed, DMCM provides a complete definition of RVE behavior storing in a Database the stress-strain response related to the equivalent damage parameter defined previously.

(a) (b)

Figure 3. Graphical representation of (a)Strain Space and (b) Stress Space for DMTS.

Figure 4. DMTS Algorithm.

Instead of only one value, corresponding to the elastic limit, for the DMCM multiple threshold surfaces are computed. The non-linear behavior of the RVE will be discretized with a series of equivalent damage levels ranging from 0 to 1. In this way, during the RVE analysis the strain increment is iteratively adapted to reach the correct strain increment (for each strain direction) that gives an iso-Damage surface. After that, FE2 analysis are no longer needed. Indeed, an additional interpolation method over the stored threshold surfaces could be performed to predict the behavior of the macrostructure at each integration point.

The figure below can provide how the proposed algorithm works with a graphical representation.

(a) (b)

Figure 5. Graphical representation of (a)Strain Space and (b) Stress Space Evolution for DMCM.

The points E and E’ represent the intersection of the direction of the stress and strain with the surfaces at damage i and ii. The points A,B,C,D and A’,B’,C’,D’ corresponds to the nearest tag that is surrounding this direction for the different level of damage.

The value of strain and stress for the points E and E’ were determined by interpolating the data from the DB of the nearest tag for each level of damage. In that way we can reconstruct the stress strain response of the RVE for any possible direction in the strain space.

The DMCM algorithm could be resumed in the figure 6:

Figure 6. DMCM Algorithm.

We will remark that, creating the DB with strain-controlled analysis, the reference surfaces in the strain space are concentric; on the other hand, in the stress space this is not obvious and depends on the geometry and the constitutive laws involved in the RVE.

The DB scheme provides at first level the tag. Then, for each tag, we will store the information of ${\textstyle \sigma }$ and ${\textstyle \lambda }$ (lambda) at different values of equivalent damage.

Figure 7. DMCM algorithm.

Considering the number of analyses needed for the creation of the stress-strain DB, which is proportional to the discretization parameter m, it is obvious that as long as the number of RVEs needed by the macroscale structure is higher than the number of Strain Histories, the advantages of the proposed method increases.

4 Reinforced composite stiffener section

In this section we analyze a reinforced composite beam comparing the three proposed methods: FE2, DMTS and DMCM. The microscale used is a long fiber reinforced composite, composed by F155 Epoxy Resin matrix and Carbon fiber (simulated as elastic material for simplicity). The geometry of the RVE was composed by 5 symmetric inclusions as it is shown in Figure 8 and the mesh used for the micro scale is composed by 272 small displacement linear elements with 305 nodes. The analysis was performed by using plane stress theory and considering 200mm as thickness for the stiffener.

The geometry and boundary conditions of the beam are also described in Figure 8. The macroscale structure is composed by 1075 triangular elements and 678 nodes subjected to a vertical displacement ${\textstyle {u}_{y}}$ . In order to simplify the analysis, in case of Full Multiscale and DMTS methods, only 84 elements are considered with double scale, using elastic homogenized properties in the rest of the structure.

Figure 8. Reinforced composite stiffener section and RVE geometry and mesh.

Table 1. Mechanical properties of composite RVE

F155 Epoxy Resin Mechanical Properties
Young Modulus	3.24e3	MPa
Poisson Ratio	0.32
Stress Traction Limit	80	MPa
Traction Fracture Energy	0.73	J/mm2
Stress Compression Limit	240	MPa
Compression Fracture Energy	2.19	J/mm2
Carbon Fiber Mechanical Properties
Young Modulus	235e3	MPa
Poisson Ratio	0.21

Below (Figure 9) are reported the sections of the elastic limit surface for d=0.01 and the stress evolution for d=0.1,0.5 for the composite. Figure 10 shows the force-displacement curves of DMTS and DMCM methods. As it can be seen, they overlap and achieve the same maximum value of force of the FM2 case.

Finally, in Figure 11 is shown the equivalent damage distribution on the stiffener and in the microstructure for the different models developed.

Figure 9. RVE Stress Evolution Surfaces.

Figure 10. Comparison Force-Displacement curves.

(a) (b) (c)

Figure 11. Equivalent damage distribution for (a) FE2, (b) DMTS, (c) DMCM.

In case of reinforce composite stiffener, computational effort decreases significantly using the optimization method proposed in this paper. Indeed, as we can see in Table 2, we obtain a speed-up of 4.7 and 661.2 for DMTS and DMCM method respectively. Despite of the 70 RVEs generated in the DMTS, against the 84 of the FE2, we are still having advantages using this method.

Table 2. Comparison of time and computational cost

Type	Time [s]	Memory [MB]	Active RVEs	Time Speedup
FE2	40,267.5	152.2	84/1075	-
DMTS	8,609.9	108.7	70/1075	4.7
DMCM	60.9	33.1	0/1075	661.2

Moreover, DMCM it is clearly faster than FM2 and DMTS as well. The linear interpolation of the RVE strain history can provide a non-linear behavior of the macrostructure, able to reproduce crack initiation and propagation.

5 Conclusions

First order multiscale homogenization method can fully describe both linear and non-linear behavior of complex microstructures and the impact they have on the macroscale response. However, the computational cost required to analyze large structures is not negligible. In the last years, many optimization methods have been developed to overcome this problem. In this paper, the authors propose two optimization techniques that can provide a significant speed-up compared to the classical FE2 without loss of accuracy on the final results. From the RVE FE analysis, the Discrete Multiscale Threshold Surface gives the linear elastic limit to determine where the 2-scale computations are required at the structural level, while the DMCM can conduct the analysis without solving the micro-scale model, since it can describe both linear and nonlinear regimes. The keys of both methods are the precomputed Stress-Strain DataBase and the equivalent damage parameter, essential to reconstruct the correct behavior of the structure.

DMTS and DMCM methods are validated by a reinforced composite stiffener, showing the capability and robustness of these techniques. The achieved speed-up of more than 600 times respect to FE2 and the accuracy of the results, justifies the time spent to obtain the Stress-Strain DataBase.