Stochastic Finite Element Modelling of Char Forming Filler Addition and Alignment – Effects on Heat Conduction into Polymer Condensed Phase

Abstract

Micro- and nano-filler particles have been considered as char-forming flame retardants for polymers. It has been shown that suitable particles may operate in the condensed phase to prevent or delay the escape of fuel into the gas phase. Good flame retardancy performance may be achieved in composites with comparatively low filler loadings. However, many candidate filler materials, such as rod-like and plate-like carbon allotrope fillers with high aspect ratio, will effectively enhance the composite’s thermal conductivity, and hence, may greatly increase heat input into the condensed phase. Moreover, anisotropy in terms of thermal conductivity must be considered when rod-like and plate-like particles are aligned, for example as a result of manufacturing processes. The presented study investigates these effects, i.e., thermal conductivity enhancement due to filler addition and alignment, using a modeling framework based on Monte Carlo simulation that was developed for predicting effective composite properties considering filler-matrix and particle-to-particle interfacial effects. A stochastic finite element analysis was performed to model rod-shaped carbon particles embedded in a polymer matrix. The chosen analysis is demonstrated to be an effective means for elucidating the effect of filler addition and alignment on the heat conduction into polymer materials containing fillers as char-forming flame retardants.

1. Introduction

The advent of engineered nanoparticles with high aspect ratio, such as graphene and carbon nanotubes (CNT), and their availability in quantities relevant for industrial production, has greatly expanded opportunities to modify polymers to meet demanding requirements in a broad range of applications. Such nano-additives have been shown to improve polymer mechanical (e.g., stiffness, strength and fracture properties) as well as physical characteristics (e.g., electrical and thermal conductivity), see e.g., [1, 2, 3, 4, 5, 6]. The same holds true in the context of flame retardancy. Among the three commonly considered flame retardant approaches—i.e., gas phase flame retardants, endothermic flame retardants, and char-forming flame retardants—nanofillers are typically active via the latter mechanism. Nanofillers operate in the polymer condensed phase where they may provide thermal insulation and a mass transport barrier that mitigates the release of fuel into the gas phase. Nanocomposites with suitable filler morphology and loading were observed to form a coherent filler network layer covering sample surfaces, which significantly reduced peak heat release and radiant heat flux [7]. In addition to the char-forming mechanism, nanofillers were found to reduce the melt flow of polymers. High aspect ratio nano-additives were reported to form jammed network structures causing melt to behave rheologically like a gel, thus inhibiting dripping of flammable material [8].

While the potential of nano-fillers to enhance flame retardancy through increased barrier properties impeding heat flux and fuel release, and altered rheological properties inhibiting flammable drips, has widely been acknowledged in the technical literature, the influence of filler addition on increased thermal conductivity and thus heat transfer into the polymer still requires further study [1]. Carbon-based fillers possess thermal conductivities that exceed those of polymers by several orders of magnitude. For example, thermal conductivity ranging from 2000 to 5000 W m⁻¹ K⁻¹ has been reported for CNT and graphene while values for typical polymers are between 0.1 and 0.3 W m⁻¹ K⁻¹ [9, 10]. Assessing and understanding thermal conduction in nanocomposites with high aspect ratio fillers is particular complex due to the inherent propensity of filler contact, alignment and agglomeration.

Besides randomly oriented and dispersed particles, polymer nanocomposites with purposely aligned particulate fillers have been created, which resulted in improved performance in a variety of applications. Nanocomposites with aligned particles have been used for the design of sensors [11, 12, 13] and high-strength modified polymers that require particle alignment in order to achieve specific anisotropic properties [14, 15, 16]. Carbon nanotubes, as a ‘one-dimensional’ high aspect ratio carbon allotrope, are especially suited to create nanocomposites with anisotropic properties, e.g., in terms of heat transfer properties [17, 18, 19, 20, 21, 22].

Determining the mechanical and physical properties using experimental methods is typically a time-consuming and costly approach. Analytical methods, on the other hand, are highly efficient for predicting effective material properties of particulate composites [23]. However, analytical methods lack accuracy when predicting properties, especially for higher filler volume fraction modified polymers. Considering these limitations, and in light of a rapidly growing number of applications involving particulate composites, experimental and analytical approaches are not sufficient to address the demands imposed by a vast field of available filler materials and fabrication parameters. Hence, alternative methods for assessing and designing filler modified composites need to be developed [24, 25, 26].

