Abstract

The composite sheet layup process involves stacking several layers of a viscoelastic prepreg sheet and curing the laminate to manufacture the component. Demands for automating functional tasks in the composite manufacturing processes have dramatically increased in the past decade. A simulation system representing a digital twin of the composite sheet can aid in the development of such an autonomous system for prepreg sheet layup. While finite element analysis (FEA) is a popular approach for simulating flexible materials, material properties need to be encoded to produce high-fidelity mechanical simulations. We present a methodology to predict material parameters of a thin-shell FEA model based on real-world observations of the deformations of the object. We utilize the model to develop a digital twin of a composite sheet. The method is tested on viscoelastic composite prepreg sheets and fabric materials such as cotton cloth, felt, and canvas. We discuss the implementation and development of a high-speed FEA simulator based on the VegaFEM library. By using our method to identify sheet material parameters, the sheet simulation system is able to predict sheet behavior within 5 cm of average error and have proven its capability for 10 fps real-time sheet simulation.

1 Introduction

Composites manufacturing methods form an integral part of the aerospace sector. Recently, the aerospace industry has gravitated towards an increase in the use of composites. The use of composites is projected to increase with a compound annual growth rate of about 10% year on year [1]. This growth in the utilization of composites has stirred a demand for automating functional tasks in the composite manufacturing processes. Out of the myriad of manufacturing methods, pre-impregnated with resin (prepreg) composite sheet layup method has been of key interest for the aerospace industry. Prepreg sheets are advantageous due to their controlled volume fraction, simple and inexpensive tooling, and ease of handling. These materials can be customized in terms of the types of fibers, weave type, and number of plies, allowing for diversification of component production and a broad range of applications [2,3]. Although the multi-component material systems delivered by prepreg sheets carry improved functionality, their manufacturing method is still susceptible to defects and quality issues. These defects can be classified as air gaps, wrinkles, bridging, etc. Currently, these layups are predominantly executed by skilled operators who place each ply (prepreg sheet) manually on the tooling and apply localized pressure using their hands or custom hand tools.

The high factor of human involvement in the layup processes introduces a potential variability and error in each subsequent layup. With an increase in demand for composite parts, it becomes crucial to increase production rates and improve quality to maintain process capability measures. This can be achieved through automation of the composite manufacturing methods. Presently, automated techniques such as automated fiber placement and automated tape layup are employed primarily for tooling which does not involve complex features. Automating a prepreg layup process introduces challenges that can only be addressed if the prepreg sheet behavior can be predicted accurately beforehand. For automated sheet layup, rectifying a defect in-process can be arduous as prepregs once draped are difficult to separate. Consequently, the separation can result in additional defects such as bridging, fiber misalignment, and impairment of previously laid up base sheets. As a result, executing a layup perfectly at once becomes a pivotal criterion for an automated prepreg layup process. This necessitates a need for a precise and low-latency simulation system that can accurately predict prepreg sheet behavior during the automated layup process. Such a simulation system acts as a digital twin of the prepreg and can be utilized to plan automation tasks.

In previous work accomplished at the Center for Advanced Manufacturing (CAM) at the University of Southern California (USC), Viterbi School of Engineering, the use of robots for performing automated prepreg sheet layup has been demonstrated [48]. Such a highly autonomous robotic cell needs estimation of the sheet parameters in situ for optimal and efficient planning. Along with accuracy, a high data transfer rate plays a critical role. In this research, the aim is to study the identification of different material parameters of a prepreg sheet, namely, bending stiffness, tensile stiffness, shear stiffness, damping stiffness, and surface density. These material parameters can help develop a FEM simulation model that can estimate and accurately simulate prepreg sheet behavior.

For developing the thin-sheet FEA model, we used the VegaFEM library [9]. In this paper, a modeling technique to accurately simulate prepreg sheets under fixed constraints is presented.1 An experimental methodology of recording sheet characteristics that can aid in the determination of key material parameters is presented as well. These recorded data are used to train a model that can evaluate the parameters with an acceptable accuracy using the Vega FEM library. The evaluated parameters are subsequently used to construct a digital twin represented by a force, damping, and mass matrix with the ability to emulate prepreg sheet behavior under external fixed constraints. The study then focuses on model evaluation and testing for different conditions. A detailed comparison of the parameter model predictions and experimental data is presented as well. We test our methodology on other materials such as cotton cloth, felt, and canvas. Additionally, we introduce a real-time sheet tracking system with depth image sensing to track the sheet for in-process validation.

The proposed in situ simulation system can aid in the development of an autonomous human–robot collaborative cell for prepreg sheet layup. In order to test the modeling concept, we developed a robotic cell setup and the sheet simulation as depicted in Fig. 1. Such a simulation system can particularly assist in the development of a grasp planning system that can aid in the handling of the prepreg sheets by robots. The potential applications of this simulation system are also discussed in detail.

Fig. 1
Left: The simulated prepreg under external forces and constraints. Right: The robotic cell with two Kuka iiwa R7 robots and one Kuka iiwa R14 robot.
Fig. 1
Left: The simulated prepreg under external forces and constraints. Right: The robotic cell with two Kuka iiwa R7 robots and one Kuka iiwa R14 robot.
Close modal

2 Related Work

Mechanical simulation of composite prepreg sheets is extensively used in the context of predicting the draping of prepreg over a mold. The sheet is represented by a mesh and fiber alignment, shear, and bending of the sheet can be predicted. A simple and computationally fast approach is to use kinematic simulations which only consider the mold geometry [11,12]. More accurate models use elasticity theory to compute the strains within the fabric [13]. Advanced finite element analysis (FEA)-based models are used to capture the deformation mechanics of the cloth [14]. In Ref. [15], a robot was used to place the flexible material on a doubly curved guiding tool. They used FEA-based models to predict the conformity. While FEA models are more accurate, they are often computationally expensive. In our work, we utilize a variant of FEA developed in the computer graphics community, providing a good trade-off between simulation speed, stability, and accuracy [16,17]. As a result, our method is not just reasonably accurate but also runs at interactive rates which permits us to perform more design iterations.

