Ef ﬁ cient meshing technique for textile composites unit cells of arbitrary complexity

Meso ‐ scale unit cell models are often used to simulate mechanical behaviour of textile composites. Apart from reliable ways to create meso ‐ scale geometries, such simulations require reliable meshing algorithms. While the former is made possible via dedicated textile pre ‐ processors or high ‐ ﬁ delity weaving simulations, the meshing remains quite problematic for complex textiles and geometries. Even though, with a lot of user input, it is possible to create very complex meshes using meshing pre ‐ processors, this approach remains infeasible for cases when a large number of models need to be analysed. This paper presents a meshing approach based on the combination of local octree ‐ re ﬁ nement with surface smoothing. This allows nearly conformal meshes to be generated for geometries of any complexity which achieve accuracy comparable to that of conformal meshes. A range of unit cells was analysed using the new approach and it was shown that the error in local stresses is within 10% of the reference solution and the average error is below 7%. It was found that the computational cost of the analysis using the new meshing tech- nique is not considerably higher than for an analysis which uses a conventional conformal mesh yet the new approach allows analysis of any geometry.


Introduction
Most of the numerical methods for modelling of composites with textile reinforcements rely on meso-scale modelling techniques. Therefore, generation of finite element (FE) meshes is an essential step in these procedures. Generic problems, arising during mesh generation, have been described in one of the seminal books on tetrahedral meshing [1] (all-hexahedral meshing still remains a very challenging problem [2]) and can be summarised as the requirement on elements to respect internal and external boundaries, have good shape and not to be too numerous. These requirements on the meshes are closely linked to the initial geometry and its features such as curvature of surfaces, closeness of the internal and external surfaces to each other and the overall size of the geometry. Meso-scale geometries of textile composites contain two of these features: high double curvature of yarn surfaces and yarns being close to or touching each other as shown in Fig. 1. The latter often triggers excessive local mesh refinement and therefore leads to a large number of elements. The aforementioned problems and also this paper focus only on the meshing itself and leave aside problems such as generation of CAD geometry and its accuracy and consistency.
A new technique presented in this paper can create a nearconformal mesh for complex geometries ensuring control over the minimum/maximum size of elements as well as their total number. The algorithm, which employs local octree refinement of voxel meshes as well as mesh smoothing, is built into the TexGen textile preprocessor [4] which makes it immediately available for many types of geometries and performs automatic assignment of local fibre orientation and local fibre volume fraction within the yarns. An example of a mesh generated with this algorithm is shown in Fig. 2.
In principle, a tetrahedral mesh generation algorithm can generate a mesh for almost any geometry. In practice, the meshing algorithms often struggle to mesh matrix domains which occur in modelling of textile composites, in particular narrow regions between two yarns, because of the limitations on the minimum size of element or on the total number of elements (from either memory restrictions or solution time). Setting a mesh tolerance equal to the minimum distance between the yarns does not always solve the problem as it can create volumes of odd shape. Alternatively, manual interventions in the meshing are not desirable as it makes the approach impractical when a large number of models need to be analysed. The meshing of textile composite unit cells has been addressed with various approaches.
Early research in modelling of textile composites used relatively simple geometry with no matrix pockets which made mesh generation relatively simple [5,6]. The next generation of models [7][8][9] was based on more realistic models which mainly generated the yarn path using a more complex curve than a trigonometric or piece-wise function. This geometry very often gave complex matrix pockets in between the yarns which were often enlarged to make the meshing easier [8,9]. However, introduction of these pockets into the model makes it less realistic not only because of the changed geometry but also because of the reduction of the yarn volume which needs to be compensated by increase of intra-yarn fibre volume fraction. These difficulties with meshing of complex textile geometries were addressed by several approaches: voxel meshing [10], mesh-free methods [11], embedded meshes (also known as the mesh superposition method) [12] as well as various manipulations with geometry to achieve a high-quality conformal mesh [13,14].
The embedded mesh approach allows the matrix mesh to be as simple as a mesh of a domain with no inclusions. A mesh of yarns is gen-erated using a simple mesh sweeping method and then linked to the domain mesh through a set of linear constraints for the nodes of both meshes. Ease of meshing comes at a cost of stress and strain discontinuity at the boundary of the two meshes and absence of an actual surface between a yarn and matrix.
Generation of high-quality conformal meshes requires the issue with the geometry in the zone where yarns approach each other to be addressed. Grail et al. [13] employed textile compaction simulations to eliminate the gap between yarns and then re-meshed each of the yarns to achieve a conformal mesh. A similar approach was used by Ha et al. [14] who generated the geometric model with allowed intersections which were then eliminated using a special algorithm to create a contact area between yarns. This contact area ensures that the yarn meshes are conformal in this area and that meshes do not contain ill-shaped elements. This geometry preparation procedure applied to a textile with yarns of regular elliptical cross-section resulted in a matrix domain which can be meshed with tetrahedral elements. In summary, both of the approaches rely on a tetrahedral meshing algo-  rithm to generate a tetrahedral mesh using surface meshes of yarns and periodic meshes for domain boundaries as input meshes. While this works in many cases, in the authors' experience, such software as Hypermesh, a commercial software for mesh generation, and TetGen [15], an open-source software, are not always capable of generating a mesh mainly due to the constraints imposed to create periodic boundary meshes. The voxel mesh approach works relatively well for simple textile geometries where voxel elements are mainly aligned with the yarns and load is applied in the yarn direction [16]. The voxel mesh accuracy degrades for geometries not aligned with the voxel mesh therefore more elements or a smoothing technique [17] are required. Increasing the mesh density in a voxel mesh is very inefficient because of the uniform element size. This was addressed by Kim and Swan [10] who implemented local refinement of voxel meshes for textile composites. It was shown that it is possible to reduce the total number of elements without reducing the accuracy of the FE results. However, this refinement still had a jagged surface and local stress concentrations owing to the nature of the voxel mesh. Other studies on voxel meshes (especially those obtained from medical imaging) employed Laplacian smoothing and other algorithms to create surfaces which would be relatively smooth and have a surface close to the original surface. This paper combines both of these ideaslocal octree refinement to approximate the boundaries of yarns as closely as possible, followed by surface smoothing to reduce the effect of stress concentration. The approach is described in detail in Section 2 which is followed by validation of the technique in Section 3. A particular focus is given in Section 4 on comparison of the proposed technique with conventional conformal meshes in terms of accuracy and computational costs. Summary and the plans for future developments of the method are given in Section 5.  2. Octree mesh refinement and mesh smoothing

