A novel computational framework for structural optimization with patched laminates

Abstract

Fiber patch placement (FPP) is a manufacturing technique for discrete variable stiffness composites. In the FPP approach, a structural component is assembled from a multitude of discrete fiber patches. However, due to the discontinuous fibers at patch edges, complex stress distributions occur. To date, a holistic FPP design framework that combines a tailored patch placement method with a dedicated mechanical model for the analysis of patched laminates does not exist. This article introduces a novel approach for the design of fiber patched laminates. It is based on the sequential placement of patches on a finite element shell mesh, using a critical element and angle selection routine in order to optimally locate and orientate fiber patches. They are added to the 3D mesh by employing a highly efficient kinematic draping algorithm. Strength-critical regions of the resulting fiber patched laminates are identified by state-of-the-art finite element analysis and extracted to a shear-lag–based mechanical submodel dedicated to the detailed analysis of patched laminates. The patch placement routine terminates once all design optimization criteria are met. The efficiency of applying optimized patch reinforcements on a continuous fiber-reinforced base laminate is demonstrated using the example of an individualized biomedical component. The work at hand presents the first patched laminate design framework combining a patch placement strategy coupled with a dedicated mechanical model. As a consequence, a substantial progress in the design of patch laminated structures is achieved.

Appendices

Appendix 1: Transformations

Principal load transformation matrix

$$ \boldsymbol{T} = \left[\begin{array}{llll} \cos(\theta)^{2} & \sin(\theta)^{2} & 2\sin(\theta)\cos(\theta) \\ \sin(\theta)^{2} & \cos(\theta)^{2} & -2\sin(\theta)\cos(\theta) \\ -\sin(\theta)\cos(\theta) & \sin(\theta)\cos(\theta) & \cos(\theta)^{2}-\sin(\theta)^{2} \end{array}\right] $$

(17)

1.1 3D vector transformation matrix

The 3D vector transformation matrix T_3D is given in (18). [N₁N₂N₃] are components of the element normal vector V_norm.

$$ \boldsymbol{T}_{3D} = \left[\begin{array}{lll} {N_{1}^{2}}(1-\cos(\theta)) + \cos(\theta) & N_{1}N_{2}(1-\cos(\theta)) - N_{3}\sin(\theta) & N_{1}N_{3}(1-\cos(\theta))+N_{2}\sin(\theta)\\ N_{1}N_{2}(1-\cos(\theta)) + N_{3}\sin(\theta) & {N_{2}^{2}}(1-\cos(\theta))+\cos(\theta) & N_{2}N_{3}(1-\cos(\theta))-N_{1}\sin(\theta) \\N_{1}N_{3}(1-\cos(\theta)) - N_{2}\sin(\theta) & N_{2}N_{3}(1-\cos(\theta))+N_{1}\sin(\theta) & {N_{3}^{2}}(1-\cos(\theta))+\cos(\theta) \end{array}\right] $$

(18)

Appendix 2: Draping on conformal maps

1.1 2.1 Determination of 3D submesh and subtriangulation

The conformal mapping tool (Choi and Lui 2015a) requires mesh preparation before it can be applied. The function only works on triangular meshes, and it requires the input mesh to be a manifold mesh with disk topology, that is, it can only have one continuous edge.

Hence, before calculating the subtriangulation, the finite element mesh has to be converted to a triangulation. That means splitting all quadrilateral shell elements that may exist in the FE mesh into triangles. The splitting is done by creating two triangles for each quadrilateral element. A mapping between the structural mesh and the triangulation is created, where each of the triangles is associated with its respective origin element of the finite element model.

The subtriangulation is found by starting at the triangle which contains the patch center location. In an iterative process, the neighboring triangles are added to the submesh, if their center is within an inclusion sphere of radius r_subtria centered at the initial triangle as can be seen in Fig. 13. The size of the inclusion sphere must be greater than the diagonal dimension of the patch to be draped and is herein chosen to

$$ r_{subtria} = \frac{3}{2} \sqrt{\left( \frac{l_{P}}{2}\right)^{2} + \left( \frac{w_{P}}{2}\right)^{2}} . $$

(19)

Fig. 13

