1 Introduction

Fibre-reinforced polymer (FRP) composite laminates have been widely used in aerospace, automotive, and wind energy industry due to their excellent material properties such as high stiffness-to-mass ratio, high strength, and light weight. Applications of FRP composite laminates to create engineering structure models fundamentally require mechanical properties as inputs. Experimental tests are ideal solutions to evaluate the mechanical properties of a composite lamina. However, it must be repeated whenever the constituents (fibre and matrix) and/or microstructure characteristics (fibre volume fraction) are altered. This procedure may, for instance, costs millions of dollars and lasts for years to generate the experimental data of mechanical properties for the design of aircraft structures [32].

To overcome the aforementioned drawbacks associated with experimental tests, various micromechanical approaches have been proposed to establish a closed-form relationship between elastic properties at the lamina scale and the elastic properties at the constituent scale. These methods fall generally into two categories, i.e., analytical and numerical methods. Analytical methods, such as the Rule of Mixture method [6], the Halpin–Tsai semi-empirical method [10], the Mori–Tanaka method [19], and the Chamis method [5], facilitate the calculation of elastic properties by a direct mathematical, empirical expression between constituent properties and elastic properties of the lamina. However, these methods have an inherent limitation in describing the stress and strain fields at microstructural scale mainly due to neglecting fibre interaction.

With the development of computing capacities, numerical methods, in particular the finite element method (FEM), have become widely used tools for studying the behaviour of composites, including inverse analysis [8, 15, 26], elastic moduli [27], failure of composite lamina [29, 33, 34], and the effective coefficients of thermal expansion [13]. Inverse analysis has been used to identify fibre mechanical properties and fibre thermal expansion coefficient, and evaluate the factors of analytical methods, and so on. [26] determined the elastic and thermal properties of graphite fibre using inverse analysis. [3] predicted fibre properties using finite element analysis of hexagonal and random representative volume element (RVE) through inverse analysis. Similarly, [14] conducted an inverse analysis to predict fibre mechanical properties. However, they used quasi-analytical gradients derived from analytical models such as Chamis or Halpin–Tsai to reduce the computational cost. [15] utilised an inverse method to identify the mechanical properties of T300 carbon fibre as well as the interphase region parameters based on a computational homogenisation approach together with experimental results and Kriging metamodelling. [8] carried out an inverse analysis in the framework of FEM to estimate the reinforcement parameter ξ of the Halpin–Tsai models which is used to calculate transverse stiffness E2. A total number of 67 FE models of 2D square, 2D hexagonal, and 3D random fibre distributions were used to obtain a new value of ξ with a high level of confidence.

Regardless of the inverse methods used or the purpose they are used for, a large number of FE analyses are required for converged solutions. However, constructing a micromechanical FE model is not a straightforward task and requires special treatments to impose boundary conditions, and generate microstructures including fibre distribution and fibre/matrix interphase, and extract outputs, etc. These complexities impose a barrier of using inverse analysis by engineers and researchers. Recently, several ABAQUS plug-ins have been developed for the ease of creating micromechanical FE models. These plug-ins were developed either by the functions available in ABAQUS or by external software. An ABAQUS plug-in named MultiMech was developed to perform multi-scale finite-element analysis (FEA) with the capability of simulating nonlinear behaviours of composites [16]. Another ABAQUS plug-in for multilevel modelling of linear and nonlinear behaviour of composite structures [7, 28]. The plug-in developed using Python scripts for analysing an RVE at microscopic level to obtain macroscopic parameters for structural analysis by user-defined FORTRAN subroutines in ABAQUS. Composite MicroMechanics (COMM) toolbox was developed in Matlab for micromechanical analysis of composites [17]. The toolbox creates an input file that can be read by ABAQUS which performs the FEA. Recently, EasyPBC plug-in was developed for ABAQUS to estimate effective elastic properties of a pre-prepared and meshed RVE [21]. While, the ABAQUS plug-in proposed by [24] is capable of generating an RVE with random fibre distribution using Random Sequential Adsorption (RSA) technique.

The aforementioned plug-ins have shown outstanding benefits and capabilities to create and simulate complex RVEs of unidirectional (UD) FRP composite lamina. However, they are designed to generate and analyse a single model. Therefore, this paper aims to develop an open-source ABAQUS plug-in named Virtual Data Generator (VDGen) that automates the time-consuming manual task requires to create a large number of virtual data for inverse analysis. The plug-in uses Latin-Hypercube (LH) sampling methods and supports the unit cell of square and hexagonal fibre arrays. In addition, the plug-in incorporates Artificial Neural Networks (ANN) to explicitly parameterise the relationship between fibre mechanical properties and fibre/matrix interphase parameters and the mechanical properties of a UD lamina. The data required here were created in advance by the plug-in and used to train the ANN model.