Octree refinement of a voxel mesh
Typically mesh refinement is applied to reduce the numerical error introduced by mesh discretisation. Experience shows that surfaces with high curvature need to be discretised with smaller elements otherwise spurious stress concentrations may occur. The refinement employed here aims to represent the geometry as closely as possible but does not aim to have a conformal mesh (where nodes fall exactly on the surfaces of the internal geometry e.g. yarns). The mesh algorithm starts with a simple voxel mesh generated for a unit cell. A refinement algorithm is then applied to each voxel. A voxel is refined by splitting into eight voxels if its nodes belong to different materials. The dimensions of the new voxels are half that of the original voxels as shown in Fig. 3. This process is repeated iteratively until a given level of refinement is achieved.
The mesh obtained using the described refinement procedure can be conveniently stored in an octree data structure, hence its name "octree mesh refinement". An open-source library p4est [18] was employed for this purpose. Storing a refined mesh as an octree also makes it simple to balance the mesh i.e. always have 2:1 transition between elements instead of an arbitrary transition as shown in Fig. 4. The mesh refinement implemented in this paper also ensures that the interface between the materials is always represented by elements with the highest level of refinement.
The octree refinement in this paper was implemented such that both the size of elements in the starting voxel mesh and the refined mesh are proportional to the size of the unit cell as 1/2 n where n is an integer number between the starting level of refinement and the maximum level of refinement. However, this condition is easy to relax if a geometry is split into arbitrary domains which are then individually meshed using the octree technique. Examples of various octree refinements are shown in Fig. 5. The meshes in this paper will be referred to by their minimum and maximum levels of refinement e.g. suffix '4-8' denotes a mesh with minimum level of refinement of 4 and maximum level of refinement of 8.
The octree mesh refinement allows the mesh to be as close to the real interface as needed without creating additional elements away from an interface. For example, achieving a 1/(2 8 ) precision for the geometry shown in Fig. 5 ('5_8' mesh of elliptic inclusion) requires 675,200 elements with the octree mesh refinement instead of 16,777,216 i.e. (2 8 ) 3 elements for a regular voxel mesh with the same minimum element size.
None of the meshes shown in Fig. 5 are conformal and have multiple hanging nodes, i.e. nodes inside a mesh domain but not connected to all the surrounding elements. Multi-point constraints are used to enable the continuity of the mesh by constraining the hanging nodes to the nodes of larger neighbouring elements. Since the octree is balanced all the transitions between elements are 2:1 and therefore the constraints for hanging nodes which fall on edges of linear hexahedral elements (e.g. C3D8 in Abaqus) can be written as [19]: For the nodes which fall on faces the constraints are: where u hanging is the vector of displacement of the hanging node and u i master are the vectors of displacement of nodes of a master element. Constraints similar to (1) and (2) can be written for quadratic elements.