Iterative subtriangulation determination. Initial triangle is shown in green, consecutive iterations are shown in light green and yellow. Triangles outside of the sphere of inclusion are shown in red

The check for the disk topology is done by starting at a vertex on the free edge, and moving along the edge, noting down all the vertices visited. Once back at the starting vertex, the number of vertices visited is compared with the total number of vertices on the edge. If the number of vertices visited is lower than the total number of free boundary vertices, it can be assumed that there are now more than one continuous edge on the subtriangulation, and thus disk topology requirement is violated. As a consequence, the subtriangulation iteration scheme is stopped and the final result is set to that from the previous loop. Figure 14 illustrates the approach.

Fig. 14

Disk topology check

1.2 2.2 Conformal mapping of subtriangulation to 2D unit disk

Once a valid 3D subtriangulation is found, the conformal mapping to a unit disk can be carried out using the method of Choi and Lui (2015b). The conformal mapping function (Choi and Lui 2015a) is based on numerical methods and therefore not completely accurate, especially close to the edges of the unit disk. Even though conformality distortions on the boundary are corrected in two steps with the aid of an iterative scheme using quasi-conformal theories, there can be some angle errors in the map. In order to minimize distortion of the draping process, two subtriangulation modification measures are implemented in the conformal mapping algorithm.

The first measure is to increase the size of the subtriangulation. The initial inclusion sphere (19) is chosen to be 50% larger than the patch diagonal length. This ensures that there is room to drape the patch. However, if, for some reason, the calculation of the conformal map of the initial subtriangulation fails, the draping function tries again, increasing the size of the radius r_subtria by 20% up to 10 times. This step is implemented to reduce the likelihood of unsuccessful draping simulations.

The second measure is the implementation of a mesh extension scheme. The accuracy of the mapping depends on how close the shape of the 3D subtriangulation is to that of a disk. Large deviations from disk shape can occur, if patch center location is close to an edge in the 3D shell mesh. To prevent excessive deviations of the subtriangulations from disk shape, the location of the patch center is checked on the flattened unit disk. If the patch center location α₀ is not within a radius of 0.1 around the center of the mapped unit disk, a mesh extension is carried out by adding artificial elements to the subtriangulation and the conformal mapping is computed again. The extension method increases the accuracy of the mapping in the area where the patch is draped, as it pushes the meshed area towards the center of the unit disc. The largest distortion at the edges is then limited to the extended mesh, which is not a part of the actual structural model. The subtriangulation extension method is illustrated in Fig. 15. It uses the elements on the mesh edge to calculate new nodes that are outside the subtriangulation. These new nodes are connected to the edge nodes, forming a set of new triangles.

Fig. 15

Subtriangulation extension scheme

1.3 2.3 Patch draping on 2D unit disk

Once the conformal mapping has been determined, patches can be draped on the 2D unit disk map, requiring the calculation of α- and β-points.

In order to find the α-points, a way to calculate distances in the two-dimensional map is needed. This approach is based on using barycentric coordinates to move between the three-dimensional subtriangulation and the two-dimensional map. The barycentric coordinate system of a triangle locates points using the vertices of the triangle. In Fig. 16a, the coordinates for a number of points on a triangle are shown. The interesting and useful property of this coordinate system is its ability to determine whether or not a point is within the boundaries of a triangle. If all the components of the coordinate are positive, the point lies within the triangle. This feature can be used to determine in which triangle a point on the two-dimensional map lies. This information can then be used to calculate the coordinates of the point in the three-dimensional subtriangulation. With this mechanism, the actual distances in three dimensions, between two points on the two-dimensional map can be calculated.

Fig. 16

Sampling of points for α-point calculation

Starting from the patch center location triangle, a geodesic line is drawn in the chosen draping direction up to the boundary of the unit disk. The line is sampled with a number of n_sample equidistant points. The location of these points is then calculated in the three-dimensional subtriangulation, using the barycentric coordinate method described above. The distance between each subsequent three-dimensional point is computed and a cumulative sum of the distance from the original point is taken. With infor- mation from the equidistant points in the two-dimensional map, and the cumulative distance in the three-dimensional map, a spline interpolation is used to construct a function that links distances on the 2D map with real distances on the 3D subtriangulation. This allows for the calculation of the corrected distances of the α-points on the two-dimensional map. An exemplary sampling is shown in Fig. 16b and c.