2 Main plug-in GUI

The concept of the plug-in arises from the need for a tool that helps to perform a large number of micromechanical FE simulations in a few simple steps. ABAQUS has different ways to increase its capabilities such as subroutines and/or adding new plug-ins. ABAQUS/CAE plug-in is one of the most powerful tools that can be used to perform pre- and post-processing via functions written in Python programming language in the kernel. The current plug-in operates through a series of user-friendly GUI commands send to the kernel to carry out tasks. The plug-in interface is shown in Fig. 1. It consists of six tab items that allow the user to navigate between them to edit input and output commands. For computational micromechanics modelling, the plug-in supports square and hexagonal unit cell fibre arrays. Despite fibres are usually randomly distributed in the matrix, it has been concluded by [31] that micromechanical modelling of the unit cell is accurate enough to predict the elastic properties of a UD lamina, while an RVE with randomly distributed fibres is essential to compute the local failure.

Fig. 1
figure 1

The graphical user interface (GUI) of the plug-in

Full size image

Figure 2 shows a typical 3D unit cell of square and hexagonal fibre arrays of the fibre reinforced composite that the plug-in supports. The mechanical properties of each constituent, i.e., fibre, matrix and interphase, can be modified in the material section. There are two ways to input the value of constituent properties, either by a single value or using a domain of lower and upper bounds (lower–upper). A material property assigned with a single value remains unchanged throughout simulations. While for others, a random value in the range of (lower–upper) is generated at each training point using Latin-Hypercube sampling technique.

Fig. 2
figure 2

3D unit cell models of regular fibre arrays: a Square and b Hexagonal

Full size image

The fibres and matrix are meshed using eight-node brick element with reduced integration (C3D8R). There were also a relatively small amount of six-node linear triangular prism elements (C3D6) due to the free meshing technique used. The interphase region is meshed with eight-node cohesive elements (COH3D8). To maintain matched meshes between the cohesive elements and the fibres and matrix elements, a suitable number of nodes are seeded to the interphase and its neighbours. The elastic behaviour of the cohesive elements is written in terms of a stiffness matrix that relates the nominal stresses to the nominal strains across the interphase. The nominal traction stress vector t consists of three components, tn, ts, tt, which represent the normal and two shear tractions, respectively. The corresponding separations are denoted by δn, δs, and δt, and the original thickness of the cohesive element is denoted by T. Then, the nominal strains can be defined as

$${\varepsilon }_{n}=\frac{{\delta }_{n}}{T}, {\varepsilon }_{s}=\frac{{\delta }_{s}}{T}, {\varepsilon }_{t}=\frac{{\delta }_{t}}{T}.$$
(1)

Therefore, the elastic behaviour of the cohesive element can be written in Eq. (1). For simplicity of computation, uncoupled behaviour between the normal and shear components is desired, so the off-diagonal terms in the elasticity matrix are set to be zero and the stiffness in the two shear directions are assumed to be equal [1, 33]

$$t=\left\{\begin{array}{c}{t}_{n}\\ {t}_{s}\\ {t}_{t}\end{array}\right\}=\left[\begin{array}{ccc}{K}_{nn}& {K}_{ns}& {K}_{nt}\\ {K}_{ns}& {K}_{ss}& {K}_{st}\\ {K}_{nt}& {K}_{st}& {K}_{tt}\end{array} \right]\left\{\begin{array}{c}{\varepsilon }_{n}\\ {\varepsilon }_{s}\\ {\varepsilon }_{t}\end{array}\right\}=K\varepsilon .$$
(2)

The plug-in imposes Periodic Boundary Conditions (PBC) on the corresponding surfaces of the unit cell to ensure the compatibility of strain and stress at the macroscale level. These consist of a series of constraints in which the deformation of each pair of nodes on the opposite surfaces of the unit cell is subject to the same amount of displacements. The PBCs are expressed in terms of the displacement vectors \({\overrightarrow{U}}_{1}\), \({\overrightarrow{U}}_{2}\), and \({\overrightarrow{U}}_{3}\) that are related to the displacements between the opposite surfaces by