A survey work done in Ref. [18] presents an overview of all the techniques and models used to simulate fabrics. An overview of recent advancements in automated composite draping was presented in Ref. [19]. Progressive drape models which are a hybrid between kinematic and FEA simulations have also been proposed [20,21]. Such methods improve the accuracy compared to the purely kinematics model and also reduce the computational expense compared to FEA. Other commonly used models are particle systems [2225] which nonetheless suffer from an accuracy loss due to their inherent particle discretization.

The parameters of the model need to be tuned to match real-world thin-sheet observations. For cloth, this can be done using Kawabata plots [26]. Prepreg composite sheets are sturdier than cloth, however, necessitating precise displacement and force measuring equipment, rendering Kawabata plots less applicable. Instead, we develop an approach which uses FEA simulation and optical motion tracking to tune the material parameters for a given sheet, without explicitly measuring any internal elastic forces.

3 Simulation Model

We use thin-shell FEA to simulate the behavior of viscoelastic prepreg materials. Our thin-shell mathematical simulation model was described in computer graphics references [16,17], and we summarize it in this chapter for completeness. We note that these models were designed as a tradeoff between computational efficiency and accuracy, and have not been previously applied to real-world structures such as composite prepreg. We stress that accuracy alone is not sufficient for a successful material optimization system, as the speed by which design iterations can be explored also greatly affects one’s ability to obtain satisfactory results. Our system runs interactively at 10 fps, which enables rapid exploration of the material parameter optimization space. We incorporate layup domain knowledge, realistic constraint-handling, and real-time tracking of sheet deformation and optimize the material properties to match real-world prepreg observations.

The prepreg is represented as a triangulated mesh with x, y, z displacements of the vertices treated as simulation degrees-of-freedom (DOFs) of the model. The internal forces as a result of bending, shear, and stretching of the mesh govern the behavior of the material. Therefore, these internal forces need to be computed based on laws of elasticity for mechanical simulations. In our approach, we also compute the bending and tensile-shear force Jacobians to improve the computational efficiency of simulations.

3.1 Tensile and Shear Forces.

Our algorithm takes the triangle elements of the mesh as an input. We parameterize the 3D surface of the material using 2D parameter space represented by parameters u and v. Corresponding 3D weft and warp direction vectors U and V of the mesh can be computed in terms of these parameters. These vectors need not be orthonormal to each other after deformation. Consider the undeformed and deformed states of a triangle element in the mesh to illustrate the concept shown in Fig. 2. The 3D vertices of the triangle Pa, Pb, Pc are computed as a function of ua, va, ub, vb, and uc, vc.

Fig. 2
Undeformed (left) and deformed (right) states of a triangle in the mesh representing the sheet. The warp and weft vectors V and U are used to compute the tensile and shear strains.
Fig. 2
Undeformed (left) and deformed (right) states of a triangle in the mesh representing the sheet. The warp and weft vectors V and U are used to compute the tensile and shear strains.
Close modal
The weft and warp vectors are represented by weighted sums of the three parametric vertices of the triangle. We can formulate a linear system of six equations: iruiui=1, iruivi=0, irui=0, irviui=0, irvivi=1, and irvi=0, where i ∈ {a, b, c}. The weights rui and rvi can be precomputed using the equations: rua = d−1(vbvc), rva = d−1(ucub), rub = d−1(vcva), rvb = d−1(uauc), ruc = d−1(vavb), and rvc = d−1(ubua), where d = ua(vbvc) + ub(vcva) + uc(vavb). The system of six linear equations is solved to obtain the vectors U and V given by Eq. (1). Viscosity of the material is given by the evolution rates or rate change of these vectors given by Eq. (2). The vectors are then used to compute the Green–Lagrange strain tensor which consists of shear and tensile strains. Rate of change of these strains are then derived. Equations (3), (4), (5), and (6) give the representations.
U=i{a,b,c}ruiPi,V=i{a,b,c}rviPi
(1)
U=i{a,b,c}ruiPi,V=i{a,b,c}rviPi
(2)
ϵuu=12(UTU1),ϵuu=12(UTU)
(3)
ϵvv=12(VTV1),ϵuu=12(VTV)
(4)
ϵuv=12(UTVVTU)
(5)
ϵuv=12(UTV+VTU)
(6)
Deriving the weft, warp, and shear components of total elastic energy of the triangle with respect to vertex position gives us the force applied at the jth vertex of the triangle given by Eq. (7).
Fj=|d|2(σuu(rujU)+σvv(rvjV)+σuv(rujV+rvjU))
(7)
Stress tensor provides the values of the stresses σij used in Eq. (7). The relationship between stress and strain tensor is given by σ = + Eε′ where E and E′ are the elastic and viscosity stiffness matrices of the material and σ and ε are the 3D stress and strain vector. We then compute the force Jacobian which is necessary for implementation and efficiency of numerical techniques. The Jacobian for ith and jth vertex where both i, j ∈ {a, b, c} is computed using Eq. (8). If viscosity is also considered, then another contribution given by Eq. (9) has to be considered as well.
FjPi=|d|2(m,n{uu,vv,uv}σmϵn(ϵmTPiϵnPj)+m,n{uu,vv,uv}σm(PiϵmTPj))
(8)
FjPi=|d|2(m,n{uu,vv,uv}σmϵn(ϵmTPiϵnPj))
(9)