With the α-points calculated for all four α-lines, the quadrants between them have to be filled with β-points. To compute them, the two adjacent α-lines are used to fill in the points, starting from the point next to the center, then moving line by line outwards until the whole quadrant is filled. To calculate each of the β-points, an iterative scheme is used as described in Section 2.3.1. This iterative scheme requires a similar distance calculation as the determination of the α-points, because in the fishnet approach, distances between the points need to be fixed to d₁ = l_P/(2n_l) and d₂ = w_P/(2n_w), respectively. The distance calculation is implemented in the same way as for the α-point calculation.

To describe the calculation algorithm, a system for numbering the points is required. The four quadrants of the patch are defined by two α-lines. Each quadrant k is made up from n_l × n_w or n_w × n_l points, respectively. These points are numbered β_{i, j, k} where i, j gives the location on the quadrant with respect to the two α-points, and k gives the number of the quadrant.

The determination of β-points is illustrated in Fig. 17. Iterations start with the definition of points P₁ and P₂ at distances d₁ and d₂, respectively, from the previously determined points β_{i, j− 1,k} and β_{i− 1,j, k}. The required initial angles ω₁ and ω₂ are found from the orientation of the vectors α_2,j − α_2,j− 1 and α_1,i − α_1,i− 1, respectively. Next, the midpoint between points P₁ and P₂ is calculated. It is then used to compute updated angles ω₁ and ω₂. The new value for ω₁ is set to be the angle between β_{i, j− 1,k} and the midpoint $\beta _{i,j,k}^{(l)}$, with (l) denoting the current iteration number l. The updated value for ω₂ is set to be the angle between β_{i− 1,j, k} and the midpoint $\beta _{i,j,k}^{(l)}$. With the updated angles, the locations of points P₁ and P₂ can be updated as well and their distance is computed. The iteration scheme is executed until the distance between points P₁ and P₂ is below a threshold 𝜖_threshold. The final point β_{i, j, k} is then set to be the midpoint between P₁ and P₂.

Fig. 17

Iteration scheme for β-point calculation

1.4 2.4 Back-mapping to 3D mesh

Once the β-points have been calculated for all quadrants, the final step in the draping process is the translation of information from the α- and β-points to the 3D triangulation, and from there to the finite element mesh.

The first step is determining which triangles in the subtriangulation belong to the patch. To do this, the grid of α- and β-points in the two-dimensional map is used. There, the midpoints of all triangles are calculated, and for each quadrilateral grid box, the triangles falling within it are identified as a part of the patch. In addition, the quadrilateral grid boxes are used to assign fiber angles 𝜃 to the triangles, using the average direction of the two sides of the grid box, which are pointing in the fiber direction. An illustration for a grid box is shown in Fig. 18.

Fig. 18

Subtriangulation grid cell. Triangle selection and local fiber angle 𝜃 computation. Triangles associated with the grid box are shown in yellow, and red arrows are used to calculate the fiber angles for the triangles in the box

This selection method is applied to all the grid cells of the α- and β-point array, yielding results like the ones that can be seen in Fig. 19. By computing the orientation of every grid box, local fiber angle distributions within the patch and shearing distortions are considered.

Fig. 19

Draping triangle selection and angle assignment. The patch is shown in green, with the α- and β-grid points in blue, and the local fiber orientations in black

At this point, it is known which triangles of the subtriangulation represent the patch and at which angles the fibers lie within each triangle. This information is translated to the 3D finite element mesh by looping through all the triangles and checking which element they belong to. At the respective elements, the patch material, patch thickness, and local fiber orientation information is added to their elemental laminate definitions.

Appendix 3: Exoskeleton hip load cases

Table 2 gives the exoskeleton hip load cases.

Table 2 Loads F and M and constraints Φ^⋆ and $\bar {u}^{\star }$ on the exoskeleton hip component for a 50.1-kg-heavy 5th percentile female and an 85-kg-heavy male patient

Appendix 4: Material properties

Table 3 gives the properties of the employed materials.

Table 3 Material properties