Stochastic analysis is one of the most reliable and recognized methods for analyzing complex problems involving many input and output parameters in the field of reliability analysis. This method can predict accurate outcomes using statistical principles. Stochastic analysis can be used in a variety of applications, e.g., financial forecasting and modeling, where numerous input and output variables need to be considered. Recently, a numerical modeling framework was developed based stochastic analysis to simulate the effective material properties of filler modified composites [27]. Specifically, a stochastic finite element analysis (SFEA) approach was employed that enabled the prediction of the effective thermal conductivity of randomly distributed and disperses spherical particles embedded in a polymer matrix. In the present contribution, aforementioned SFEA approach was adopted to predict the effective thermal conductivity of a polymer matrix containing randomly oriented or aligned rod-shaped filler particles. The particle geometry was adapted to mimic CNT. The study described herein investigates the effect of filler addition and alignment on heat transfer into polymer composites in the context of fire-retardant materials.

2. Stochastic finite element analysis framework

The nature of stochastic analysis, and thus the presented modeling approach, requires performing numerous iterations in order to calculate the effective thermal conductivity of a polymer matrix containing a rod-shaped filler. The SFEA algorithm described in [27] was adopted and employed for the present study. This algorithm provides a framework for connecting a customized stochastic analysis with a parametric finite element analysis (FEA). In this manner, the process of applying uncertainty to input variables, and creating and solving FEA models is automated. The modeling framework, which uses several scripting languages, is briefly summarized in the present section. The interested reader is referred to [27] for additional details on the modeling approach. Figure 1 depicts a flowchart outlining the main domain, i.e., the elemental structure and connections, of the framework’s various modules.

Figure 1.

Stochastic finite element analysis framework [27].

The main domain was developed in Visual Basic for Applications scripting language (VBA; Microsoft, Redmond, Washington, USA). Input parameters are provided via the ‘Front End’ module. The set of required input parameters comprises (i) the modeling domain that is defined by the size of the considered cubic representative volume element (RVE); (ii) a set of filler volume fractions that are to be analyzed; (iii) the material properties for the filler and matrix; (iv) the filler particle size distribution; (v) information on boundary effects, i.e., particle-to-particle and particle-to-matrix interfacial thermal resistance (ITR) as well as a threshold gap size that defines direct contact between particles and particles to the RVE boundary; (vi) details for the FEA mesh generation; and (vii) details regarding the model output acceptance criteria required for statistical analyses, i.e., standard deviation and variance. The input parameters are transferred to a database with an appropriate management system (‘DBMS’), which holds and communicates input and computed data between the various modules.

The Monte Carlo simulation (MCS) module, also developed using VBA scripting language, performs two subprocesses, i.e., the random number generator (RNG) and the FEA modeling. Using the algorithm depicted in Figure 2, the MCS module retrieves needed input parameters from the database, performs iteratively the SFEA, computes statistical data (standard deviation and variance) after each iteration, and finally stores results back into the database. The MCS module repeats the modeling subprocess until the acceptance criteria defined in the database are satisfied. Once results converge according to the criteria specified, the MCS module determines the effective thermal conductivity (by calculating an average value). The MCS module repeats the above processes until all specified filler volume fractions are analyzed.

Figure 2.

Monte Carlo simulation (MCS) module algorithm [27].

The RNG module was developed in the numerical computing environment MATLAB (MathWorks, Natick, Massachusetts, USA), which has pseudorandom number generating capabilities. This module retrieves the input data defining the RVE size and the particle size distribution from the database. The RNG module sequentially creates sets of random numbers for anchor points in Cartesian coordinates as well as vectors required for generating the position and orientation of rod-shaped filler particle geometries, respectively. The RNG module also performs a collision detection using a geometrical model to avoid particles interfering with each other as well as with the RVE surfaces. When detecting interference, the RNG rejects the most recently generated particle. In the geometrical model, the rod-shaped particle geometry is represented by a series of spheres (see Figure 3). The distance between each sphere associated with the most recent rod-shaped particle and all preexisting sphere geometries in the RVE, and the RVE boundaries, is evaluated to discern a particle collision. While representing rod-shaped particles using a series of spheres is only an approximation, it provides an expedient means for performing the collision detection algorithm. For the analyses presented herein, as series of 200 spheres was used to represent rod-shaped particles for interference detection.