The stiffness component governs how stress–strain relationship affects the forces acting on the vertices of the triangle. The new position of the vertices is then found by taking into account internal and external forces during the simulation. The force acting due to bending stress is also superimposed with the tensile and shear forces to improve simulation accuracy. We will now discuss the method to compute the bending force.

3.2 Bending Forces.

The bending force is computed as a function of the hinge angle between two adjacent triangles in the mesh; we adopt the mathematical model of Ref. [16]. Figure 3 shows the two normals n1 and n2 corresponding to each triangle and to the most three neighboring triangles in the mesh. The angle between these normals or the bending angle is the hinge angle θ. Consider the total bending energy Eb=iψ(θi) as a function of θ summed over all possible hinges i of the mesh. The function ψ is an application specific function of the bending angle θ. We can obtain the bending force by differentiating the energy with respect to the vertex position x as F(x)=iψ and the hessian can be obtained as H(x)=iψHess(θi)+ψθiTθi, where Hess(θi) is the second-order derivative of θi with respect to x. In this work, the function ψ(θ) is given by the equation
ψ(θ)=k(2tan(θ2)2tan(θ¯2))2,θ(π,π)
(10)
where k is a constant dependent on material properties and θ¯ is the angle at the rest configuration. Details on how to compute the gradient and the hessian of the bending energy analytically are given in Ref. [16]. At each time-step, we update the bending forces and the Jacobian of the forces based on the vertex positions.
Fig. 3
(Left) Two adjacent triangles in the mesh and the bending angle θ between them. (Right) Three neighboring triangles for a triangle under consideration are shown.
Fig. 3
(Left) Two adjacent triangles in the mesh and the bending angle θ between them. (Right) Three neighboring triangles for a triangle under consideration are shown.
Close modal

4 Estimating Sheet Parameters

4.1 Overview.

In this section, we discuss the methodology of estimating material parameters for composite prepregs, cotton cloth, felt, and canvas. The experimentation was accomplished in three stages: data acquisition, data and parameter training/testing, and parameter optimization. Data acquisition was conducted within a physical environment through a guided manipulation of the various material sheets and collecting point cloud data of the position of each sheet as it moved. The training data optimization was conducted in a purely computational space, repetitively running the acquired data through an optimizing simulation sequence to characterize optimal parameter outputs. The outputs were then run through the simulation testing sequence to produce the final results. This can be seen in Fig. 4. This method was completed for all four material types separately. The parameter estimation approach will be discussed comprehensively in Sec. 4.3.

Fig. 4
Process overview: the initial state of the sheet is defined as the sheet configuration under initial boundary conditions. The releasing/released sheet state is defined as the sheet behavior after releasing one of the boundary conditions. After conducting physical experiments, initial mesh and observed data in the form of a mesh of the sheet are obtained from two sheet states, respectively. These data are further fed to the optimizer to acquire computed parameters.
Fig. 4
Process overview: the initial state of the sheet is defined as the sheet configuration under initial boundary conditions. The releasing/released sheet state is defined as the sheet behavior after releasing one of the boundary conditions. After conducting physical experiments, initial mesh and observed data in the form of a mesh of the sheet are obtained from two sheet states, respectively. These data are further fed to the optimizer to acquire computed parameters.
Close modal

4.2 Acquisition of Training and Testing Data.

The data collection was conducted using a Hexagon RS5 laser scanner and contact probe attached to the Romer absolute arm (87-Series). The class 2M laser scanner is hand operated and generates point cloud data. All scanning equipments are shown in Fig. 5. The sampling filter can be manually set to optimize the percentage of points recorded and exposure time based on the light in the sampling environment and color of the component being scanned to reduce noise and capture the target locations. The settings can be seen in Table 1.

Fig. 5
(a) Romer absolute arm, (b) laser scanner, (c) contact probe, and (d) contact probe variety
Fig. 5
(a) Romer absolute arm, (b) laser scanner, (c) contact probe, and (d) contact probe variety
Close modal
Table 1

Settings used on the Hexagon RS5 laser scanner

SettingsValue
Max capture rate752,000 points/s
Percentage25%
Exposure time200 μs
Point spacing0.052 mm
Line width130 mm
Sampling rate51 Hz
SettingsValue
Max capture rate752,000 points/s
Percentage25%
Exposure time200 μs
Point spacing0.052 mm
Line width130 mm
Sampling rate51 Hz

The contact probe was calibrated and operated using a TESA TKJ 3 mm Ruby Ball Probe. The probe and laser were subject to accuracy specifications, designated by the manufacturer in Table 2.

Table 2

Accuracy specifications for laser scanner and contact probe

EquipmentAccuracy
Laser scanner0.028 mm
Contact probe0.046 mm
EquipmentAccuracy
Laser scanner0.028 mm
Contact probe0.046 mm

4.2.1 Sheet Preparation.

To prepare each sheet for scanning, quarter-inch markers were placed at even intervals along the sheet in a 17 × 17 grid, for a total of 289 markers. The markers used were white 3M double-sided foam tape squares, cut to the correct size. The markers were raised from the surface, allowing for easier detection by the scanner. One of the sample materials used was white in color and the markers had to be colored in black to be picked up by the scanner, otherwise the color difference was considered optimal for scanning purposes with the settings previously mentioned. The four locations of the clamps obscured the markers and were recorded separately via the contact probe for use within the simulation as a fixed point (Fig. 6).

Fig. 6
Sheets were prepared by attaching 289 quarter-inch markers to each sheet
Fig. 6
Sheets were prepared by attaching 289 quarter-inch markers to each sheet
Close modal

4.2.2 Observed Data Generation.