Mesh smoothing
After octree mesh refinement the refined meshes still exhibit one of the main disadvantages of voxel meshes as they have jagged surfaces which create artificial stress concentrations. The problem can be resolved by applying a smoothing algorithm to the surface nodes. Possible smoothing algorithms include Laplacian smoothing and Laplacian Eigen-decomposition [20]. This work uses Gaussian filtering method for smoothing of meshes as described e.g. by Taubin [21]: where 0 < λ < 1 is a weight and L p i ð Þ is a Laplacian operator whose simplest definition can be written as, essentially, a weighted distance between the node p i and N neighbouring nodes p j :  Table 1); c. Conformal mesh (Mesh 3 in Table 1); d. Octree-refined mesh (minimum level of refinement 4, maximum -7); e. Octree-refined mesh (minimum level of refinement 6, maximum -7); f.  The smoothing operator defined by Eq. (4) is the simplest definition as it assumes equal weights between the nodes. The operator can be applied iteratively with coefficient λ taking different values at each step.
Smoothing algorithms have been applied to smooth voxel meshes before but always starting with a regular voxel mesh which is often far from an accurate representation of the boundaries. This leads either to excessive distortions of elements and geometry or to the significant reduction of the geometry volume. Applying the smoothing algorithm to an octree refined mesh resolves the second problem as the elements on the surface are small and their reduction will not result in a significant change of volume. Excessive distortion of the elements is avoided by imposing a threshold on how much a node can move from its original position.