Figure 3.

Schematic of series of spheres representing rod-shaped particles.

The process performed by the RNG module can be controlled to yield both randomly oriented and aligned rod-shaped particles within the RVE. In the case of aligned particles, constraints are imposed on the vector indicating particle orientation. As shown in Figure 4, after generating the set of random numbers for each particle, these data are stored in the database in a tabulated format for later use by the FEA module. The RNG module iteration is terminated when the required filler volume fraction is reached.

Figure 4.

Random number generator algorithm [27].

It should be mentioned that the RNG module also has the ability of creating particles that conform to a given size distribution (whilst, this feature was not utilized in the present study). The interested reader is referred to [27] for details on the algorithm that produces particles obeying a certain size distribution, and the effect that different size-ordered particle addition has on computational performance.

The FEA module was developed as a fully customizable parametric FEA platform in ANSYS Workbench (Version 19, ANSYS Inc., Canonsburg, PA, USA) using IronPython scripting language, which enabled applying uncertainties to input parameters required for performing the FEA simulation. This platform consists of a model generation environment, i.e., ANSYS DesignModeler, and a model solution environment, i.e., ANSYS Mechanical, which enable creating the parametric geometry and the finite element model, respectively. JAVA scripting language was used to automate the process of reading input data from the database (i.e., RVE dimensions and particle anchor points and orientation vectors) and creating particle geometries in the model generation environment. The three-dimensional geometry thus created is transferred to the model solution environment for further analysis. Similar to the model generation environment, the model solution environment also uses JAVA scripting language to automate the FEA process. The model solution environment retrieves further inputs parameters from the database, including material properties, information on boundary effects and conditions, and mesh generation parameters, and then constructs the finite element model for each model iteration. After performing the analysis, the FEA results are saved to the database in tabulated format for further statistical analysis. Figure 5 illustrated the algorithm for the FEA module.

Figure 5.

FEA module algorithm [27].

3. Steady-state numerical modeling

The developed SFEA framework was employed to estimate the effective thermal conductivity under steady-state conditions of filler modified composites with randomly and aligned rod-shaped particles embedded in a polymer matrix. The rod-shaped particles mimic CNT embedded in epoxy polymer in order to elucidate the effect of filler addition and alignment in the context of fire-retardant materials. Table 1 shows the CNT longitudinal and lateral thermal properties and volumetric mass density, which were adopted from [9, 28, 29]. The mean particle diameter was set to 2.85 nm with a constant particle aspect ratio of 56, i.e., only a single particle size was utilized in this study to limit the parameter space affecting the results. The RVE size was set 200 nm.

Material	Epoxy	CNT
Density [kg m^−3]	1250	1600
Thermal conductivity, longitudinal [W m^−1 K^−1]	0.25	3500
Thermal conductivity, lateral [W m^−1 K^−1]	0.25	1.5

Table 1.

Volumetric mass density and thermal conductivity of polymer and CNT [9, 28, 29].

Since CNT have anisotropic thermal properties it is not possible to define their thermal conductivity using global coordinates. Hence, an algorithm was developed in JAVA scripting language that provides dedicated local Cartesian coordinates for each rod-shaped particle. As depicted in Figure 6, the local coordinates (x,y,z) have their origin at one end of a particle with the x-direction aligning with the particle’s longitudinal axis. For the case of aligned particles the components describing the vector for each particle’s major axis (x) were constrained to remain within upper and lower bounds. For the present analyses, these constraints correspond to a maximum possible angle of approximately 8.5° between a particle’s major axis (x) and the global (RVE) X-direction.

Figure 6.

Illustration of local coordinate systems for rod-shaped particles.

Three-dimensional ten-node quadratic tetrahedral thermal solid elements, i.e., SOLID87, were used to generate the finite element mesh for both the particles and matrix. This element type, which provides one degree of freedom (temperature), is recommended for meshing irregular geometries. The latter characteristic is desirable in the present context, given that the rod-shaped particles constitute geometries that typically are difficult to mesh. As an example, Figure 7 depicts the meshing generated for the matrix (left-hand side) and randomly distributed and aligned rod-shaped particles occupying the RVE devoid matrix (right-hand side).

Figure 7.