These data were collected during two trials for each material tested. Each edge of a sheet was grasped using four vertically fixed clamps, allowing the sheet to rest suspended between them, as shown in Fig. 7. The clamping locations were chosen at differing and unique distances from adjacent corners. The collection period for each trial featured a five-stage setup, with each stage introducing a new modification to sheet positioning by moving the clamp with respect to the room, while retaining the clamping position on the sheet. The translational movement of the clamp was chosen arbitrarily, but all movements had a change in distance of no more than 350 mm.

Fig. 7
Clamps fixed the sheet configuration at four points. This is defined as the initial state.
Fig. 7
Clamps fixed the sheet configuration at four points. This is defined as the initial state.
Close modal

The first stage required placing the sheet under the initial boundary conditions of the first stage, resulting in the sheet suspended in a relaxed, horizontal position. The next three stages involved isolated movement of only one grasping location to a new position, with a general increase in the z-axis, coinciding with x- and y-axis movement toward the center of the sheet. The last stage for every trial was a release of the clamp, allowing the material to settle into a hanging position. Within each stage, each clamp changed location only once and was allowed to settle into position, until no visible movement could be detected. At the beginning of each stage, after clamp relocation had occurred, the new clamp locations were collected via the contact probe and the sheet was scanned via the laser scanner. This procedure can be observed in Fig. 8.

Fig. 8
Labeled sequence of stages depicting sheet position movements during one trial
Fig. 8
Labeled sequence of stages depicting sheet position movements during one trial
Close modal

The point cloud data collected from each stage were processed using pc-dmis and exported as an XYZ file for post-processing. From each stage, each sheet contained thousands of data points, with each marker averaging 300 points. The simulator requires only one point from each of the markers, so each XYZ file was processed through blender software using an original python script to isolate the center point of each marker and export them as 3D points, relative to the world coordinate system of the scanner base. A comparative image of the same sheet is shown in Fig. 9, demonstrating the reduction from the point cloud marker clusters to single points.

Fig. 9
Point cloud clusters (left) versus single point vertices (right) on one sheet
Fig. 9
Point cloud clusters (left) versus single point vertices (right) on one sheet
Close modal

On successful construction of the initial mesh, the boundary conditions of the sheet are changed by releasing one of the clamps. Note that the remaining clamps should not be moved to maintain consistency in the entire process. The sheet state after changing the boundary conditions is defined as released state. Figure 8 shows the difference between the initial state and the released state.

Since no mesh is required to generate from the released state, the point cloud data of such a state are clustered to represent each marker. The data set acquired in this process is further defined as the observed data. Figure 9 shows the point cloud data clustering process.

4.3 Model Parameter Estimation

4.3.1 Sheet Simulation System.

The crucial element of the parameter estimation process is the composite sheet simulation system. Figure 10 gives the block diagram for the proposed simulation system. The simulator system utilizes the initial mesh as geometric information input and applies the model parameters to construct the composite sheet model.

Fig. 10
Simulation process overview: the simulator used the material parameters and mesh information to predict sheet behavior under specified conditions
Fig. 10
Simulation process overview: the simulator used the material parameters and mesh information to predict sheet behavior under specified conditions
Close modal

The model parameters are categorized into two types: (1) material parameters and (2) integrator parameters. Material parameters consists of the surface density and the internal force parameters, which are tensile stiffness, shear stiffness, and bending stiffness. On the other hand, integrator parameters include damping stiffness and damping mass.

Once the sheet model is constructed [27], the numerical integrator applies the boundary conditions and external forces, such as gravity, to the sheet model and solves the deformation equation [28].

After the predicted mesh is generated, the prediction error is obtained by comparing the predicted mesh and the observed data. The prediction error, E, is a function of model parameters, P, initial mesh, M, and observed data, O. Algorithm 1 is used to calculate the prediction error function.

Algorithm 1

1: Let dO

2: MpreSimulate(P,M)

3: for allddo

4:  vFindClosest(d,Mpre)

5:  disp|vd|

6: end for

7: ErrorAverage(disp)+0.5*Maximum(disp)

Let d be a data point in the observed data, O. After getting the predicted mesh, Mpre, by applying P and M to the simulator, we query each d to find the closest vertex, v, on Mpre. The prediction error is then defined as the average displacement between v and d plus half of the maximum displacement across the entire mesh. The prediction error is used for parameter optimization, which will be discussed in Sec. 4.3.3.

4.3.2 Parameter Boundary Selection.

The model parameter identification uses a nonlinear optimizer to compute the optimal parameters for simulating the composite prepreg. However, the optimizer is not required to compute all model parameters. Some of these parameters can be measured directly.

As mentioned in Sec. 4.3.1, model parameters comprise of material parameters and integrator parameters. Since integrator parameters are not related to model construction, we set damping stiffness and damping mass to 1.0 and 0.0, respectively. The surface density can be measured directly by scaling the sheet and dividing it by the surface area. Therefore, the remaining parameters, tensile stiffness, shear stiffness, and bending stiffness, are the parameters that require optimization.

To ensure effective performance of the optimizer, appropriate upper bound, lower bound, and initial values for the parameters are required. Figure 11 shows the input–output diagram for parameter boundary selection. Tensile stiffness and shear stiffness are sampled into three categories: (1) 500, (2) 5000, and (3) 50,000. For bending stiffness, the parameter is sampled into (1) 0.01, (2) 0.001, and (3) 0.0001.

Fig. 11
Input–output diagram for parameter boundary selection
Fig. 11
Input–output diagram for parameter boundary selection
Close modal