Model of an elliptic cylinder inclusion
An elliptic cylinder inclusion in a matrix has been selected as the simplest example to demonstrate the validity of the approach. The cross-section of the inclusion was an ellipse with major and minor axes equal to 0.3 mm and 0.2 mm, respectively. The matrix was modelled as a cube with a side equal to 1.0 mm. The materials of both the inclusion and matrix were isotropic with the Young's moduli equal to 200 and 100 GPa, respectively. The Poisson ratio for both materials was 0.3. The inclusion and the matrix are assumed to be perfectly bonded at the interface. The boundary conditions applied to the model were uniform displacements 0.1 mm and −0.1 mm applied at two opposite sides of the cube. Other sides of the cube were free of any loads or boundary conditions. The model was meshed using conformal (within Abaqus/CAE), voxel, octree and smoothed octree meshes. Illustrations of meshes are given in Fig. 6. Details of these types of meshes are given in Tables  1-3. In all cases, C3D8 elements were used for the analysis. Only one set of smoothing parameters has been used in this paper (λ = 0.3, 50 iterations).
Three conformal meshes were generated to give a reference solution which can be used as a benchmark for solutions obtained with other meshes. Convergence of the meshes was determined using the following criteriaconvergence of maximum von Mises stresses and normalised root mean square deviation (NRMSD) of local stresses along three paths shown in Fig. 7. The NRMSD was calculated as the root mean square deviation normalised by the mean value: where y i are the values of stresses in the reference model, b y i are the values of stresses in the model of interest and E Á ð Þ is mean value of a vector.
The maximum stresses and NRMSD along three selected paths in the model converged fast with the conformal mesh refinement as shown in Table 1. The solution obtained using the most refined conformal mesh consisting of approximately 3.1 × 10 6 elements was selected as a reference solution. The maximum stress error in a model, b σ max ; is written as a deviation from the value obtained with the reference mesh, σ max , i.e. b σ max À σ max ð Þ =σ max . The second most refined mesh   has the maximum stress error of 2.56% while other metrics were all below 1%. Two voxel meshes were compared against the selected reference mesh. The comparison is given in Table 2. The maximum stress error was found to be up to 30% but it should be pointed out that this error is the result of stress concentration at one or several elements. The metrics for average deviations decrease below 2% with increase of voxel mesh density, as expected, but still do not achieve the accuracy of a conformal mesh with a similar number of elements. While it is possible to increase the number of elements further, it becomes more computationally expensive to do so.
The mesh convergence of new meshing approaches was recorded using the same metrics as for conformal and voxel meshes. The octree-based meshes were refined by reducing the maximum element size (increasing the minimum level of refinement from 4 to 7) and reducing the minimum element size (increasing the maximum level of refinement from 7 to 8). The results of the comparison are given in Fig. 8 and Table 3. Reducing the maximum element size reduces the difference between the model and the reference solution but this effect is relatively small. In comparison, both the reduction of the minimum element size and smoothing have a strong effect on the difference between the model and the reference solution. The most refined model (minimum refinement -7, maximum refinement -8) has the NRMSD along three selected paths within 0.5% from the reference solution. However, the maximum stress error still remains as high as 8.7% even for the smoothed octree mesh. It is expected that further refinement would reduce the error but would be impractical because of the computational costs. Fig. 8 also shows total computational costs required to generate each of the meshes and obtain a corresponding solution.
The costs of obtaining a solution with the reference conformal mesh is also shown in Fig. 8. It is important to note that obtaining solutions with octree or smoothed octree meshes was faster than with the conformal mesh. Stress fields obtained with the generated meshes are shown in Fig. 9. A magnified view of the same stress fields in some of the meshes is shown in Fig. 10.
It should also be noted that calculation of NRMSD involves several interpolations which can introduce some errors. The first interpolation occurs within Abaqus/CAE when points for comparison are queried at the intersections of the path with the elements. The second interpolation is performed to sample the points from two different meshes to the same grid.