Finite element mesh of matrix (left) and aligned rod-shaped particles occupying the RVE devoid matrix (right).

As demonstrated in [27], it is essential to model the ITR between particles and the matrix as well as between particles that are in contact with each other in order to achieve a model that realistically captures effective thermal conductivities for different filler loadings. In this study the particle-to-matrix thermal contact conductance (TCC) was adopted from literature [30, 31, 32] as 10⁸ W m⁻² K⁻¹. Also, the direct particle-to-particle heat transfer threshold was set to approximately 1 nm. Note that implementing this threshold is necessary since the employed particle collision algorithm prevents true direct particle-to-particle contact. ITR and particle-to-particle thermal contact was implemented using three-dimensional 6-node quadratic surface-to-surface elements, i.e., CONTA174 and TARGE170. For details on the chosen approach to model contact phenomena the readers is referred to [27].

Thermal boundary conditions were applied to the RVE to perform the stead-state thermal analysis and calculate effective thermal conductivities. A temperature 22 and 32°C were defined on opposite sides of the RVE, respectively, with the remaining surfaces considered adiabatic. Note that for the case of aligned filler particles, the alignment direction is referred to as the longitudinal direction, as opposed to the two lateral directions that are defined in a Cartesian coordinate system. Three sets of boundary conditions were also applied for randomly distributed filler particles for determining the effective thermal conductivities along the global Cartesian coordinate directions (X,Y,Z), thus enabling the assessment of isotropy. The applied boundary conditions create a temperature gradient and thus a heat flux between the warm and cold RVE surfaces, as illustrated in Figure 8. The thermal conductivity of nodes K_i located in the warm surface was determined using Eq. (1).

$$K_i=\frac{Q_i\times l}{T_2-T_1}$$ E1

where Q_i is the calculated numerical total heat flux at the ith node located on the warm side, $l$ is the RVE length, and T₁ and T₂ correspond the temperature on the warm and cold surface, respectively. Consequently, Eq. (2) yields the effective thermal conductivity, K_eff.

$$K_{\mathrm{eff}}=\frac{{\sum }_{i=1}^n{K}_i}{n}$$ E2

where $n$ is total number nodes on the warm surface of the RVE.

Figure 8.

Illustration of longitudinal (left) and lateral (right) heat flux between opposing warm and cold RVE surfaces.

4. Results and discussion

The described modeling framework was employed to calculate effective thermal conductivities of composites with randomly distributed particles that were either aligned or had random orientations. The studied filler volume fractions were 2.0, 4.0, 7.5 and 10%. Note that successful mesh generation becomes challenging for high filler aspect ratios, and hence, the particle aspect ratio was limited to 56 in the current study. While this value is comparatively low for CNT it does represent actual (multiwall) CNT structures as indicated in [33]. Moreover, modeling filler volume fractions exceeding 10% was found to demand excessive computational effort, and hence, analyses were limited to 10% filler volume fractions and below. Note that effective enhancement of flame retardant properties was ascertained in CNT-polymer composites that were significantly below the set 10% limit, see e.g., [7].

A convergence study was performed for a composite with randomly distributed and aligned particles at a filler volume fraction of 4.0%. To assess the sensitivity of the computed effective thermal conductivity to mesh refinement, numerical analyses were performed at different levels of mesh refinement. The results are depicted in Figure 9. It was observed that changing mesh density from ∼266,000 nodes to ∼376,000 nodes created a change in the effective thermal conductivity result of only less than 3.5%. Consequently, in order to maintain computational efforts within reasonable bounds, mesh generation was controlled to remain below 400,000 nodes.

Figure 9.

Convergence study for composites with 4% filler volume fraction.

Figure 10 shows an example of a model with randomly distributed and randomly oriented particles. Ideally, particle spacial distributions for this case should result in isotropic thermal conductivity properties. To further investigate this hypothesis, a study was performed in which effective thermal conductivities were computed along the global RVE directions, i.e., X,Y,Z coordinates. This analysis was completed for the chosen set of filler volume fractions (2.0, 4.0, 7.5 and 10%). Corresponding effective thermal conductivity results for a single model are included in Table 2. The given data indicates that thermal conductivity values were essentially isotropic for models with lower filler loading (i.e., 2%) while for increasing filler volume fractions a mild level of anisotropy was sometimes observed. For that reason, the average effective thermal conductivity was computed from the three Cartesian coordinate directions and subsequently used for comparing composites with randomly oriented particles with aligned filler composites.