After getting three samplers for each parameter, we shuffle them and get 27 candidates to test for parameter feasibility. All candidates that cause system failure are eliminated, and those who have the smallest prediction error, or are comparable to the smallest one, are highlighted. Candidate {5000, 5000, 0.001} has the best overall performance among 27 candidates, and therefore, it is selected to be the initial values set for the parameter optimization. Then, the initial values are used as the median for the parameter boundary. Thus, the upper boundary for the {tensile stiffness, shear stiffness, bending stiffness} is {9000, 9000, 0.0011}. The lower boundary for the parameter set is {1000, 1000, 0.0009}.

4.3.3 Optimization Algorithm.

The optimization library used in this work is NLOPT [29], an open-source library for nonlinear optimization. The selected algorithm is improved stochastic ranking evolution, a global-gradient-free optimization algorithm [30]. Figure 12 shows the block diagram for the optimization process.

Fig. 12
Parameter optimization process
Fig. 12
Parameter optimization process
Close modal
The parameter optimizer utilizes the training sets and initial parameters as inputs and calculates the prediction errors for each sampler. The goal for the optimizer is to find the parameters that minimize the sum of the errors from the training set.
{m1,,m6}M{o1,,o6}OPoptargminPi=1k(Ei(P,mi,oi))
(11)

The optimization problem can be expressed as Eq. (11). Let M be initial meshes in the training set, which contains mesh data, mi. O is the observed data set, which contains observed data oi. Recall that the prediction error, Ei, is defined in Algorithm 1. Since five samples are used for model training, k is set to 5. The optimizer tries to find the parameters that reduce both the average displacement difference and the maximum difference between the predicted mesh and the observed data for each training set. This error value was set to (average error + 2* Max error). Once the error converges, the identified parameters, Popt, are used to simulate the composite sheet in real-time.

5 Results

5.1 Experimental Specifics.

Four materials are considered in this work: prepreg composite sheets, common cotton cloth, felt, and canvas. The composite sheet was supplied by Boeing Inc. and came as 3 ft × 4 ft sheets, while all other sheets were purchased locally and were 2 ft × 2 ft in dimension, shown in Fig. 13. As discussed previously, the prepreg composite sheet contains viscoelastic properties due to the resin contents, while all other sheets contain no resin or viscous content. The elastic materials all have a standard/uniform weave. Details of the elastic materials are shown in Table 3.

Fig. 13
Sheets used for four material samples. Left: elastic fabric materials. Right: viscoelastic prepreg material.
Fig. 13
Sheets used for four material samples. Left: elastic fabric materials. Right: viscoelastic prepreg material.
Close modal
Table 3

Elastic material measurements

MeasurementClothCanvasFelt
Side1611 mm613 mm615 mm
Side 2611 mm613 mm612 mm
Thickness0.1 mm0.65 mm1.8 mm
Mass42.1 g168.9 g66.9 g
Surface density0.1126 kg/m20.4509 kg/m20.0018 kg/m2
MeasurementClothCanvasFelt
Side1611 mm613 mm615 mm
Side 2611 mm613 mm612 mm
Thickness0.1 mm0.65 mm1.8 mm
Mass42.1 g168.9 g66.9 g
Surface density0.1126 kg/m20.4509 kg/m20.0018 kg/m2

The composite prepreg material provided by Boeing came with manufacturer density and thickness specifications (refer Table 4). The elastic materials came with no specifications and thus required density calculations and thickness measurements, shown in Table 3.

Table 4

Viscoelastic material specifications

MeasurementComposite sheet
Side 11.17 m
Side 20.975 m
Thickness0.3 mm
Poisson’s ratio0.3
Surface density0.3 kg/m2
MeasurementComposite sheet
Side 11.17 m
Side 20.975 m
Thickness0.3 mm
Poisson’s ratio0.3
Surface density0.3 kg/m2

The physical experiment process described previously was repeated twice for each material type for a total of eight trials. In the case of the composite sheet, two different sheets were used, as it was determined that the combination of experimental movements and exposure to air may degrade the sheet, possibly providing poor results. The elastic materials all used the same sheet twice. Of the two trials for each material, the data collected from one trial were used for training purposes, while the other was used for testing purposes.

5.2 Sheet Parameter Estimation.

Estimation of sheet parameters was accomplished through the simulator, as previously discussed. The internal parameters were considered to be in one state throughout the simulation procedure, but research has shown that linear elastic woven fabrics show a variance in many of their internal parameters due to the anisotropic behaviors of the fabric [31]. As such, manual manipulation of the parameters from the initial internal parameters was needed to find an ideal starting range for the optimization. The initial parameters were found experimentally and were then run through the optimizer, providing lowest error results.

Tables 5 and 6 highlight the training and testing performance of the simulator in terms of the average error and maximum error for the training and testing data sets. Table 5 gives the results for the training of fabric materials and Table 6 gives the results for the testing of fabric material.

Table 5

Fabric material training data

Fabric material training
MaterialAvg error (cm)Maximum error (cm)
Felt2.37.8
Cloth1.15.5
Canvas1.46.0
Fabric material training
MaterialAvg error (cm)Maximum error (cm)
Felt2.37.8
Cloth1.15.5
Canvas1.46.0
Table 6

Material testing data

Fabric material testing
MaterialAvg error (cm)Maximum error (cm)
Felt2.68.2
Cloth1.65.6
Canvas1.56.2
Fabric material testing
MaterialAvg error (cm)Maximum error (cm)
Felt2.68.2
Cloth1.65.6
Canvas1.56.2

Figures 14, 15, and 16 show the configuration of stage 5 for the felt, cloth, and canvas sheets, respectively. In each figure, the upper image is the observed position photograph, the middle image is the mesh generated from the scanned point cloud data, and the lower image is an image of the simulation mesh at the end of the stage.