Textile composites
Section 3.1 showed the applicability of the octree refinement and smoothing techniques for the relatively simple case of an elliptical cylinder inclusion. However, the main motivation of the suggested approach is to apply it to various textile composites. For that purpose, unit cells of a plain weave, a sheared plain weave and an orthogonal 3D woven composite were generated using Tex-Gen. Geometry parameters of the unit cells are given in Tables 4 and 5 and the unit cells are shown in Fig. 11. To allow generation of acceptable conformal meshes all of the unit cells had an artificial gap between yarns.
Similar to the example with an elliptical inclusion, the unit cells of textile composites were meshed with conformal, voxel, octree and smoothed octree meshes. The conformal meshes were generated using TexGen and linear C3D4 elements were used for all of them. C3D8R elements were used for other meshes. Examples of the meshes are shown in Fig. 12. The properties of matrix and yarn materials are given in Table 6.
Convergence studies of conformal meshes were performed for each of the textiles. The convergence was judged using several metricsstress in a selected cross-section and NRMSD (Eq. (5)) along selected paths. The paths are shown schematically in Fig. 13. As for the example in Section 3.1, the converged solutions obtained with conformal meshes were selected as the reference solutions for each of the unit cells.
The following boundary conditions were applied to all unit cells: uniform displacements of 0.1 mm and −0.1 mm were applied to two opposite sides in the x-direction to simulate an extension in this direction, top surface of the unit cell was constrained in the out-ofplane direction.
Details of convergence studies of the conformal and voxel meshes are given in Appendix A. As in Section 3.1, the maximum stress error converges fast for the conformal meshes but NRMSDs for some paths usually remain above 2%, which can be explained by there being a more complex stress state in the unit cell of a textile composite than in the problem with a single inclusion. The voxel meshes converge more slowly and show up to 11% error in the maximum stress error and NRMSDs. These errors can be reduced by further refining the voxel mesh but it would result in large meshes which are prohibitively expensive to use in a FE analysis.
The solutions obtained with octree and smoothed octree meshes were compared to the reference solutions and detailed results are given in Tables 7-9 for plain weave, sheared plain weave and orthogonal 3D weave unit cells, accordingly. It was found that reducing the maximum element size (increasing the minimum level of refinement) improves accuracy of the solutions more than was observed in Section 3.1. Decreasing the minimum element size (increasing the maximum level of refinement) also has a significant effect on the accuracy. Smoothing also increases the accuracy of the solutions. In summary, deviation less than 7% from the reference solution can be achieved for all these examples. Maximum stress error and NRMSD along path '1' for all unit cells are shown in Fig. 14. The graphs also show the difference between the reference mesh and a conformal mesh with lower refinement. Distributions of von Mises stresses in longitudinal and transverse yarns of the plain weave composite for different levels of refinement are given in Fig. 15.
Total analysis cost for each of the performed simulations is also shown in Fig. 14. It can be seen that reducing the minimum element size, which results in larger meshes, results in growth of the analysis time. The cost of generating an octree mesh and performing analysis using this mesh is always lower than for an equivalent smoothing due to obvious extra operations required for the latter. The actual time to run the simulations in Abaqus were almost identical for both types of meshes as it mainly depends on the number of nodes and constraints rather than on a type of the mesh.

Discussion and conclusions
The meshing algorithm presented here can be used for automatic generation of FE meshes of unit cells of woven composites of arbitrary complexity. Combination of an octree voxel mesh and mesh smoothing algorithm helps to create meshes which have accuracy comparable to equivalent conformal meshes. The octree refinement procedure, which is used to create the mesh, approximates the internal geometry of the    the proposed method was benchmarked against conventional conformal meshes. It was found that the smoothed octree meshes can estimate local stresses along selected paths with an average deviation of less than 7% when compared to the converged conformal meshes. The error for the maximum stress was estimated to be within 10% for selected cross-sections. Obviously, there might be cases where the presented method does not result in low error e.g. in areas where geometrical features were not meshed with sufficient resolution or the smoothing algorithm did not eliminate stress concentrators. However, it is still expected that average stresses are within an acceptable range and can be achieved within reasonable time cost. Total analysis time, i.e. the time required to generate a mesh and to perform analysis using these meshes, was also compared. The time to generate octree and smoothed octree meshes increases almost exponentially with reduction of the minimum element size (increasing maximum level of refinement). However, for the meshes with lower Fig. 11. Unit cells of selected textiles: non-sheared plain weave, sheared plain weave, and orthogonal 3D woven (from left to right).  level of refinement (as e.g. equal to 7, which was reported here) the time to generate the mesh is often lower than the time required to generate a conformal mesh to achieve similar precision. For the meshes with a high level of refinement, the total analysis time is quite high because of the time required to generate the mesh and to pre-process extra constraints for the hanging nodes. This time will be improved by employing parallel execution of some parts of the code. It also should be mentioned that the precision of the solutions obtained with smoothed octree meshes can be reduced because of ill-shaped elements that can occur from the element smoothing at some boundaries. This can be improved by converting the octree mesh into a tetrahedral mesh with no hanging nodes [22,23] as shown in Fig. 16. These tetrahedral meshes can be post-processed after boundary smoothing to reshape the elements by moving internal nodes to improve element quality.
In summary, the presented smoothed octree meshes offer automatic meshing of any complex geometries for which conformal meshes are not always possible. The resulting mesh is not conformal but follows the boundary between materials with a selected precision. The octree refinement of the meshes reduces the computational cost which is advantageous in analysis of large unit cells. It was shown that the smoothed octree meshes can achieve precision comparable to conformal meshes.    Table 9 Comparison of solutions for the composite with orthogonal 3D weave reinforcement.