Figure 10.

Example of randomly distributed and randomly oriented filler particles.

Filler volume fraction [%]	2.0%	4.0%	7.5%	10%
Thermal conductivity, X; [W m^−1 K^−1]	0.556	0.781	1.08	1.29
Thermal conductivity, Y; [W m^−1 K^−1]	0.525	0.676	1.14	1.26
Thermal conductivity, Z; [W m^−1 K^−1]	0.525	0.778	1.13	1.14
Average thermal conductivity; [W m^−1 K^−1]	0.536	0.745	1.12	1.23

Table 2.

Average and directional thermal conductivities of composites with randomly distributed and randomly oriented filler particles.

As explained previously, the MCS module of the modeling framework performs numerous iterations for each filler volume fraction and subsequently computes the effective thermal conductivity and stores these data in the database. The MCS module repeats this process until specified acceptance criteria are satisfied. For the presented study, the analysis process was terminated after 100 iterations for each of the set filler volume fractions. (Alternatively, a threshold for the unbiased standard deviation or variance could be defined as a termination criterion.) The effective thermal conductivity for a certain filler volume fraction was then computed from the mean of the results stored in the database. Statistical analyses were also performed on the data in order to assess the quality of the employed stochastic process. Data plots for specific volume fractions suggest that data is normally distributed, as shown in Figure 11 by the normalized probability density of effective thermal conductivity data for the direction lateral to filler alignment in a composite with 4.0% filler volume fraction. Data were computed for normality tests for each volume fraction, including the data mean, median, skewness and kurtosis. Tables 3 and 4 list corresponding results for the longitudinal and lateral directions of composites with randomly distributed and aligned particles. These results indicate that the effective thermal conductivity data obey normal distributions, e.g., mean and median of effective thermal conductivity results were practically identical (differences are less than 0.04%). Similar to results presented in [27], appropriate randomness of computed data was thus ascertained. Graphs showing normally distributed effective thermal conductivity data for the chosen filler volume fractions are plotted for the longitudinal and lateral case in Figures 12 and 13, respectively.

Figure 11.

Normalized probability density of effective thermal conductivity data for the direction lateral to filler alignment and 4.0% filler volume fraction.

Filler volume fraction [%]	2.0	4.0	7.5	10.0
Number of iterations; [/]	100	100	100	100
Mean value; [W m^−1 K^−1]	0.775	1.027	1.318	1.489
Median value; [W m^−1 K^−1]	0.774	1.031	1.325	1.486
Standard deviation; [W m^−1 K^−1]	0.019	0.032	0.037	0.038
Variance; [W m^−1 K^−1]	3.7 × 10^−4	1.0 × 10^−3	1.3 × 10^−3	1.4 × 10^−3
Skewness; [/]	1.197	−0.159	−0.672	0.083
Kurtosis; [/]	−0.079	−0.539	0.607	0.128
95% confidence value; [W m^−1 K^−1]	0.0038	0.0064	0.0073	0.0076

Table 3.

Results for longitudinal effective thermal conductivity and statistical analyses.

Filler volume fraction [%]	2.0	4.0	7.5	10.0
Number of iterations; [/]	100	100	100	100
Mean value; [W m^−1 K^−1]	0.262	0.276	0.306	0.332
Median value; [W m^−1 K^−1]	0.262	0.276	0.306	0.331
Standard deviation; [W m^−1 K^−1]	0.001	0.001	0.003	0.004
Variance; [W m^−1 K^−1]	1.6 × 10^−6	2.6 × 10^−6	9.4 × 10^−6	2.3 × 10^−5
Skewness [/]	−0.819	−0.267	−0.180	0.216
Kurtosis [/]	2.410	0.031	−0.124	0.344
95% confidence value; [W m^−1 K^−1]	0.0002	0.0003	0.0006	0.0009

Table 4.

Results for lateral effective thermal conductivity and statistical analyses.

Figure 12.

Normalized probability density of effective thermal conductivity data for the direction of filler alignment and different filler volume fractions (VF).

Figure 13.

Normalized probability density of effective thermal conductivity data for the direction lateral to filler alignment and different filler volume fractions (VF).