Fig. 14
Trial 1 stage 5 results of felt sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Fig. 14
Trial 1 stage 5 results of felt sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Close modal
Fig. 15
Trial 1 stage 5 results of cloth sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Fig. 15
Trial 1 stage 5 results of cloth sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Close modal
Fig. 16
Trial 1 stage 5 results of canvas sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Fig. 16
Trial 1 stage 5 results of canvas sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Close modal

Figures 17 and 18 depict the initial and final positions for the composite sheet trials respectively. Within the initial stage shown in Fig. 17, consistency can be seen between the photo, observed data, and simulation data.

Fig. 17
Trial 1 stage 1 results of composite sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Fig. 17
Trial 1 stage 1 results of composite sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Close modal

Within the final stage, there was consistent error across all material types, however this error was higher than all other stages indicating that this position was the most difficult for the simulator to emulate. Figure 18 depicts this position and displays results from the observed data and simulated data.

Fig. 18
Trial 1 stage 5 results of composite sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Fig. 18
Trial 1 stage 5 results of composite sheet: (upper) observed position image, (middle) generated mesh, and (lower) simulated mesh
Close modal

The error from the training data is shown in Table 7. The optimal parameters were a bending stiffness of 4.77682 × 107 N/m2 and a shear stiffness of 2942.93 N/m. The optimization process began with initial values based on the data sheet supplied by Boeing in Table 8.

Table 7

Composite training and testing data

Composite sheet material
SheetAvg error (cm)Maximum error (cm)
Training 11.610.2
Testing 11.413.9
Testing 22.213.6
Composite sheet material
SheetAvg error (cm)Maximum error (cm)
Training 11.610.2
Testing 11.413.9
Testing 22.213.6
Table 8

Initial optimizer parameters

ParameterInitial value
Sheet thickness0.0003 m
Poisson’s ratio0.3
Sheet density0.3 kg/m2
Shear stiffness3.3 × 102 N/m
Bending stiffness1.1 × 108 N/m2
ParameterInitial value
Sheet thickness0.0003 m
Poisson’s ratio0.3
Sheet density0.3 kg/m2
Shear stiffness3.3 × 102 N/m
Bending stiffness1.1 × 108 N/m2

5.3 Sheet Simulation.

The results obtained from the material parameter estimation are then utilized to develop a simulation system that can predict the behavior of the material under varied constraints. This section focuses on one of the potential applications of such a high-fidelity simulator. In composite prepreg sheet layups, as discussed earlier, it is important for the sheet to be held at appropriate locations such that it does not result in any potential defects during the layup. A sheet simulation system that can predict the appropriate locations for grasping the sheet can be instrumental in process planning. The key elements of such a simulation system would be (1) a sheet simulation model generated using the estimated material parameters and (2) a real-time sheet tracking system that can generate a mesh of the sheet in real-time at high frame rates, and that can then be used for comparison with simulated result. To demonstrate the feasibility of such a system, an experiment was designed using the carbon fiber reinforced epoxy sheet provided by Boeing Inc.

The methodology for material parameter estimation proposed earlier was used to generate an appropriate model for the Boeing composite sheet. The simulation system was built using the VegaFEM library [9]. Figure 19 shows the block diagram of the proposed real-time sheet tracking and simulation process.

Fig. 19
Real-time sheet simulation process
Fig. 19
Real-time sheet simulation process
Close modal

The dimensions of the composite sheet used in this experiment were 4 ft × 3 ft, which required a multi-camera system to track the entire sheet. The real-time sheet tracking system proposed in this study consists of three Realsense D415 sensors. The entire sheet is captured by fusing the red-green-blue-D feed from all the three cameras. Color-based filtering techniques are employed to achieve filtering of the prepreg sheet from rest of the scene. Re-sampling is performed on resulting points to obtain a uniform distribution of the filtered points. The face normals are then recomputed using these points. Surface reconstruction is performed in scale space by implementing advancing front surface reconstruction [32]. A scale space describes the point set at a dynamic scale, and this additional dimension allows us to control the degree of smoothness required for the reconstruction [33]. After post-processing, a surface mesh with around 6000 triangles was obtained.

To test the system performance, we constrained the composite sheet in four different states. Initially, the sheet is constrained at four grasping locations. A human and two 7-DOF robotic manipulators KUKA iiwa R7 are used to apply the fixed constraints. The real-time sheet tracking system captures the mesh and grasping locations for the corresponding state. In the next step, one of the grasping locations is released. The sheet is allowed to settle to its minimum energy state and the real-time sheet tracking system captures the sheet in the respective state. This exercise is repeated for three other states by changing the grasping locations and capturing the sheet behavior. These data are then compared to the predictions generated by the simulation system. Figure 20 contrasts the four states in which the composite sheet was suspended and the corresponding simulator predictions.

Fig. 20
Real-time sheet simulation result: the first row shows the predicted mesh from the simulator. The second row shows the actual behavior of the composite sheet.
Fig. 20
Real-time sheet simulation result: the first row shows the predicted mesh from the simulator. The second row shows the actual behavior of the composite sheet.
Close modal

The proposed simulation and real-time sheet tracking system can play a pivotal role in predicting optimal grasping locations for the draping of composite sheet. These grasping locations can then be used to deploy a human robot collaborative cell where the robots can aid the human in the draping process by holding the sheet appropriately [34]. Furthermore, the sheet tracking data can act as a rectifying input for the grasping location in case of sub-optimal simulator predictions. Such a conjunctive system can aid in streamlining the prepreg layup process and achieve an overall higher degree of automation.

6 Conclusion