$$\left\{\begin{array}{c}\overrightarrow{u}\left(0,{x}_{2},{x}_{3}\right)-\overrightarrow{u}\left({L}_{1},{x}_{2},{x}_{3}\right)={\overrightarrow{U}}_{1}\\ \overrightarrow{u}\left({x}_{1},0,{x}_{3}\right)-\overrightarrow{u}\left({x}_{1},{L}_{2},{x}_{3}\right)={\overrightarrow{U}}_{2}\\ \overrightarrow{u}\left({x}_{1},{x}_{2},0\right)-\overrightarrow{u}\left({x}_{1},{x}_{2},{L}_{3}\right)={\overrightarrow{U}}_{3}\end{array}\right.,$$
(3)

where L1, L2, and L3 are the lengths of the unit cell along with three orthogonal directions, respectively. PBC requires matching nodes on opposite sides of the unit cell. Hence, elements of equal size are assigned to the edges of the unit cell to ensure periodic mesh required for PBC.

The Output tab allows the user to select appropriate results that suit the work. The macroscopic normal and shear strain components are calculated by

$${\varepsilon }_{ii}=\frac{{\overrightarrow{U}}_{i}}{{L}_{i}}\;\left(\mathrm{normal\; strain}\right),$$
(4)
$${\varepsilon }_{ij}=\frac{{\overrightarrow{U}}_{i}}{{L}_{j}}+\frac{{\overrightarrow{U}}_{j}}{{L}_{i}}\;(i\ne j)\;(\mathrm{shear\; strain}).$$
(5)

The macroscopic stress is calculated as

$${\sigma }_{ij}=\frac{\sum {F}_{i}}{{A}_{j}},$$
(6)

where Fi is the resultant force on the ith surface which represents the reaction force at a reference point where the displacement is applied, and A is the area of the surface. Therefore, Young’s modulus, Poisson’s ratio, and shear modulus are, respectively, calculated from Eqs. (7), (8), and (9)

$${E}_{ii}=\frac{{\sigma }_{ii}}{{\varepsilon }_{ii}}$$
(7)
$${G}_{ij}=\frac{{\sigma }_{ij}}{{\varepsilon }_{ij}}$$
(8)
$${\nu }_{ij}=-\frac{{\varepsilon }_{ij}}{{\varepsilon }_{ii}}.$$
(9)

The flowchart of the pre‑ and post-processing procedure of the plug-in is described in Fig. 3. The user is to define the analysis of data and the required outputs, as given in Fig. 3. When the ‘OK’ or ‘Apply’ button is clicked, the plug-in creates three files to be called in the subsequent steps. ExperimentNAME.dat file contains all input commands which are given by the user in the GUI interface. These commands are vital to creating the FE model. Also, this file provides an opportunity to modify the input commands using exact plug-in keywords in the file. ~ Temp.txt is dedicated to storing the experiment name the user wants to run. Material properties values created by Latin-hypercube sampling are stored in a tabular format in ExperimentNAME.csv file, which makes them easier to read and process by external software. Once these files are created, the Python script (Execute.py) can be submitted for analysis from ABAQUS command environment using either dos (‘abaqus cae script = Execute.py’) or dos (‘abaqus cae noGUI = Execute.py’) command. It is strongly recommended to execute it using the latter in which ABAQUS/CAE runs commands in Execute.py without the added expense of running a GUI display. Whichever option adopted to perform the analysis, the plug-in continuously provides the user with useful information, e.g., number of jobs done, number of jobs remain and approximate time to complete. Upon running, the plug-in instantly creates two files to store outputs after completion of each job and to record errors that occurred during the analysis.

Fig. 3
figure 3

Flowchart of the plug-in

Full size image

3 Numerical example: prediction of effective elastic properties

To validate the newly developed plug-in, the effective elastic properties of carbon fibre/epoxy (T300/PR-319) and glass fibre/epoxy (E-Glass/MY750) were determined and compared with EasyPBC plug-in developed by [21] as well as the experimental data [12]. The mechanical properties of the fibre, matrix, and interphase are given in Table1. Since T300 carbon fibre is classified as a transversely isotropic material, the elastic properties highlighted by an asterisk * symbol in the table are obtained by applying the following relations:

Table 1 Elastic properties of fibre, matrix [12] and interphase [15]
Full size table
$${G}_{f12}={G}_{f13}, { G}_{f23}=\frac{{E}_{f2}}{2(1+{\nu }_{f23})}.$$
(10)

The E-glass is considered as isotropic material. It is important to note that the input data for the interphase region are not accurately known as they are difficult to measure from simple laboratory experiments. However, an initial stiffness Ki of 105 GPa/mm is used in [2, 18, 25, 30] to simulate the elastic behaviour of the RVE model. In this paper, the elastic parameters from [15] are used to for the interphase as an approximation.

Table 2 shows the comparison of the predicted effective elastic properties determined by VDGen and EasyPBC plug-ins. It is noted among the prediction results that VDGen provides reliable results that are identical to those from EasyPBC. However, an obvious discrepancy exists between experimental results and those predicted by the two plug-ins. This is mainly due to the inaccurate parameters used for the interphase region.

Table 2 Comparison of predicted effective elastic properties determined by VDGen, EasyPBC and experimental data [12, 15]
Full size table

Figure 4a shows the stress contours of a loaded unit cell under transverse and Fig. 4b shows the stress contours for in-plane shear loading conditions. It can be seen the periodic stress contours distribution which is additional verification of the PBC.

Fig. 4
figure 4

Stress contour plot of the unit cell under a Transverse loading and b In-plane shear loading

Full size image

4 Application: identifying fibre and interphase parameters using ANN

4.1 Background

In this section, a machine learning (ML) technique Artificial Neural Network (ANN) is used to construct a relationship between fibres and interphase parameters and effective elastic properties of the lamina. It is inspired by the animal brain’s structure and function, which learns from former examples. ANN consists of three main layers: an input layer, one or more hidden layers, and an output layer. Each layer has several neurons which are responsible for transmitting weight and biases (equivalent to chemical and electric signals in the animal brain) between two layers. Figure 5 illustrates a typical structure of single neurons where each input (x) comes from the previous layer multiplied with its individual weight of the connection (wi) and then summed up with biases [23]. Then, this sum is composed with activation function (f), resulting in another vector (a) as

Fig. 5
figure 5

A typical one neuron structure of the ANN

Full size image
$$a=f\left({x}_{i}.{w}_{i}+b\right).$$
(11)

Another key step of ANN is a defined objective function that is to be minimised during the training process. Mean Squared Error (MSE) and Sum Squared Error (SSE) among others are examples of functions used to assess the network’s behaviour by measuring the errors between the output and the target. The errors are reduced through tuning the values of the weight and biases by the so-called back-propagation. Back-propagation is widely recognised as a powerful tool for the training of ANN very efficiently. Several algorithms have been proposed to address the slow convergence associated with the back-propagation. However, it is quite difficult to decide which algorithm is more computationally efficient as it depends on many factors. Readers may refer to a comparative study carried out by researchers to evaluate accuracy and convergence time for different algorithms [4, 9].

4.2 ANN model to identify micro-parameters

ANN model is developed to identify the micro-parameters, e.g., fibre and fibre/matrix interphase parameters. The relationship between micro- and macro-properties in UD lamina is fairly complex and nonlinear. Moreover, the number of micro-parameters to be identified is usually more than the number of macro-properties, which makes the ANN a complicated task. Therefore, to ease the training process, the micro-parameters are set to be the input layer of the neural networks model and the macro-properties are of the output layer. However, the calculation of optimal micro-parameters becomes difficult when they are in the input layer as it is not possible to obtain an analytical inverse response solution with the ANN model that has multiple neurons in the hidden layer. This issue is overcome using trained ANN to enlarge the dataset. Details of model building are explained in the following section.

4.2.1 Model building

The whole procedure of the fibre and interphase parameters identification using ANN is illustrated in Fig. 6. Firstly, a total of n = 500 FE models were created by VDGen using the procedure outlined in Sect. 3. In each model, a random value of the parameters to be identified is created by LH sampling within the range given in Table 3. The remaining fibre properties were obtained by applying the transversely isotropic material relationships

Fig. 6
figure 6

Fibre and interphase parameters identification through ANN

Full size image
Table 3 Ranges of fibre and interphase elastic properties
Full size table
$${E}_{f3}={E}_{f2}, {v}_{f13}={v}_{12f},{{ G}_{f13}={G}_{f12}, G}_{f23}=\frac{{E}_{f2}}{2(1+{\nu }_{f23})}.$$
(12)

Ef1 remained unchanged in all samples and its value was 230GPa. For all samples, the fibre volume fraction is 60% and matrix properties are given in Table1. By the end of this phase (Step 1), a dataset of 500 samples containing the inputs (x = [Ef2, νf12, νf23, Gf12, Ti, Knn, Kss]) and the targets (t = [E11, E22, ν12, ν23, G12]) required to train ANN is obtained (Fig. 7).

Fig. 7
figure 7

ANN scheme for fibre and interphase parameters identification

Full size image

In Step 2, the ANN was built and trained using the dataset created in the previous step. The training, testing, and validation of the ANN model were conducted by MATLAB R2015a software. By MATLAB default, 70%, 15%, and 15% of the original dataset were used for training, validation, and testing, respectively. Nonlinear tangent sigmoid and linear functions were employed as the activation functions in the hidden layers and the output layer, respectively.

$$\begin{gathered} f\left( x \right) = 2/\left( {1 + e^{ - 2x} } \right) - 1\quad tansig\left( {\text{tangent sigmoid activation function}} \right) \hfill \\ f\left( x \right) = x\quad purelin\left( {\text{linear activation function}} \right). \hfill \\ \end{gathered}$$
(13)

Selecting a best representative ANN structure plays an important role in output prediction. In this study, two hidden layers were used and the number of neurons in each hidden layer was changed until the best possible prediction was obtained. Initially, the number of neurons in the first hidden layer (nh1) was set to 20 then increased by one, whereas the total number of neurons in both hidden layers (nh) was retained at 100. Usually, the data are randomly divided into three subsets (training, validation, and testing), and different initial weight and bias values are used in each time the neural networks are trained. As a result, different neural networks trained for the same problem may give different outputs for the same inputs. In this study, therefore, 20 runs were performed on each ANN architecture to ensure inclusion of different data for each subset. MSE, which is the average squared difference between the output (y) vectors and the target (t) vectors, was used to compute the difference and back-propagated though the networks to update weights and biases

$$MSE=\frac{1}{n}\sum_{i=1}^{n}{\left(y-t\right)}^{2}.$$
(14)

Levenberg–Marquardt (LM) back-propagation algorithm was adopted in this study. This function uses back-propagation scheme to update weights and biases according to Levenberg–Marquardt optimization algorithm, which can accurately achieve results with fewer data comparing with its counterparts.

To decide which is the best ANN taking into account the fact of dividing the input data into three main subsets (training, validation, and testing) during the training process, the sum of correlation coefficient (R-value) between the target and the output of the entire dataset was used to attain the optimal ANN structure

$$\mathrm{max}\,R(x)={R}_{v}^{E11}+{R}_{v}^{E22}+{R}_{v}^{v12}+{R}_{v}^{v23}+{R}_{v}^{G12},$$
(15)

where \({R}_{v}^{E11}, {R}_{v}^{E22},{R}_{v}^{v12}, {R}_{v}^{v23},\mathrm{ and}\, {R}_{v}^{G12}\) are the regressions of the dataset of E11, E22, ν12, ν23, and G12, respectively. The optimal ANN was then used to extend the dataset and creating new N samples required to accurately identify fibre and interphase parameters. This was conducted using LHS to create N random samples of the inputs within the range given in Table 3. These samples are then passed through the trained ANN to obtain the output. LH sampling method is appropriate for this work as it ensured a full coverage of the input sample space. A denser space is constructed by selection more points, which results in more reliable results. The optimum carbon fibre and interphase parameters were selected based on the smallest value obtained from Eq. (16)

$$\mathrm{min}\,F(x)=\frac{\left|{E}_{11}^{t}-{E}_{11}^{ANN}\right|}{{E}_{11}^{t}}+\frac{\left|{E}_{22}^{t}-{E}_{22}^{ANN}\right|}{{E}_{22}^{t}}+\frac{\left|{v}_{12}^{t}-{v}_{12}^{ANN}\right|}{{v}_{12}^{t}}+\frac{\left|{v}_{23}^{t}-{v}_{23}^{ANN}\right|}{{v}_{23}^{t}}+\frac{\left|{G}_{12}^{t}-{G}_{12}^{ANN}\right|}{{G}_{12}^{t}},$$
(16)

where \({{E}_{11}^{t}, E}_{22}^{t}, {v}_{12}^{t}, {v}_{23}^{t}\mathrm{ and }\,{G}_{12}^{t}\) are the experimental effective elastic properties given in Table 2. \({{E}_{11}^{ANN}, E}_{22}^{ANN}, {v}_{12}^{ANN}, {v}_{23}^{ANN},\mathrm{ and }{G}_{12}^{ANN}\) are corresponding effective elastic properties predicted by the ANN.

Finally, the identified parameters were used to calculate the effective elastic properties using the plug-in (Step 3) and results were compared with experimental data.

4.3 Results and discussion

The ANN which can predict fibre and fibre/matrix interphase parameters from micromechanical FE modelling dataset is designed. The number of total samples created by micromechanical FE modelling to train the ANN is 500. 350 of them is randomly assigned for training, 75 for validation set and the rest used for testing. The ANN is built and validated as explained in Sect. 4.2.1.

Table 4 presents some architecture samples used to verify the performance of the ANN in terms of R-value. Since 20 runs are performed within each ANN structure to ensure the use of different data for each subset and due to the space restrictions of the paper, the table shows only the best performance run from the 20 runs. It can be seen from the table that the maximum R-value is attained when 42 neurons at the first hidden layer and 48 neurons at the second hidden layer are used.

Table 4 Verification cases of the ANN architecture in terms of R-value for each output element
Full size table

Figure 8 shows the regression graphs for the target and the output of the verification data set only. This figure shows the closeness among the output data predicted by the ANN and the target data obtained from FEM. The dashed line in each subfigure represents the perfect results when the outputs equal targets. It can be seen that E11, E22, ν12 and ν23 are well predicted by ANN with an R-value between 0.91 and 0.99 and a regression slope (m) between 0.81 and 0.96. G12 is slightly less well predicted comparing to other effective elastic properties with an R-value and a regression slope of 0.86 and 1.08, respectively. Hence, the selected ANN is capable of providing a good correlation between the target and the output.

Fig. 8
figure 8

ANN predictions versus FEM results for a E11, b E22, c ν12, d ν23, and e G12

Full size image

After training the ANN, the selected model with the highest R-value is used to generate N samples. It is found that 10,000 samples of the random input parameters generated by the LH sampling method are sufficient to produce a dense space. These new input data (Ef2, νf12, νf23, Gf12, Ti, Knn and Kss) are then processed by the ANN to obtain and the outputs (E11, E22, ν12, ν23, and G12). The closest point of the new output to the experimental data is found through Eq. (16). The corresponding carbon fibre and interphase parameters of the closest point are given in Table 5.

Table 5 Identified parameters for fibre and interphase region
Full size table

Finally, the identified fibre and interphase parameters (micro-parameters) obtained from the ANN are used as input for the FE model to predict the effective properties of the UD lamina, i.e., macroscale level properties. The effective elastic properties calculated by the FE model using the identified parameters are given in Table 6 (second column). The table also shows a comparison between these properties and those predicted by the ANN and the experimental data. It can be seen that the effective elastic properties agree well with the experimental data with a maximum error of about 6%. This error is mainly due to the using fixed value for E11 in the training of ANN.

Table 6 Macro-level properties predicted by ANN and FE model using identified parameters and experimental data [15]
Full size table

5 Conclusions and future improvement

An ABAQUS® plug-in (VDGen) has been developed for generating virtual data for identifying the uncertain materials’ properties in unidirectional (UD) lamina. In combination with artificial neural networks (ANNs), the data generated from the plug-in enable the determination of the relationship between fibre mechanical properties and fibre/matrix interphase parameters at microscale and the mechanical properties of a UD lamina at macroscale. Application of the plug-in to a T300/PR-319 UD lamina has shown very good agreement between the predictions and the experimental data when using the identified constituent properties.

A few improvements of the plug-in should be considered in the future. The current plug-in is designed to support the square and hexagonal unit cell fibre arrays. Micro-mechanical FE modelling of randomly distributed fibres in the matrix is essential when studying the failure of the composite lamina. However, using random fibre distribution causes arbitrary meshing condition on opposite RVE edges [20, 35]. Further development of the plug-in to support RVEs with randomly fibre distribution and capable of generating periodic mesh on the opposite edges is under investigation by the authors.

At this stage, the plug-in is only designed to calculate the effective elastic properties of a lamina. We aim to develop it further, so that it will be capable of conducting failure analysis under uniaxial, biaxial, and multiaxial loading conditions.

The effect of fibre shape has recently been subjected to intensive studies by means of computational micromechanics [11, 22]. The current core Python scripts of the plug-in will be developed further, so that RVEs with different fibre shapes can be automatically generated.