The effective thermal conductivity data calculated by SFEA framework can be considered continuous random variables, and hence, it is recommended to calculate the probability of occurrence of an explicit effective thermal conductivity within an identified interval. This calculation can be performed using Eq. (3).

$$P\left(a\le {\rm X}\le b\right)={\int }_a^{b}f\left(\chi \right)\mathit{d\chi }$$ E3

where $P$ is the probability of an event of explicit effective thermal conductivity within the interval $a$ and $b$ ; ${\rm X}$ and $f\left(\chi \right)$ are correspondingly a continuous random variable and the probability distribution function. A cumulative distribution function (CDF) can be computed from Eq. (3) for each of the selected filler volume fractions. Corresponding CDF graphs are depicted in Figures 14 and 15 for the longitudinal and lateral cases of aligned filler composites, respectively.

Figure 14.

Cumulative density function of effective thermal conductivity data for the direction of filler alignment and different filler volume fractions (VF).

Figure 15.

Cumulative density function of effective thermal conductivity data for the direction lateral to filler alignment and different filler volume fractions (VF).

Finally, the modeling approach was used to achieve the objective of the study, that is, assessing the effect of filler addition and alignment on heat transfer into polymer composites in the context of fire-retardancy. Figure 16 depicts average effective thermal conductivity results for different filler volume fractions for the cases of randomly oriented (isotropic) and aligned rod-shaped particles mimicking CNT. For the aligned filler morphologies, effective thermal conductivity results are shown for the direction of filler alignment (longitudinal) and the corresponding lateral direction. Despite the fact that the modeling approach employed generic material properties and simplifying assumptions for the CNT geometry, the data is in satisfactory agreement with data published in the technical literature (e.g., [32]). Experimentally characterized CNT-polymer composites involve a wide range of material properties and fabrication routes. Notwithstanding these differences, the modeled thermal conductivity in the order of 1.0 W m⁻¹ K⁻¹ for isotropic composites with filler loadings approaching 10% are shown to be realistic.

Figure 16.

Effective thermal conductivity of randomly oriented and aligned rod-shaped particles embedded in epoxy polymer matrix.

The relative increase in effective thermal conductivity for the different filler volume fractions is depicted in the graph in Figure 17 for the cases of randomly oriented (isotropic) and aligned particles (longitudinal and lateral). These data clearly demonstrate that heat transfer into the polymer can greatly be reduced when filler particles are aligned parallel to the surface of a component. For example, while thermal conductivity in an isotropic and aligned filler composite was found to respectively increase almost sixfold and fivefold over the matrix for 10% filler loading, the lateral thermal conductivity in the aligned filler composite rose only by a factor of 1.3. In the context of fire-retardancy it can therefore be concluded that aligning CNT and other high aspect ratio carbon allotrope fillers parallel to the surface of a polymer component may provide an effective means for alleviating heat input into the material while enabling the desired mass transport barrier for mitigating fuel release into the gas phase, and reduced peak heat release and radiant heat flux.

Figure 17.

Relative increase in effective thermal conductivity of randomly oriented and aligned rod-shaped particles embedded in epoxy polymer matrix.

5. Conclusions

A stochastic finite element analysis framework was employed to simulate the effective thermal conductivity of randomly distributed rod-shape particles mimicking carbon nanotubes embedded in a polymer matrix. Particles were either randomly oriented or aligned, creating isotropic or anisotropic thermal conductivity behavior, respectively. The modeling framework that is based on Monte Carlo simulation considers filler-matrix and particle-to-particle interfacial effects.

The numerical study indicated that the effective thermal conductivity is greatly enhanced for aligned filler composites in the alignment direction and isotropic filler modified composites. However, in the direction lateral to filler alignment the increase in thermal conductivity is only modest. Therefore, in order to limit heat input into the material, CNT and other high aspect ratio carbon allotrope fillers may be aligned parallel to the surface of a polymer component. In this manner, the flame-retardancy effectiveness of filler modified polymer composites can be improved, while providing a mass transport barrier that lessens the release of fuel into the gas phase, peak heat release and radiant heat flux, all of which were previously described in the technical literature.

Acknowledgments

The authors gratefully acknowledge Hi-Kalibre Equipment Ltd. for providing software and computing resources for this research. This work was further supported by the Natural Sciences and Engineering Research Council of Canada, Grant number: RGPIN-2016-04650.