In this work, we presented a methodology to predict material parameters to generate an energy-based model of various materials. We showcased how these models aid in simulating the material behavior. In addition, we discussed the implementation and development of a high-speed thin-shell simulator based on the VegaFEM library. The procedure followed for data collection for model training and testing is highlighted as well. The inputs to the parameter model are elaborated along with the boundary conditions pertaining to fixed constraints and external forces. Validation of the developed model for various test cases is presented and discussed in depth.

This work has provided detailed insights into the simulation of multi-material composite material-based systems, particularly prepreg sheets along with other non-viscoelastic materials. Emphasis has been laid on predicting behavior of prepreg composite sheets under different types of constraints and on simulating how the sheet will behave under these conditions. The prediction capabilities of our system open new avenues in the domain of automated composite layup processes, particularly in planning for the development of completely autonomous and human–robot collaborative cells. This system can be employed to plan and execute a defect-free layup using robotic assistants, using the simulator to generate optimal locations to grasp the prepreg sheets during a layup process. In the future, we aim to explore the capability of employing this simulator in trajectory and motion planning of collaborative robots mainly used for manipulation of prepreg sheets, thus enabling the automation of composite layup tasks.

Footnote

1

This paper is based on the work presented by the authors at the 2021 ASME Manufacturing Science and Engineering Conference [10].

Acknowledgment

This work is supported in part by National Science Foundation Grant No. 1925084. Opinions expressed are those of the authors and do not necessarily reflect opinions of the sponsors.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The authors attest that all data for this study are included in the paper. Data provided by a third party are listed in Acknowledgment.

References

