Impact Statement
Engineers frequently use Finite Element models to predict the performance of complex structural systems.Distortions of the structure’s geometry from fabrication and manufacturing processes play a significant role in determining ultimate strength and structural stability of the system. Existing methods for quantifying distortions heavily rely on assumptions that limit their practicality. 3D point cloud data facilitates capturing the actual geometry of a distorted structure, potentially overcoming these limitations. This work presents a process for incorporating point cloud measurements of distortions into FE models though a Gaussian process interpolation.We demonstrate how this approach can yield strength estimates that are comparable to established methods while being applicable to a wider range of applications and quantifying a performance envelope for risk analysis.
1. Introduction
1.1. Background
The safe operation of structures requires inspection and evaluation to establish the actual state of the structure. Inspection results can then be used to conduct calculations to ensure that load capacity is sufficiently greater than load demand. Inspections typically occur either when the structure is degraded by instantaneous damage such as blast or impact, or on a regular basis to monitor gradual damage from fatigue or environmental exposure. With the growing availability and capability of advanced monitoring and inspection techniques that provide full-field 3D data, there is increasing interest in using quantitative 3D inspection data to update finite element (FE) models to reflect the actual state of the structure in the form of a digital twin. This updated FE model can provide meaningful information about current serviceability and remaining service life that an idealized model cannot provide.
The specific interest here is the use of point cloud data collected during ship structural inspections to update FE models to introduce initial imperfections for compressive ultimate strength computations, motivated by prior efforts in large-scale civil structures, for example, Zhao et al. (Reference Zhao, Tootkaboni and Schafer2015), Castellazzi et al. (Reference Castellazzi, D’Altri, de Miranda and Ubertini2017), Kasra et al. (Reference Kasra, Ali, Sara and David2018), Artese and Nico (Reference Artese and Nico2020), and Cabaleiro et al. (Reference Cabaleiro, Riveiro, Conde and Sánchez-Rodríguez2020). Such calculations require the introduction of initial imperfections in order to obtain realistic results as well as numerical convergence. Moreover, for many structures the details of imperfection shapes and magnitudes have a significant impact on the results as explored by prior studies in this area, for example, Smith (Reference Smith1975), Amlashi and Moan (Reference Amlashi and Moan2008), Benson et al. (Reference Benson, Downes, Dow, Grigoropoulos, Samuelides and Tsouvalis2009), Estefen et al. (Reference Estefen, Chujutalli and Guedes Soares2016), Ringsberg et al. (Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021), and Yi et al. (Reference Yi, Noh, Lee, Seo and Paik2021). Initial imperfections can be computed using a superposition of mode shapes determined via modal analysis or directly specified using distortion patterns based on prior surveys and engineering judgment, most commonly formulated as Fourier series. Both of these methods do not offer a direct way in which to account for realistic distortion scenarios that occur to in-service structures such as fabrication deviations or in-service damage.
Utilizing 3D Point cloud data to perform FE model updating provides initial distortion patterns and eliminates assumptions about the spatial form of initial distortions by using the actual shape of the structure. This data can be generated through a variety of techniques including laser scanning, structured light scanning, or photogrammetric reconstruction (Puliti et al., Reference Puliti, Montaggioli and Sabato2021). As a collection of Cartesian data points that depict the surface of a structure, point clouds can provide crucial information about the full 3D deformation field of target structures that is accurate to a millimeter or even submillimeter scale, for example, Gordon and Lichti (Reference Gordon and Lichti2007) and Jafari et al. (Reference Jafari, Khaloo and Lattanzi2017). This accuracy can be harnessed for high-fidelity representation of the out-of-plane deformations in FE models.
This article investigates an approach to FE model updating using interpolated point cloud data. This approach uses a kriging interpolator that accurately maps unstructured point cloud data to an FE model while inherently quantifying measurement uncertainty for envelope structural analysis. The applied methodology addresses the inherent technical limitations of spectral and modal methods. The approach is applied to nonlinear FE calculations of ultimate strength and collapse of a full-scale ship deck panel, or grillage, where both test data and point cloud data of the deformed shape are available. While the application is specific, the grillage is broadly representative of stiffened steel plate structures and the methods and conclusions are readily extended to other structural systems.
1.2. Prior work
Methods to utilize point cloud data in FE models can be broadly divided into methods that generate either a mesh or mesh-able geometry from the point data and methods that map the data to an existing FE model. Methods seeking to generate a mesh directly from point cloud data are based on using the points as the vertices of the mesh elements. This mesh is then incorporated with the solid model that represents the geometry of the structure or structural component. Poisson surface reconstruction or Delaunay triangulation are common algorithms for such approaches (Tucci et al., Reference Tucci, Bartoli, Betti, Bonora, Korumaz and Korumaz2018; Mugnai et al., Reference Mugnai, Lombardi, Tucci, Nocentini, Gigli, Fanti, Brumana, Pracchi, Rinaudo, Grimoldi, Scaioni, Previtali and Cantini2019; Bartoli et al., Reference Bartoli, Betti, Bonora, Conti, Fiorini, Kovacevic, Tesi and Tucci2020), along with voxelization-based algorithms (Castellazzi et al., Reference Castellazzi, D’Altri, Bitelli, Selvaggi and Lambertini2015, Reference Castellazzi, D’Altri, de Miranda and Ubertini2017; Guan et al., Reference Guan, Lin, Wang, Zhou, Luo and Zheng2020). Techniques that approximate the surface of a structure include mathematical approaches such as the use of b-splines or polynomial surfaces (Bertolini-Cestari et al., Reference Bertolini-Cestari, Invernizzi, Spanò and Mallamaci2012; Conde et al., Reference Conde, Ramos, Oliveira, Riveiro and Solla2017; Cuartero et al., Reference Cuartero, Cabaleiro, Sousa and Branco2019; Xu and Neumann, Reference Xu and Neumann2020; Xu and Yang, Reference Xu and Yang2020; Batar et al., Reference Batar, Tercan and Emsen2021).
While using point cloud data to directly generate FE meshes is a viable approach, it suffers from several shortcomings. One of the challenges in this mesh generation approach is having a varying point density in the point cloud data across the surface of the target structure. Commonly, FE models require adjusting the size of mesh elements in areas of high curvature to improve the accuracy of the results, which is challenging with a typically irregular point cloud. In addition to point density variances, another challenge involves having a point density that is altogether too high or too low. While a high point density is helpful for creating a fine mesh, it comes with a trade-off of a high computational cost, with the opposite being true for a low point density. Another challenge is that point clouds are typically line of sight and as such the nonmeasurable zones must be meshed separately and merged with the newly generated mesh. An analysis of this trade-off reveals that the use of mesh generation techniques requires careful planning to ensure a high enough point density for a given application, along with potentially complex data preprocessing and subsampling to ensure that the density is appropriate for the intended computational use (Wałach and Kaczmarczyk, Reference Wałach and Kaczmarczyk2021).
Likewise, even though surface representation approaches for FE modeling perform adequately in many scenarios, the approximated surfaces lose potential information provided by the point cloud data. As the point cloud is fundamentally high-resolution geometric information, the transformation from point cloud to 3D surface introduces uncertainties into the modeling process that can potentially mask small deformations. Also, while various algorithms (such as NURBS and polynomials) use the point data as control points to define the mathematical surfaces, there are still residual errors between the best-fit surface and the actual data points. The nature of these residual errors cannot be explicitly quantified, compounding model uncertainty. These same concerns exist for voxel-based approaches, as an octree level resulting in large voxels can mask important geometric information contained in a point cloud.
1.3. Contribution of this work
Established approaches for updating FE models to account for distortions using either modal analysis or assumed distortion patterns rely on assumptions regarding geometry that often do not hold true for realistic structures. The contribution of this work is a new pathway to incorporate point cloud data into FE modeling that does not rely on such assumptions. While there is prior work in using a point cloud to either represent the surface geometry or directly build the FE mesh, the presented approach uses a Gaussian process on the point cloud data to define and map the initial distortions of the FE model. This approach avoids the complexities of surface fitting or mesh-generating algorithms, and it provides equitable results compared to benchmark methods for maritime infrastructure assessment (Ringsberg et al., Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021). Additionally, it is unique compared to the methods used by benchmark participants in that it gives explicit uncertainty quantification in the form of its Gaussian variance estimate. This variance estimate ultimately provides decision-makers and infrastructure owners with quantifiable risk that is fundamentally tied to the uncertainty in the point cloud data. As point cloud data can be captured at different times during the service life of a given structure, this process can create a performance envelope of key parameters to support life cycle analysis (Okasha and Frangopol, Reference Okasha and Frangopol2010; Cheng and Frangopol, Reference Cheng and Frangopol2021).
The remainder of this article is organized as follows: Section 2 provides an overview of the methodology; Section 3 reviews a physical test specimen and associated FE model used for experimental analysis; Section 4 provides the experimental results with accompanied discussion; and Section 5 presents the conclusions.
2. Methodology
The overall goal of this methodology is to update FE models based on initial geometric distortions quantified through point cloud analysis (Figure 1). Interpolation of a local neighborhood of point cloud data centered on predefined FE mesh nodes allows for local smoothing that is suitable for numerical analysis. Gaussian process regression, commonly referred to as kriging, is employed for this interpolation, as it provides explicit uncertainty quantification in the form of a variance estimate. Given a variance interval at each nodal location, a range of deformation fields can be evaluated, resulting in an envelope analysis of mechanical performance and ultimate hull panel strength. This performance range can ultimately enable decision-makers to understand the risk associated with the current state of the structure in question, and it can better inform them of the potential necessity for repairs or remediation.
Figure 1. Overview of the methodology.
2.1. Data collection and outlier removal
Common sensors for collecting point cloud data include laser scanners, trackers, structured light scanners, and digital cameras. Regardless of the type of sensor used to create the point clouds, the data will have some number of outliers that should be removed and the methods applied herein are readily utilized regardless of the data collection method. Outlier data points are classified based on the local neighborhood structure of the data, as seen in Rusu et al. (Reference Rusu, Marton, Blodow, Dolha and Beetz2008). The average distance to the k-nearest neighbors in the defined neighborhood size is calculated for each data point in the point cloud, and the mean,
$ \mu $
, and standard deviation,
$ \sigma $
, of these distances are determined. Then, a data point is deleted when its average distance to its local neighbors is greater than the maximum allowable distance according to the following equation:
$$ \max\;\mathrm{allowable}\ \mathrm{distance}=\mu +\alpha \sigma . $$
The values for k and
$ \alpha $
are chosen empirically such that approximately 1% of the data points in the point cloud are removed. This approach, therefore, removes points that are far away from its neighbors while still preserving the geometric features of the point cloud (Rusu et al., Reference Rusu, Marton, Blodow, Dolha and Beetz2008).
2.2. Registration and deformation quantification
After data collection and outlier removal, it is necessary to register the point cloud data into a common reference frame. Situations for registration include the alignment of several smaller point clouds from multiple scanners of a common target into one global coordinate axis system (Sánchez-Aparicio et al., Reference Sánchez-Aparicio, Riveiro, González-Aguilera and Ramos2014), or the alignment of a point clouds of the same scene in different conditions (such as a structure in different load states) for change detection (Jafari et al., Reference Jafari, Khaloo and Lattanzi2017). A common algorithm for such registration is the iterative closest point algorithm (Besl and McKay, Reference Besl, McKay and Schenker1992), which iteratively minimizes the Euclidean distances of point correspondences between data sets through the rotation and translation of one point cloud against the other.
As the first step, the algorithm finds the nearest neighbor point correspondences between the points in the reference point cloud P and a data point cloud Q (Figure 2). In order to reduce the search complexity, kd-tree space partitioning is used by Friedman et al. (Reference Friedman, Bentley and Finkel1977). The algorithm then calculates the Euclidean distances of the correspondences between all points
$ {p}_i $
and
$ {q}_i $
, where
$ p\in P $
,
$ q\in Q $
, and
$ i=1,\dots, n $
correspondences. The ICP then iteratively minimizes a cost function to determine the 3 × 3 rotation matrix, R, and 3 × 1 translation vector, t, to register Q to P:
$$ E\left(\mathbf{R},\mathbf{t}\right)=\underset{R,t}{\arg \hskip0.1em \min}\sum \limits_{i=1}^n\parallel {p}_i-\left(\mathbf{R}{q}_i+\mathbf{t}\right){\parallel}_2. $$
Figure 2. Examples of point correspondences between two point clouds.
In cases where point cloud data of an entire structure (or component) is collected from a single laser scanner, it only needs to be aligned into a local coordinate axis system, commonly done by orienting along a best-fit plane. This requirement can typically be fulfilled based on the settings of the scanner without a need for further processing.
Given the registration of the point cloud data, out-of-plane deformations are characterized in the deformation quantification step. The registration in the previous step is conducted so that the deformations are quantified in a single coordinate axis that is oriented in the out-of-plane direction. In such a case, the deformation of each point is defined as the Euclidean distance orthogonal to the best-fit plane of the point cloud (Figure 3). This characterization is equivalent to a vector of directed Hausdorff distances (Huttenlocher et al., Reference Huttenlocher, Klanderman and Rucklidge1993) in a case where the reference point cloud is perfectly planar and all Euclidean distances to the points in the distorted cloud are in the surface normal direction of the planar points. This arrangement facilitates the interpolation process discussed in the following subsection.
Figure 3. Deformation quantification using directed Hausdorff distances. Euclidean distances measured between data in Q and corresponding data in the best-fit plane quantify the out-of-plane deformations.
2.3. Interpolation and uncertainty quantification
While the density of point cloud data has the potential to result in data points located very near mesh nodes, the unstructured nature of point clouds requires interpolation to cover these nodal locations reliably and uniformly. Here, kriging interpolation is used to make deformation predictions at the mesh node locations. While kriging was originally developed for geostatistical scenarios (Cressie, Reference Cressie1993), its use has expanded and can be seen in such areas as chemical and pollutant mapping in environmental sciences, as well as mechanical stress and fatigue predictions (Teixeira et al., Reference Teixeira, Nogal, O’Connor, Nichols and Dumas2019; Xu et al., Reference Xu, Bechle, Wang, Szpiro, Vedal, Bai and Marshall2019; García-Nieto et al., Reference García-Nieto, García-Gonzalo, Puig-Bargués, Duran-Ros, Ramírez de Cartagena and Arbat2020). Additionally, kriging uses the spatial autocorrelation of the underlying data to make predictions. This is relevant because it ensures that the interpolation process maintains the proper curvature of the initial distortions of the target structure.
Kriging comes in many forms, and this work uses ordinary kriging, which assumes an unknown but constant trend in the data set. This contrasts to other forms such as simple kriging and universal kriging, which use a known constant trend value or a trend function, respectively, in their modeling assumptions. Researchers have made comparisons of the performance of these kriging models, and they have shown that universal kriging generally provides better results in cases where the underlying trend is understood and can be adequately captured (Zimmerman et al., Reference Zimmerman, Pavlik, Ruggles and Armstrong1999; Altman, Reference Altman2000; Gunes et al., Reference Gunes, Cekli, Rist, Mason, Hill and Mönch2008). Within the context of structural deformations, it is possible to estimate the underlying trend of point cloud data in such cases as elastic flexure of beams. However, because there are many cases that would preclude an accurate trend model, such as damage scenarios, plastic deformations, or initial distortions from fabrication, ordinary kriging is selected for its performance on data sets with unknown trends. For a full derivation of the equations for ordinary kriging, the reader is referred to Cressie (Reference Cressie1993) and Bailey and Gatrell (Reference Bailey and Gatrell1995). However, for the purposes of this work, the primary equations are presented here.
Given the deformation observations, Z, of each data point in the point cloud, ordinary kriging makes the following model assumption:
$$ Z(s)=\mu +\varepsilon (s) $$
at locations
$ {s}_i $
for
$ i=1,2,\dots, n $
. Here,
$ \mu $
is the unknown but constant trend, while
$ \varepsilon $
is the spatially based error function. In order to make a prediction of the deformation value at a target location,
$ {s}_0 $
, ordinary kriging takes a weighted average of the surrounding n observations such that:
$$ \hat{Z}\left({s}_0\right)=\sum \limits_{i=1}^n{\lambda}_iZ\left({s}_i\right), $$
where
$ \hat{Z}\left({s}_0\right) $
is the prediction and
$ {\lambda}_i $
are the kriging weights. In ordinary kriging,
$ {\sum}_{i=1}^n{\lambda}_i=1 $
. The process seeks to minimize the mean square prediction error by solving a system of linear equations based on the partial derivatives of the Lagrangian function. The kriging weights are incorporated into the equality constraint with the Lagrangian multiplier.
In order to determine the kriging weights, the spatial autocorrelation of the data is incorporated into the Lagrangian function. It is accounted for through the semivariogram, which is:
$$ \gamma (h)=\frac{1}{2n}\sum \limits_{i=1}^n{\left(Z\left({s}_i+h\right)-Z\left({s}_i\right)\right)}^2, $$
where h is defined as the Euclidean distance between the target and the corresponding observation, while
$ \gamma $
is the semivariogram value. Thus, the system of equations in terms of the semivariogram results in the following:
$$ \lambda ={\Gamma}^{-1}\gamma, $$
where
$$ \lambda =\left[\begin{array}{c}{\lambda}_1\\ {}\vdots \\ {}{\lambda}_n\\ {}m\end{array}\right], $$
$ m\equiv $
Lagrangian multiplier
$$ \Gamma =\left[\begin{array}{cccc}{\gamma}_{11}& \cdots & {\gamma}_{1n}& 1\\ {}\vdots & \ddots & \vdots & \vdots \\ {}{\gamma}_{n1}& \cdots & {\gamma}_{nn}& 1\\ {}1& \cdots & 1& 0\end{array}\right], $$
$$ \gamma =\left[\begin{array}{c}{\gamma}_{10}\\ {}\vdots \\ {}{\gamma}_{n0}\\ {}1\end{array}\right]. $$
Solving equation (6)) provides the kriging weights to substitute into equation (4)) to calculate the prediction at a target location. Finally, the variance at the target location is determined by the following:
$$ {\sigma}^2={\lambda}^T\gamma =\sum \limits_{i=1}^n{\lambda}_i{\gamma}_{i0}+m. $$
This variance interval is used to quantify the uncertainty as seen in the next step.
2.4. Deformation field representation
After interpolation, each interpolated point location is comprised of a Gaussian distribution of the distortions at the given mesh node location. In this work, three deformation fields are derived from the prediction intervals of the distortions as each mesh node. First, the mean value of each prediction interval is used to define an “expected” distortion field. This resulting “expected” distortion field could be determined through simpler interpolation methods if uncertainty quantification was not desired. However, for this work, two additional distortion fields are extracted. One distortion field uses a determined lower bound value from each prediction interval and represents a smaller distortion field. The other field uses the upper bound value to represent a larger distortion field. The independent evaluation of each of these distortion fields, therefore, provides an envelope analysis of the given structure’s behavior that can quantify risk.
A convenient metric for this is the z-score, defined as the deviation from the distribution mean that is normalized by the standard deviation. Although other methods could be used, this metric is easily interpretable to a wide audience (as typically encountered when briefing decision-makers and infrastructure owners). Also, different confidence intervals can then be easily selected and implemented to determine a suitable range of values for quantifying the uncertainty and risk. The selection of a z-score results in three separate distortion fields: one based on the mean kriging prediction, one based on the upper bound of the confidence interval, and one based on the lower bound, as previously mentioned.
Multiple approaches were investigated to determine a suitable z-score for deformation field representation. The first of these approaches included sampling from the kriging intervals at each prediction location to determine a mean z-score. Then, Monte Carlo simulation was conducted with sampling of the distribution tails to determine suitable uncertainty bounds for further analysis. A second method followed the same approach but replaced the mean z-score of the sampled distortion field with the Mahalanobis depth (Liu et al., Reference Liu, Parelius and Singh1999; Liu and Wu, Reference Liu and Wu2017). This metric was chosen as a measure of centrality because, unlike the mean z-score, it incorporates the covariance structure of the given distortion field. However, these alternatives resulted in such tight uncertainty intervals (low z-scores) that it was decided to revert to the approach using the 95% confidence interval. Nevertheless, they could certainly be used in other situations and may be more appropriate given a different data set.
2.5. Importing into FE models
FE mesh modification is accomplished by calculating nodal displacements for each node in accordance with the interpolated distortion fields. The updated nodal positions match the distortion field locations and therefore more accurately match the actual curvature of the target structure. This is possible because the best-fit plane location established during the registration process matches the undistorted node locations of the FE model. Thus, the distortion fields in the FE model are defined directly by the output from the interpolation and uncertainty quantification process.
The nodal displacements are implemented based on a nodal shift of the mesh as opposed to a load step in the FE model. Additionally, parabolic shape functions are used such that the elements do not remain flat after the node positions have been moved to match the curvature of the physical infrastructure.
While it is possible to directly import the raw point cloud data into the FE model without interpolation, and this approach was taken during preliminary laboratory-scale studies, evaluation of strain fields between the FE models using the raw data (i.e., not interpolated), numerical inaccuracies and outliers in the strain field were observed. These observations motivated the approach presented here.
3. Experiment and Data
3.1. Grillage and test description
A full-scale ship grillage collapse experiment is utilized as the application example, as seen in Figure 4. The grillage is fabricated from welded steel plates and tee-stiffeners in both the longitudinal and transverse directions with overall dimensions of 7.5 m × 2.5 m. Stiffening in the longitudinal direction consists of three stiffeners and a single large girder; transverse stiffening consists of four transverse frames spaced 1,829 mm apart. Plating consists of two A36 steel plates welded together in the transverse direction near the midsection of the assembled structure. The grillage was outfitted with heavy end plates and side panels to allow it to be mounted into the testing apparatus. The welded plates had different thicknesses of slightly under 8 mm (Section 2) and slightly over 6 mm (Section 4), respectively. The plates are welded together in Section 3.
Figure 4. (a) Photograph and (b) diagram of the grillage. Both are reprinted from Nahshon et al. (Reference Nahshon, Shilling, Ringsberg and Darie2021) with labels added. The camera in (b) shows the photographer location in (a).
Ultimate strength testing was performed using the Carderock Division of the Naval Surface Warfare Center’s grillage test fixture, located in West Bethesda, MD. In this text fixture, one end of a grillage is fixed while the other is compressed under displacement control using 10 evenly spaced hydraulic actuators. Vertical motion of the panel is restricted via steel rod tie-downs along each side. A test panel is outfitted with numerous strain gauges that were monitored to minimize eccentric loading and monitor buckling response. Actuator displacements are tracked in order to generate an overall load-shortening curve for a grillage. Following test runs 25 and 50% of the expected ultimate capacity, a full collapse test is conducted. For the panel used in the context of this article, the peak strength of the grillage was 6.59 MN, with elastic–plastic buckling observed. Failure occurred in Section 2 of the test specimen as lateral torsional buckling of the girder, followed by local flange buckling of the stiffeners.
This same grillage experiment has been extensively evaluated as part of a round-robin benchmark study and full details of the grillage geometry and experimental setup are provided in Ringsberg et al. (Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021).
3.2. Point cloud data collection and preprocessing
Point cloud data was collected on the plates, the midspan of the web, and edges of each longitudinal stiffener flange using a Faro Vantage laser tracker. The point cloud data was preprocessed using CloudCompare (CloudCompare, 2015) software. Outlier data points were removed using the statistical outlier removal approach identified in Section 2 using a local neighborhood size of
$ k=100 $
nearest neighbors and
$ \alpha =3 $
. The approach removed slightly less than 1% of the data points from the point cloud and resulted in a point density that was approximately 1 point for every 160 mm2. This removal approach slightly underestimates the recommended number of outliers in prior work, but was deemed acceptable based on the scanning sensor and idealized environment (Rusu et al., Reference Rusu, Marton, Blodow, Dolha and Beetz2008).
The point cloud was registered with the longitudinal and transverse axes oriented in the xz-plane and all out-of-plane deformations along the y-axis. The registration and deformation quantification of the plates was established based on a best-fit plane to the scan data (Figure 5). The scan data of the webs and flanges of the stiffeners were not collected in the same reference coordinate axis as the plate data because the original purpose for its collection did not require it. To overcome this issue, the flange and web scan data was registered such that its respective centroid location matched that of the centroid location of the undistorted mesh node locations from the FE mesh.
Figure 5. Initial out-of-plane distortions of the plate for (a) measured point cloud data, (b) kriging-based interpolated data at FE mesh node locations, and (c) spectral-based data at FE mesh node locations. Note that for (c), the maximum distortion in the point cloud data is used to define a maximum distortion level, distortions are reduced in the outer bays, and distortions are not applied to the outermost bays.
4. Numerical Modeling
4.1. Meshing and solver details
All numerical modeling was performed using Abaqus 2019 FE software (Dassault Systèmes, 2019). The geometry of the FE model (Figure 6) includes the plate sections, the longitudinal and transverse stiffeners, end caps of the transverse stiffeners, and side plates. The geometry was defined based on various measurements of the components, including ultrasonic measurements of the plate thickness. Displacements on the side plates and end caps of the transverse stiffeners were fixed in the vertical direction to mimic the tie downs. Boundary conditions at the longitudinal ends of the grillage were fixed except for displacement in the axial direction at the loaded end.
Figure 6. Finite element model of the grillage panel.
The FE mesh was prepared using the AutoHull script developed by Benson et al. (Reference Benson, Downes, Dow, Grigoropoulos, Samuelides and Tsouvalis2009). AutoHull prepares the geometry and mesh from a spreadsheet definition of the plates and stiffeners in such a manner as to readily allow the specification of initial distortions and residual stresses (Benson et al., Reference Benson, Downes, Dow, Grigoropoulos, Samuelides and Tsouvalis2009). 8-node reduced integration shell elements that included five integration points through the element thickness (S8R) are used throughout. Guidelines in DNV GL AS (2013) recommend a mesh density based on expected half wavelengths of the buckled shape to properly capture buckling and local plastic deformations. As such, this model used a mesh size of 30 mm (slightly finer than what was recommended) resulting in roughly 100,000 nodes and 35,000 elements.
Material properties were defined using multi-linear stress–strain curves based on materials testing of coupons generated from the materials of the plates and stiffeners. Weld-induced residual stresses at the interface between the stiffeners and the plate sections were also included for portions of the analysis (explicitly identified in Section 4 of this work) and were implicitly assumed to be 75% of the material yield strength.
A dynamic nonlinear implicit solution solver was used to compute grillage response. This solver was selected in order to utilize inertia to regularize unstable behavior and obtain high-quality solutions without the solution noise observed in explicit methods or convergence issues associated with arc-length methods. The solver performed a forward-time integration using an implicit solution scheme utilizing the backward Euler method. In order to maintain a quasi-static response despite conducting the calculation with inertial effects dynamic effects must be damped. This is accomplished over every solution step using automatic settings available within the ABAQUS solver (*Dynamic, application = QUASI-STATIC). Sufficiently small steps must be utilized such that damping levels are sufficient to prevent introducing inertial effects into the computed solution. A convergence study was performed demonstrating that increasing the number of steps did not alter the force-displacement curve. For analyses where residual stresses are included, a static step without applied displacements or loads is run prior to the collapse analysis to ensure that the calculation is in an initial state of equilibrium. Further details on the solver may be found in Dassault Systèmes (2019).
Point cloud-based initial distortions were applied to the model using ABAQUS’s capability to read initial nodal displacements from a text file is utilized (*Imperfection). Where comparisons are made to a Benchmark solution, harmonically defined analytical initial distortion fields were applied using AutoHull functionality. The mode content and magnitudes of the Benchmark initial imperfection are based on experimental measurements. Full details are provided in Benson et al. (Reference Benson, Downes, Dow, Grigoropoulos, Samuelides and Tsouvalis2009). In both methods, no stress or strain is introduced as nodes are moved as a preprocessing step.
4.2. Interpolation and uncertainty quantification
The distortions of the mesh nodes were determined as follows: out-of-plane plate distortions in the y-direction, flange distortions in the y-direction, stiffener web distortions in the x-direction, and web tilt of the stiffeners. The coordinate axis information for the point cloud data and FE model follow Figure 6.
The point cloud data of the plate sections was interpolated onto the associated mesh node locations of the FE model based on a local neighborhood of 100 nearest neighbors. This neighborhood size was determined empirically to maintain a balance between capturing the curvature of the point cloud data and computational expense. Other neighborhood structures were not explored for this work and remain a topic for future study. Next, the distortions of the mesh nodes along the edge of each flange were interpolated based on a local neighborhood of 10 data points on the respective edge of the corresponding flange. The neighborhood was decreased to capture only the local curvature of the flange for each interpolation point. Upon completion of the interpolated values at mesh nodes along the edges of the flange, the distortions for the interior nodes were determined through linear interpolation between the associated edge nodes. The same approach was applied to the web nodes to capture the off-axis curvature (about the y-axis) of the web. Lastly, the web tilt (calculated based on measured geometries of the point cloud data) was interpolated along each unique cross-section of web nodes, and this tilt was applied to the corresponding points to capture the curvature of the data points about the z-axis (long axis).
This process was applied to create three separate distortion fields. One field was based on the mean values of the kriging prediction intervals, and the other two distortion fields were determined using the upper and lower bounds of the 95% confidence interval (z-score equal to
$ \pm 1.96 $
) at each prediction location. This prediction interval was selected after the alternative approaches using a mean z-score and Mahalanobis depth centrality metrics yielded z-score values of 0.1 and 1.0, respectively. Based on the results discussed in Section 4, the 95% confidence interval was chosen because provided a wider range of distortion values, and was therefore more illustrative.
It should be emphasized that the high quality of the point cloud data resulted in tight variance intervals. This necessitated the use of the 95% confidence interval to illuminate differences between the three distortion fields. While high-quality data will lead to small variance intervals for the distortion predictions, it naturally follows that the distortion fields will be close together.
4.3. Evaluation metrics
First, the suitability of the methodology for characterizing the initial distortions was investigated by observing the interpolated point data of the plate in the Fourier domain. In order to accommodate the requirement for uniform grid spacing, additional interpolation nodes were added at the centroid location of each mesh element. This was done strictly for the evaluation in Fourier space. This evaluation is beneficial because it enables the visualization, and ideally delineation, of both the nature of measured distortions and the level of point cloud measurement noise.
Next, the behavior of the FE model with initial imperfections mapped through this methodology is compared against the behavior of the Benchmark solution using analytical field initial imperfections. The load-end shortening (force-displacement) curves from these two models are compared against each other and against the load-end shortening curve from the physical test. Compared metrics include the ultimate capacity, failure mode, and failure location.
Last, the load-end shortening curves of the three distortion fields of the methodological FE model are compared against each other. The purpose of this comparison is to highlight the output of the uncertainty quantification scheme. This scheme is quantified by the associated range of ultimate capacity values determined from the three distortion fields defined by the kriging intervals in the uncertainty quantification step.
5. Results
5.1. Characterization in the Fourier domain
Figure 7 illustrates the measured distortions of the plate, represented in a 2D Fourier space. This representation provides a characterization of the distortions and sensor noise in the interpolated data. The signals with lower amplitudes characterize either the mode shapes that do not significantly contribute to structure behavior or the residuals between the interpolated data and the true surface (or both). Such information could be used in other approaches for modeling imperfections in these applications. The dominant frequencies are low and characterize the global curvature of the plate. But there is also higher frequency energy, distributed along the z-axis.
Figure 7. Plate distortions in the Fourier domain (DC signal removed).
5.2. Load-end shortening curve comparison
Figure 8 shows a comparison of the resulting load-end shortening curves of the FE simulations using initial distortions developed from the kriging scheme and a scheme based on a spectral analytical function. These results are shown in comparison to the results from the physical test. The spectral analysis function-based scheme serves as a benchmark approach based on the work of Ringsberg et al. (Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021). Note that all force-displacement curves are measured at the loaded end of the deck panel. Close agreement is obtained for peak load for both analyses as compared to the test data. The analysis using kriging-based initial distortions exhibits a stiffer response coupled with a sharper drop-off of load.
Figure 8. Load-end shortening curves.
The higher stiffness that both FE results exhibit is attributable to either differences in material behavior or assumptions regarding residual stresses. Material modeling errors are quite possibly due to documented sensitivities in capturing the transition of linear elastic to nonlinear plastic material behavior in the context of buckling analyses (Schafer et al., Reference Schafer, Li and Moen2010). Regarding residual stresses, it is possible that significant differences in actual versus assumed residual stress patterns exist as no measurements were collected. In order to quantify the effects of residual stresses, the model was run without residual stress (Figure 9). This approach follows the study of residual stress effects by Ringsberg et al. (Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021). While FE models can certainly provide notably different results based on modeling choices for the boundary conditions, it is expected that the boundary conditions used in this model provide a quality representation of the boundary conditions in the physical experiment (see Section 3.1). Further information about the boundary conditions of the model are seen in Ringsberg et al. (Reference Ringsberg, Darie, Nahshon, Shilling, Vaz, Benson, Brubak, Feng, Fujikubo, Gaiotti, Hu, Jang, Paik, Slagstad, Tabri, Wang, Wiegard and Yanagihara2021).
Figure 9. Load-end shortening curves from updated FE model, with and without residual stress.
While the kriging-based model demonstrates a higher stiffness than the spectral-based model, it should be noted that the kriging approach can be used in situations where the spectral approach cannot. While the spectral-based approach is adequate for scenarios involving newly constructed panels with minor surface distortions, the kriging approach would be a more feasible option in cases of large distortions that are potentially found in existing or damaged infrastructure. This approach would be beneficial for the evaluation of residual capacity of such infrastructure.
5.3. Ultimate strength and failure mode representation
A motivating goal of this study was to determine the difference in ultimate capacity of the deck panel predicted by FE analysis using imperfection mapping of point cloud data. Table 1 presents a summary of ultimate capacity and relative error of the FE prediction using the kriging as well as the Benchmark methodologies. Both approaches show strong agreement with the ultimate capacity of the physical test specimen, albeit at a noticeably shorter displacements. This is noteworthy due to the behavior of the system in the physical test reaching its ultimate capacity after slightly more than 5 mm difference of displacement than those of the simulations, whose responses reflect more traditional collapse behavior.
Table 1. Comparison of ultimate capacities.
This difference in response is likely attributed to system behavior in the physical panel, as the force-displacement curves of the simulations demonstrate a steeper elastic slope (a surrogate for higher stiffness) and more abrupt drop to a residual capacity after buckling when compared to the physical deck panel. As this panel transitions from linear to nonlinear behavior, prior research has demonstrated the modeling sensitivity to material properties (Schafer et al., Reference Schafer, Li and Moen2010). Thus, a closer match to the exact response of the physical panel using these modeling techniques is unlikely, as the FE method is inherently stiff in collapse scenarios compared to measured physical behaviors. Nevertheless, the ultimate capacity of the physical system is closely met by the ultimate capacity from the simulated models with small relative errors. It is also noteworthy that even though the load-end shortening curve from the results of the methodology show a noticeable difference from that of the benchmark, the ultimate capacities are closely matched.
Unlike the predicted ultimate capacity, the location of the buckling failure in the physical test is not accurately replicated by the FE simulation. The failure of the physical panel occurred in the second section of the panel, but the failure in the FE simulation of the methodology occurred in the fourth section (Figure 10). For the benchmark simulation, failure occurred in the third section.
Figure 10. Failure mechanism and location comparison between (a) the physical test (reprinted from Nahshon et al. (Reference Nahshon, Shilling, Ringsberg and Darie2021)) and (b) the finite element model simulation. Camera in (b) shows photographer location in (a).
The plate was thicker in Section 2 (just under 8 mm) than it was in Section 4 (just over 6 mm), which would suggest that failure would be more likely in Section 4. The failure mechanisms of the FE model based on this methodology included lateral torsional buckling of the girder with buckling in the plate, and local flange buckling of the other stiffeners, all of which occurred in Section 4. While these behaviors occurred in the physical test, as well, they did not occur in the same location.
The initial out-of-plane distortions of the plate were more significant in Section 2 when compared to Section 4. The fact that failure did not occur in the section with the largest initial distortions indicates that, in this particular case, a combination of additional factors such as unanticipated material behavior or load cell distribution contributing significantly to the behavior of this complex system.
5.4. Uncertainty quantification
When evaluating the point cloud uncertainty quantification, the results show a tight interval that is not readily evident when observing the entirety of the displacement range (Figure 11a). This relatively tight interval suggests that the point cloud data is of high fidelity with low noise, also indicated by the low kriging variance of the prediction intervals in the interpolation results.
Figure 11. Simulation results with uncertainty bounds for (a) the entire performance range, (b) a portion of the elastic behavior region, (c) ultimate capacity, and (d) a portion of the plastic behavior region.
A closer inspection of the force-displacement curve in the linear elastic region (Figure 11b) reveals symmetrical uncertainty bounds. This supports the implementation of the uncertainty quantification approach for this type of system behavior because the uncertainty captures a performance range that is bifurcated by the performance resulting from the deformation field defined by the mean kriging results of the interpolation process. In this figure, the upper bound of the confidence interval refers to the upper bound of the 95% confidence interval from the kriging prediction, and it displays below the results from the expected simulation performance (solid line) while the lower bound displays above the expected performance.
This relationship happens because the higher values from the kriging intervals at the node locations result in larger distortions, which is synonymous with higher amplitudes of initial distortions that are present in the deck panels. The resulting behavior shows a decrease in the stiffness of the deck panel system, thus resulting in a force-displacement curve that is lower than the curve based on the mean kriging predictions. The inverse is true for the results from the lower bound of the kriging intervals.
For uncertainty quantification of the ultimate capacity (Figure 11c), the values reported in Table 2 are taken from the maximum reported value of the simulation using the deformation fields defined by the lower bound of the kriging interval and the mean value of the kriging interval.
Table 2. Range of ultimate capacity based on uncertainty quantification.
Figure 11d illustrates how envelope analyses of initial distortions do not always follow consistent trends. This shift in the envelope system response occurs post-buckling, and so it is likely due to localized plasticity in the simulation. What this envelope behavior highlights is that localized distortions in a structural system can result in complex and unpredictable system-level response, particularly when simulating elasto-plastic behavior.
6. Conclusion
This work presented a methodology for using point cloud data in FE model updating. Gaussian process regression (kriging) is used to interpolate distortions quantified in a point cloud at nodal locations in an independently developed FE mesh. This approach reduces numerical inaccuracies that can occur from the direct meshing of point cloud data, as has been done in prior work. This approach also allows for direct uncertainty quantification, as the variance intervals from the kriging interpolation provide the information necessary for envelope structural analysis. This capability can provide a new approach to risk-based structural assessment.
Additionally, this approach is compared to a spectral-based approach for mapping initial distortions into FE models. While the kriging approach includes complexities of capturing spatial autocorrelation in the data, its benefit is realized in situations involving the evaluation of in-service infrastructure that may have more significant and unpredictable distortion patterns than those distortions that can be captured in the spectral approach.
The approach was demonstrated using point cloud data collected from a full-scale grillage panel representative of a ship hull, including initial distortions. Using the presented updating methodology, the results were used to update a FE simulation and were compared against the experimental results and benchmark results. The updated FE model showed strong agreement with the benchmark results and captured the ultimate strength of the deck panel with a relative error of less than 2% when compared to the experimental test. Additionally, it demonstrated the same failure mechanisms as the physical test, including lateral torsional buckling of the main girder followed by local flange buckling of the remaining stiffeners, and then buckling of the plate section. However, while the failure mechanism was the same, the simulation did not capture the correct location of the initial local buckling failure in the panel.
Future avenues of research include the exploration of different interpolation schemes and uncertainty quantification metrics. The kriging process requires the determination and modeling of the spatial autocorrelation of the data, and a more automated and less empirical approach to this process is warranted. However, this approach is currently planned for evaluating FE models of press-fit shafts to determine its suitability for infrastructure components with more significant distortions that those seen in this work.
Author contribution
Conceptualization: W.G., K.N., D.L.; Data curation: W.G., K.N.; Formal analysis: W.G., K.N.; Funding acquisition: K.N., D.L.; Investigation: W.G., K.N., K.A.; Methodology: W.G.; Project administration: W.G., D.L.; Resources: K.N.; Software: W.G., K.N., K.A.; Supervision: D.L.; Validation: W.G., K.N., K.A., D.L.; Visualization: W.G., K.N.; Writing—original draft: W.G.; Writing—reviewing and editing: W.G., K.N., K.A., D.L.; All authors approved the final submitted draft.
Competing interest
The authors declare that they have no competing interests exist.
Data availability statement
Some data used during the study (laser scanner point cloud data) are available in a repository online in accordance with funder data retention policies (Nahshon et al., Reference Nahshon, Shilling, Ringsberg and Darie2021) available from https://doi.org/10.17632/bsh4gmg8s5.1. Some data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.
Funding statement
Portions of this work were performed by W.G., K.A., and D.L. and were supported by a grant from the Office of Naval Research (No. N00014-18-1-2014) and under the Structural Reliability Program overseen by Dr. Paul E. Hess III. K.N. was fully supported by funding under this same program. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Office of Naval Research.
Open access
Comments
No Comments have been published for this article.