1.
Mazumdar
,
S.
,
Pichler
,
D.
,
Benevento
,
M.
,
Seneviratine
,
W.
,
Liang
,
R.
, and
Witten
,
E.
,
2020
, “
American Composites Manufacturing Association 2020 State of the Industry Report
,” http://compositesmanufacturingmagazine.com/2020/01/2020-state-of-the-industry-report/
2.
Fleischer
,
J.
,
Teti
,
R.
,
Lanza
,
G.
,
Mativenga
,
P.
,
Möhring
,
H.-C.
, and
Caggiano
,
A.
,
2018
, “
Composite Materials Parts Manufacturing
,”
CIRP Annals
,
67
(
2
), pp.
603
626
.
3.
Barile
,
C.
,
Casavola
,
C.
, and
De Cillis
,
F.
,
2019
, “
Mechanical Comparison of New Composite Materials for Aerospace Applications
,”
Comp. Part B: Eng.
,
162
, pp.
122
128
.
4.
Malhan
,
R. K.
,
Kabir
,
A. M.
,
Shah
,
B.
,
Centea
,
T.
, and
Gupta
,
S. K.
,
2018
, “
Automated Prepreg Sheet Placement Using Collaborative Robotics
,”
North America Society for the Advancement of Material and Process Engineering (SAMPE)
,
Long Beach, CA
,
May 21–24
.
5.
Malhan
,
R. K.
,
Kabir
,
A. M.
,
Shah
,
B.
,
Centea
,
T.
, and
Gupta
,
S. K.
,
2018
, “
Hybrid Cells for Multi-Layer Prepreg Composite Sheet Layup
,”
IEEE International Conference on Automation Science and Engineering (CASE)
,
Munich, Germany
.
6.
Malhan
,
R. K.
,
Kabir
,
A. M.
,
Shah
,
B.
,
Centea
,
T.
, and
Gupta
,
S. K.
,
2019
, “
Determining Feasible Robot Placements in Robotic Cells for Composite Prepreg Sheet Layup
,”
ASMEs 14th Manufacturing Science and Engineering Conference
,
Erie, PA
,
June 10–14
.
7.
Malhan
,
R. K.
,
Kabir
,
A. M.
,
Shah
,
B. C.
, and
Gupta
,
S. K.
,
2019
, “
Identifying Feasible Workpiece Placement With Respect to Redundant Manipulator for Complex Manufacturing Tasks
,”
IEEE International Conference on Robotics and Automation (ICRA)
,
Montreal, Canada
,
May 20–24
.
8.
Malhan
,
R. K.
,
Shembekar
,
A. V.
,
Kabir
,
A. M.
,
Bhatt
,
P. M.
,
Shah
,
B.
,
Zanio
,
S.
,
Nutt
,
S.
, and
Gupta
,
S. K.
,
2021
, “
Automated Planning for Robotic Layup of Composite Prepreg
,”
Robot. Comput. Integr. Manuf.
,
67
, p.
102020
.
9.
Sin
,
F. S.
,
Schroeder
,
D.
, and
Barbič
,
J.
,
2012
, “
Vega FEM Library
,” http://www.jernejbarbic.com/vega
10.
Chen
,
Y.
,
Joseph
,
R.
,
Kanyuck
,
A.
,
Khan
,
S.
,
Malhan
,
R. K.
,
Manyar
,
O. M.
,
McNulty
,
Z.
,
Wang
,
B.
,
Barbic
,
J.
, and
Gupta
,
S. K.
,
2021
, “
A Digital Twin for Automated Layup of Prepreg Composite Sheets
,”
Manufacturing Sciences and Engineering Conference
,
Virtual
,
June
.
11.
Aono
,
M.
,
Breen
,
D. E.
, and
Wozny
,
M. J.
,
1994
, “
Fitting a Woven-Cloth Model to a Curved Surface: Mapping Algorithms
,”
Comput. Aided Des.
,
26
(
4
), pp.
278
292
.
12.
Wang
,
J.
,
Paton
,
R.
, and
Page
,
J.
,
1999
, “
The Draping of Woven Fabric Preforms and Prepregs for Production of Polymer Composite Components
,”
Compos. Part A: Appl. Sci. Manuf.
,
30
(
6
), pp.
757
765
.
13.
Aono
,
M.
,
1990
, “
A Wrinkle Propagation Model for Cloth
,”
CG International’90
,
Springer
, pp.
95
115
.
14.
Collier
,
J. R.
,
Collier
,
B. J.
,
O’Toole
,
G.
, and
Sargand
,
S.
,
1991
, “
Drape Prediction by Means of Finite-Element Analysis
,”
J. Text. Inst.
,
82
(
1
), pp.
96
107
.
15.
Krogh
,
C.
,
Glud
,
J. A.
, and
Jakobsen
,
J.
,
2019
, “
Modeling the Robotic Manipulation of Woven Carbon Fiber Prepreg Plies Onto Double Curved Molds: A Path-Dependent Problem
,”
J. Compos. Mater.
,
53
(
15
), pp.
2149
2164
.
16.
Tamstorf
,
R.
, and
Grinspun
,
E.
,
2013
, “
Discrete Bending Forces and Their Jacobians
,”
Graph. Models
,
75
(
6
), pp.
362
370
.
17.
Volino
,
P.
,
Magnenat-Thalmann
,
N.
, and
Faure
,
F.
,
2009
, “
A Simple Approach to Nonlinear Tensile Stiffness for Accurate Cloth Simulation
,”
ACM Transactions on Graphics
,
28
(
4
), pp.
1
16
.
18.
Ng
,
H. N.
, and
Grimsdale
,
R. L.
,
1996
, “
Computer Graphics Techniques for Modeling Cloth
,”
IEEE Comput. Graph. Appl.
,
16
(
5
), pp.
28
41
.
19.
Elkington
,
M.
,
Ward
,
C.
, and
Sarkytbayev
,
A.
,
2017
, “
Automated Composite Draping: A Review
,”
SAMPE SEATTLE 2017
,
Seattle, WA
,
May 22
.
20.
Hancock
,
S. G.
, and
Potter
,
K. D.
,
2006
, “
The Use of Kinematic Drape Modelling to Inform the Hand Lay-Up of Complex Composite Components Using Woven Reinforcements
,”
Compos. Part A: Appl. Sci. Manuf.
,
37
(
3
), pp.
413
422
.
21.
Sharma
,
S. B.
, and
Sutcliffe
,
M. P. F.
,
2003
, “
Draping of Woven Fabrics: Progressive Drape Model
,”
Plast. Rubber Compos.
,
32
(
2
), pp.
57
64
.
22.
Breen
,
D. E.
,
House
,
D. H.
, and
Getto
,
P. H.
,
1992
, “
A Physically-Based Particle Model of Woven Cloth
,”
Visual Comput.
,
8
(
5–6
), pp.
264
277
.
23.
Breen
,
D. E.
,
House
,
D. H.
, and
Wozny
,
M. J.
,
1994
, “
A Particle-Based Model for Simulating the Draping Behavior of Woven Cloth
,”
Text. Res. J.
,
64
(
11
), pp.
663
685
.
24.
House
,
D.
, and
Breen
,
D.
,
2000
,
Cloth Modeling and Animation
,
CRC Press
,
Boca Raton, FL
.
25.
Provot
,
X.
,
1995
, “
Deformation Constraints in a Mass-Spring Model to Describe Rigid Cloth Behaviour
,”
Graphics Interface '95
,
23
(
19
), pp.
147
147
.
26.
Kawabata
,
S.
,
1980
,
The Standardization and Analysis of Hand Evaluation
, 2nd,
Textile Machinery Society of Japan
,
Osaka, Japan
.
27.
Baraff
,
D.
, and
Witkin
,
A. P.
,
1998
, “
Large Steps in Cloth Simulation
,”
SIGGRAPH '98: Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques
, pp.
43
54
.
28.
Sin
,
F. S.
,
Schroeder
,
D.
, and
Barbič
,
J.
,
2013
, “
Vega: Nonlinear FEM Deformable Object Simulator
,”
Comput. Graph.
,
32
(
1
), pp.
38
50
.
29.
Johnson
,
S. G.
, “
The NLopt Nonlinear-Optimization Package
,” http://github.com/stevengj/nlopt
30.
Runarsson
,
T. P.
, and
Yao
,
X.
,
2005
, “
Search Biases in Constrained Evolutionary Optimization
,”
IEEE Trans. Syst. Man Cybernet. Part C: Appl. Rev.
,
35
(
2
), pp.
233
243
.
31.
Volino
,
P.
,
Magnenat-Thalmann
,
N.
, and
Faure
,
F.
,
2009
, “
A Simple Approach to Nonlinear Tensile Stiffness for Accurate Cloth Simulation
,”
ACM Trans. Graph.
,
28
(
4
), pp.
1
16
.
32.
Da
,
F.
, and
Cohen-Steiner
,
D.
,
2020
,
Advancing Front Surface Reconstruction
, 5.0.2 ed.,
CGAL Editorial Board
.
33.
Digne
,
J.
,
Morel
,
J. M.
,
Souzani
,
C. M.
, and
Lartigue
,
C.
,
2011
, “
Scale Space Meshing of Raw Data Point Sets
,”
Comput. Graph. Forum
,
30
(
6
), pp.
1630
1642
.
34.
Manyar
,
O. M.
,
Desai
,
J. A.
,
Deogaonkar
,
N.
,
Joseph
,
R. J.
,
Malhan
,
R. K.
,
McNulty
,
Z.
,
Wang
,
B.
,
Barbič
,
J.
, and
Gupta
,
S. K.
,
2021
, “
A Simulation-Based Grasp Planner for Enabling Robotic Grasping During Composite Sheet Layup
,”
IEEE International Conference on Robotics and Automation (ICRA)
,
Xi’an, China
,
May 30–June 5
.