The complex geometry of structural components result s in uneven stress distribution in structures under loading. The regions with stress concentrations often act as critical locations from where damages can initiate and propagate under different types of loading . Accurate estimation of stress distribution , especially at such critical loca tions , is vital for a more reliable prediction of possible damage prognosis or remaining useful life (RUL) estimation of the structure or a component . Traditional sensing methods often provide response measurements at a few localized points and demand sensor deployment in large numbers to get a distributed response. This may yet be sparse for the steep changes in stress at critical locations. In this study, we propose a hybrid data + model -based submodeling (HDMS) method to achieve a refined estimate of distributed structural response in and around the critical locations . The HDMS method uses just the measured response on the pre -selected boundaries around the location s of interest, as input to drive the corresponding submodel of a structural component or a connection , given it s geometry and material properties are known . The performance of the HDMS method in response estimation is demonstrated through two vertically loaded plates , one with two holes connected by a slit and the other with two wide slits, respectively. The refined response estimated by HDMS could determine asymmetric response, nonlinear behavior , and permanent set at the critical locations with an aver age error less than 50 𝜇strain at higher load stages, making HDMS a versatile method for refined response estimation . KEYWORDS Hybrid submodeling , Response estimation, Digital Image Correlation, Finite Element Submodeling, Force -based submodeling , Full -field response 1 Introduction Health monitoring and prognosis are the two most important steps of the condition -based maintenance (CBM) method , which is widely accepted as an economical and efficient maintenance method compared to predictive and preven tative maintenance methods in many fields1,2. In the CBM approach , the maintenance decisions often depend on t he remaining useful life of the structure o r structural component .3– 5. In this context, the remaining useful life is synonymous with remaining fatigue life , especially in metallic structures. Complex geometric features in structures such as holes, slits, a sharp change in geometry, welded joints, connections, etc., act as critical locations with s tress concentrations, from where cracks can initiate and propagate under fatigue or other kinds of loading. Statistical methods traditionally use nominal stress or local strain approaches to estimate the fatigue life before the crack initiation at a critic al location . Each of these statistical methods require s structural response information under varying loading as a primary variable for fatigue life prediction6. The nominal stress and local strain -based methods are total life approaches that use the cyclic stress range (𝑆) and plastic or total strain range (𝜖) A preprint - July 1 4, 2020 2 with cycle counts (𝑁) to estimate the fatigue life using rain-flow counting methods and 𝑆−𝑁 and 𝜖−𝑁 curves respectively . The Miner’s rules exten ds this approach to variable amplitude loading7. The accuracy of the cycle count depends on the accuracy of peaks in the stress or strain response history. The peak stress locations in a structure and stre ss distributions around critical locations can vary with the changing load scenarios. Accurate estimation of stress distribution in and around the critical locations under varying load scenarios as experienced in real structures is of interest to estimate th e fatigue life of the structure. Response at desired locations of a structure can be measured or estimated through direct and indirect sensing methods, respectively. Direct sensing methods can only be used to measure local phenomena limited to the locatio ns of sensor deployment. This requires deployment of a large number of sensors to gauge the response in at all possible critical locations with a dense network of sensors. Distributed strain sensors12,13, distribute d optical fiber sensors14, Digital image correlati on (DIC) methods15–18, self -sensing materials , and different types of sensing skins19–22 are employed for distributed response measurement purposes . These methods are relatively new, expensive to scale up their deploy ment to real structures and manage the recorded data . Even with such sensing methods, the critical locations with peak stress can be just inaccessible to physical sensing23. Indirect and virtual sensing methods are developed for full-field or localized response estimation using available response measurements from other locations. Limongelli24 estimated acceleration response along the heig ht of building applying spline interpolation to measured responses of optimally located sensors. Response data from selected locations were used in conjunction with finite el ement (FE) model to obtain full -field dynamic strain response , applying a variety of methods such as Finite Element Approximation and Local Diffuse Approximation25, Empirical Mode Decomposition26, modal superposition27, real-time data expansion , and finite elemen t stress -strain transformation equations28. Several methods are developed for optimal sensor placement instead of dense sensor networks to reduce redundant data and improve observed information 8–11. Petersen et al .29 estimated full -field dynamic response of a long span floating bridge employing a joint- input state estimation algorithm and Kalman filter using acceleration measurements. Wang et al.30 built Zernike moment descriptors to estimate full -field data as a function of load by training with Zernike polynomials and full -field strain data. Kullaa31 studied analytical virtual sensing by applying the Craig - Brampton sub-structuring method only to the desired part of the structure . Kullaa32 introduced a combined empirical and analytical virtual sensing m ethod for full -field dynamic response estimation through a Bayesian approach. Lall et al.33 were the first to d rive a submodel with sensor information for load -free board assemblies subject to impact. The velocity components from digital image correlation measurements at the cohesive joint were used as a boundary condition at the attachment degree of freedom to drive the submodel by displacement -based submodeling to estimate the response . There is a shortage of method s that can estimate refined stress/strain distribution , including nonlinearity around a critical location for any given loading on the actual structure. This study proposes an HDMS method to estimate spatially refined response distribution near critical locations of structures . The HDMS method uses sensor measurements on the pre -selected boundaries of a region of interest in a structure to drive the submodel by applying force -based submodeling procedure . The following are the main objectives and contributions of the paper: • To propose a novel f ramework of the HDMS method , which estimates refined response near critical locations without the knowledge of global ge ometry and actual loads on the structure • Successful implementation of HDMS as a force -based submodel of a subregion , driven by similar or a combination of different types of sensor measurements on the subregion boundary away from critical locations A preprint - July 1 4, 2020 3 • Experiemntal verification of performance of the HDMS method through two cases of critical locations in plate structures • To demonstrate the successful estimation of both linear and nonlinear behavior near critical locations using the HDMS method The remainder of the paper is organized as follows. First, FE submodeling and its variations are discussed. Second, the HDMS framework is explained in detail. Third, experimental verification procedure and use of the Digital Image Correlation (DIC) method for measurement of boundary response are described . Fourth, a case study for the application of the HDMS method for response estimation around holes in a plate structure in tension connected by a slit is discuss ed. Fifth, another case study of response estimation around a slit in a plate loaded in tension using HDMS is discuss ed. Sixth, the study is summarized , and conclusions drawn from this study are presented. 2 Hybrid data + model -based submodeling 2.1 Submodeling Submodel is an independent, refined FE model of a region of interest containing features such as holes, fillets, contact interfaces or connections in a global FE model of a structure (Figure 1 (a)). The critical features of the sub-region that are simplified in the global model , such as geometric details can be defined more accurately through submodeling . Submodeling allows using a finer mesh , as shown in Figure 1 (b), to obtain refined results in the critical region , which regulates the discretization errors of the coarse mesh FE models, with a considerably reduced computational burden compared to the full FE model of fine mesh34. Despite increasing computational capabilities, submodeling is still relevant . Submodeling is employed in a variety of challenging tasks such as simple stress concentrations35, 3D crack propa gation in complex welded structures38, wheel /rail rolling contact analysis36, analyzing the mechanical condition of a central buckle of a long -span suspension bridge37, fatigue analysis of welded rib to floor -beam connections of an orthotropic bridge deck39 and more40. Figure 1 (a) FE Submodeling, (b) Mesh refinement in a submodel A submodel is driven by the response of the globa l model measured at the submodel ’s cut-boundary , as depicted in Figure 1 (a). A finite element submodel can be defined as follows : 𝐾𝑠𝑈𝑠=𝑃𝑠, (1) where 𝐾𝑠 is the full stiffness matrix of the submodel composed of element matrices 𝑘𝑒𝑠, 𝑈𝑠 is the response vector of the submodel composed element degrees of freedom 𝑢𝑒𝑠, and 𝑃𝑠 is the distributed load vector of A preprint - July 1 4, 2020 4 the submodel composed of distributed load on the faces of each element 𝑝𝑒. The element -wise finite element model of the submodel is written , as shown in Equation 2: 𝑘𝑒𝑠𝑢𝑒𝑠=𝑝𝑒𝑠, (2) where 𝑝𝑒𝑠=0 for internal elements of a submodel. There are two primary methods of submodeling. When the submodel is driven by the nodal displacements of the global model at the cut-boundary of the submodel region , it is called displacement -based or node -based submodeling. When stresses or traction from the integration points of boundary surface elements are used to drive the submodel , it is called stress -based or force -based submodeling . Let, 𝑢𝑔1,…,𝑢𝑔𝑛 are the displacements at nodes 1,…,𝑛 of the global model on a cut-boundary. The displacement boundary conditions 𝑢𝑠1,…,𝑢𝑠𝑚 of the 𝑚 submodel nodes 1,…,𝑚 that lie on the same cut -boundary between 1 and 𝑛 are determined using a an interpolation function with 𝑢𝑔1,…,𝑢𝑔𝑛 at the 𝑛 nodes as the control points . The displacements are used as essential boundary conditions in Equation 1 via 𝑢𝑒𝑠 in Equation 2 to obtain a refined response of the submodel. Likewise, in the case of force -based submodeling, the interpolation functions are built from strain vectors 𝜀𝑔1,…,𝜀𝑔𝑘 containing normal and shear strains sampled at the integration points of global model elements 1,…,𝑘 on the cut - boundary. The strain vectors 𝜀𝑠1,…,𝜀𝑠𝑙 at the integration points of 𝑙 submodel elements 1,…,𝑙 on the same cut-boundary between 1 and 𝑘 are obtained from the global strain vector in terpolation functions. Assuming linear behavior of the at the cut -boundaries, normal (𝜀𝑛) and shear (𝜀𝑠ℎ) strains interpolated to the submodel boundary elements are converted to normal (𝜎𝑛) and shear (𝜎𝑠ℎ) surface traction on the boundary element faces using the following constitutive Equation : [𝜎𝑛 𝜎𝑠ℎ]=[𝐸0 0𝐺](𝜀𝑛 𝜀𝑠ℎ), (3) wher e E is Young’s modulus , and G is the shear modulus. The normal and shear traction values obtained from Equation 3 are used to bu ild 𝑝𝑒𝑠 which is in turn used to bu ild the distributed load vector 𝑃𝑠 on the boundaries of the submodel and a refined response from the submodel 𝑈𝑠 is obtained. A conventional submodel analysis requires both the global model and the local submodel where the global model is solved first , and the results at the boundaries are used to drive the submodel . Initially developed as a finite element zooming method41, the displacement -based submodeling has be en extensively studied for application in the analysis of stress concentrations35,42,43. The nonlinearity and fracture mechanics associated with c ritical locations were also explored in various studies34,44,45. The cut-boundary displacements are valid only if the stiffness of the refined submodel does not change46. In such cases , surface -based submodeling is also a feasible solution except that it gives rigid body displacements when no dis placement boundary conditions are defined . This issue can be resolved using mi xed boundary conditions or by inertia relief , which applies a set of rotational and translational accelerations at chosen points such as the center of gravity to counterbalance the unbalanced external forces. Various commercial finite element software packages such as ABAQUS 2017 , MSC/NASTRAN 2009 , and ANSYS 2020 provide options for automatic inertia relief , which is commonly used in the analysis of unsupported structures or systems like aerospace structures50. However , it is not advisable to compute submodel displacements in these cases. 2.2 Hybrid data +model -based submodeling framework The proposed HDMS method is a hybrid submodeling method that integrates the finite element submodel and the measured structural response at the submodel boundary to estimate a distributed refined response around critical locations in real structures. Unlike the regular submodeling methods, the HDMS method does not require both global and local models and the actual loads acting on the actual structures. A refined response can be estimated using only the FE model of the submodel region with known material A preprint - July 1 4, 2020 5 properties , geometry , and sparse boundary response measurements from the actual structure . A finite element submodel of the subregion with critical locations is driven by the normal and shear stress es or displacements measured from the actual structure along the submodel boundary , using sensors instead of the boundary response of the global fini te element model (see Figure 2). The submodel mesh size near the critical region is chosen such tha t a tradeoff is achieved between discretization error and the computational time while maintaining a good aspect ratio of the submodel elements . This is further demonstrated through the case studies presented in the following sections . The cut-boundaries of the submodel should be chosen at a sufficient distance from the cr itical locations or any other local effects such as stress concentrations in the actual structure . The local effects coinciding with the boundaries would result in additional driving forces on the submodel , that would otherwise not affect the response at a critical location in the actual structure47. It has b een established by St.Venenat’s principle that the submodel stress distribution near stress concentrations can be improved to a certain limit by keeping the submodel boundaries away from the stress concentration51, and this has been prove n for displacement -based submodeling around a stress concentration52. As a result, i t is important to deploy the sensors (or select measurement locations) sufficiently away from the stress concentrations inside or outside of the submodel. T he sensors are located on submodel boundaries that are selected to avoid any steep stress gradients or material nonlinearity. Narvydas a nd Puodziuniene46 demonstrated that force - based submodelin g could give satisfactory results under moderate global model mesh density. In the case of HDMS, the global mesh density corresponds to the spatial density of sensor measurements on the boundary of the submodel region. Both force -based submodeling and displacement -based submodeling are suitable to employ in HDMS depending on the type of response available on the submodel boundaries . The nodal displacements or element stresses from the global model are to be interpolated according to the refined mesh of the submodel boundary with the least possible boundary error. Cubic spline interpolation is employed to minimize the boundary condition error, as suggested by Cormier et al.35. Let 𝑧1,…,𝑧𝑞 be the locations of 𝑞 strain or displacement sensors on one of the cut -boundaries of the HDMS submodel. The strain sensor measurements 𝜀𝑧1,…,𝜀𝑧𝑞 and displacement measurements 𝑢𝑧1,…,𝑢𝑧𝑞 are used to build respective piecewise cubic spline interpolation functions 𝑆(𝑧) (Equation 4). A one-dimensional cubic spline is const ructed of 𝑞−1 piecewise 3rd-degree polynomials 𝐶𝑖(𝑧) passing through 𝑧1,…,𝑧𝑞 acting as control points. 𝑆(𝑧)= { 𝐶1(𝑧), 𝑧1≤𝑧≤𝑧2⋯ 𝐶𝑖(𝑧), 𝑧𝑖≤𝑧≤𝑧𝑖+1 ⋯ 𝐶𝑞−1(𝑧), 𝑧𝑞−1≤𝑧≤𝑧𝑞, (4) where 𝑆(𝑧𝑖) are the sensor measurements. The coefficients of each of the polynomials can be estimated through the sensor measuremen ts at control points. The respective functions are further employed to determin e interpolated force boundary conditions at the submodel cut -boundary elements 1,…,𝑙 or displacement s at submodel cut-boundary nodes 1,…,𝑚 of the refined HDMS submodel cut -boundary respectively with cubic -spline interpolation. The interpolated strain s are converted to stress ( Equation 3) and applied to each of the boundary element faces as traction loading. The element s on the FE submodel cut-boundaries and their corresponding element faces in ABAQUS 2017 are identi fied according to their location s. The submodel can now be solved as an independent model driven by the applied traction boundary loading or essential boundary conditions through Equation 1. A preprint - July 1 4, 2020 6 Submodeling is a fundamental boundary value problem of a bounded body loaded at the surface with traction , displacement , or a combination of both . Strain (𝜀) measurements on the boundaries are cheaper and more feasible to obtain compared to the displacement measurements on the submodel boundaries in real structures . In such a case , in the absence of any displacement boundary conditions, the HDMS becomes a traction boundary value problem. If the global model has a closed -form solution within the theory of elasticity, we can apply stresses from the boundary of an interior region as surface traction to a corresponding submodel analysis . According to the uniqueness theorem of elastostatics, a suitabl y refined submodel shall converge to a unique exact solution , limited by the accu racy of the boundary conditions35. That is for a body subjected to a set of mixed boundary loadings ; there can be one an d only one set of stress distribution 53–55. Kirchhoff’s theorem provides a theoretical confirmation of the uniqueness of the solution of the boundary value problem by stating that, ‘In the traction boundary value problem , there is a uniqueness to within a rigid body displacement. ’ The traction from sensor measurements on the cut - boundaries , when applied on FE submodel boundaries, does not necessarily balance out . There exist residual unbalanced forces causing rigid body displacements that can be eliminated by using inertial relief , as explained in Section 2.1. Figure 2 Hybrid data+ model -based submodeling 3 Experimental validation The proposed HDMS method is implemented and validated in this work through two experimental case studies of plate structures containing different critical features . In the first case study, the stress distribution is estimated around the holes connected by a slit in a steel plate under tensi le loading . In the second case study, the stress distribution is determined around the slits acting as critical locations in a steel plate under static tensile loading. 3.1 Data Acquisition for HDMS The response of the actual structure measured at the cut -boundaries sho uld be spatially dense enough to minimize the error due to boundary conditions driving the HDMS . In this study, combinations of spatially dense DIC and strain gauges are used to acquire submodel boundary response from the plate structures and also verify the response estimated by HDMS around the critical locations. . The strain measurements from DIC are used in this study instead of displacement measurements to establish that HDMS can be applied to real structures where obtaining strain measurements is cheaper and more feasible than displacement measurements , especially when the submodel boundaries are inaccessible . Here, t he DIC A preprint - July 1 4, 2020 7 measur ements at submodel cut -boundaries are used as discrete sensors. It should be noted that the availability of spatially dense DIC measurements only allowed us to cut the cost of multiple strain sensors in the study and validate the HDMS method but not essential for implementation. In the case of the plate with two slits , the submodel boundary res ponse data necessary to perform HDMS is acquired entirely through DIC. In the case of the plate with holes connected by a slit , the boundary response data is obtained with a combination of strain gauge s and DIC measurements. The following section briefly describes the technique of DIC and its set up for acquiring surface response. DIC is a non-contact optical method to measure 3D coordinates for the evaluation of deformation that is being increasingly employed in various industries , including Civil Engineering56. DIC technique tracks the changes in grey value patterns in small neighborhoods called subsets on the surface of the deforming body comparing to an initial reference image to measure the 2D and 3D deformation vector fields and strain maps . Random unique pixel blocks with a range of contrast and intensity levels required for precise measurement is achieved by spray painting the test specimen’s surface with a random stochastic pattern of black dots on a contrast white painted background57. The images are obtained in special LED lighting with high-resolution cameras , as shown in Figure 6 (a). The position of sensor setup for the test surface and the definitions of measuring distance and measuring volume are illustrated in Figure 3 . The camera resolution , focal length , and distance from the measured surface are chosen based on the measuring volume. The camera aperture and exposure are selected according to the desired frame rate and lighting conditions . Such f lexibility makes the DIC technique suitable for recording measurements for wide range s of speeds58 and measuring volume sizes56,59,60. For both experimental studies in the paper , test object images at different stages of deformation are captured by 12M P sensors on an Adjustable base ARAMIS 3D measuring system. The system has two cameras set at an angle (𝛼) of 24° fitted with 24 mm focal length lenses , and the aperture in both came ras is adjusted to the intensity of the LED lighting of the system and desired measur ement speed. The images acquired are processed in DIC software GOM Correlate developed by the Trilion Quality Systems to estimate the strain field61. Figure 3 Top view of a representative DIC sensor set up with respect to monitored object A preprint - July 1 4, 2020 8 4 Case study 1: Plate with holes connected by a slit This section describes the proposed theory of the HDMS method through the case study of a plate with holes connected by a thin slit and subjected to uniaxial tensile loading . The test description, implementation of the HDMS method , and the results refined response estimate at edges of holes are presented in the following sections. 4.1 Test description The steel plate is 304.8 mm (12 in) wide, 609.6 mm (24 in) long , and 12.7 mm (0.5 in) thick , as shown in Figure 4 (a). A 76.2 mm (3 in) long horizontal slit is located at the cent er of the plate connecting two holes of 2.5 in diameter. The setup is intended to mimic fatigue cracks in bridge steel members arrested by drilling holes. The plate is clamped at top and bottom , restricting vertical and out of plane displacements at the restraints under tensile loading , as shown in Figure 6 (a). The plate is loaded and unloaded between a maximum of 204.72 kN (46 kips) and a minimum of 20. 47 kN (4.6 kips) in tension following the profile shown in Figure 5 . A preliminary analysis of the global FE model of the full plate under a static uniaxial tensile load is conducted to understand the qualitative stress distribution of the plate. The plate is modeled with solid eight -node linear elem ents with reduced integration in ABAQUS 2017 . A mesh size of 0.05 in ele ments is used to model edges of the plate holes. The front and back surfaces of the plate model in contact with each clamp are constrained in all degrees of free dom and coupled with ABAQUS kinematic coupling and connected to reference points located at the midpoints of top and bottom surfaces of the plate . The plate is loaded symmetrically through the reference points. The constraints enforced by kinematic coupling mimicking the cl amps can be seen as dark blue square patches of low vertical strains in Figure 4 (b). Analysis of the global FE model of the full plate under uniaxial load shows high tensile stresses on the edges of the holes at mid -thickness of the plate, as shown in Figure 4 (b), making these the critical locations. As a result, strain gauges SG1, SG2 are placed in these locations with high chances for crack initiation and growth with further fatigue loading. Additional strain gauges are located on the top and bottom boundaries of the submodel. SG3 and SG5 are at the midpoint s of top boundary on the front and back faces of the plate , respectively . Likewise , SG4 and SG6 are located on the bottom boundary , as shown in Figure 4 (a). The plate ’s strain response distribution on one side (front) is monitored using the gauges and DIC equipment complementing each other . Figure 6 (a) shows the test setup with DIC cameras perpendicular to the sp eckled plat e surface. The strain gauges SG1, SG2, SG3 , and SG4, are also seen in respective locations covered by masking tape. The images of the plate under loading were captured at a frequency of 2 Hz throughout loading and unloading. These images are post -processed to estimate the strain distribution on the plate sur face. Figure 6 (b) shows the colormap of vertical strain, 𝜀𝑦 distribution in the region within the measuring volume on the surface of the plate facing the DIC cameras . The color bar on the right side of the colormap shows the corresponding scale of micro -strain (𝜇strain ) values. It is observe d that the highest strain values are close to the outer edges of the holes, which is the expected critical locations with stress concentrations. The discontinuities in the 𝜀𝑦 distribution colormap correspond to the areas where speckle pattern is invisible to the cameras. These include structural discontinuities holes , the slit , and the regions covered by strain gauges , masking tape , and wires . The strain distribution recorded by DIC clearly shows a pattern similar to the Global FE plate model , as shown in Figure 4 (b). A preprint - July 1 4, 2020 9 Figure 4 (a) Plate with holes connected by a slit (all dimensions shown in mm ), (b) Strain distribution in global FE model of the plate with holes under uniaxial tension, (c) Location of strain gauges and submodel cut -boundaries Figure 5 Uniaxial tensile load profile for the plate specimen with holes at ti me instants chosen for HDMS analysis A preprint - July 1 4, 2020 10 Figure 6 (a) Test setup for the plate with holes monitored by DIC cameras and strain gauges , (b) Vertical strain distribution of the plate with holes measured by DIC 4.2 Implementation of HDMS Sub model A 3D FEM submodel of 304.8 mm (12 in) × 304.8 mm (12 in) × 12.7 mm (0.5 in) is built with slit connecting holes located at the center using ABAQUS 2017 as shown in Figure 7 . The submodel region and the corresponding cut -boundaries are shown in Figure 4 (c) as a shaded region on the global plate sketch. Elastic and nonlinear submodels are generated using A572 grade steel properties with a yi eld stress of 380 MPa (55 ksi) . Solid eight -node linear elements with reduced integration are used. A mesh size sensitivity analysis is conducted with a uniformly distributed boundary load of 0.35 kN mm-1 (2 kip in-1) to find a suitable mesh size around the edges of the holes while keeping the aspect ratio and input file runtime low. The input file run time, strain response a t SG1, and SG2 locat ions for submodels with different mesh sizes are compared. Table 1 displays the runtime, strain at SG1 location, corresponding % changes , and solid element aspect ratio near the holes. It is observed from Table 1 and Figure 8 (a) and (b) that the strain values at SG1 change at a constant rate with reducing mesh size around the holes . The input file runtime for each mesh size is observed to increase gradually until 1.27 mm (0.05 in) and suddenly increases by 480.15 % at 0.635 mm (0.025 in), for only a 2.9 % change in 𝜀𝑦. A mesh size of 1.27 mm (0. 05 in) is used around the edges of the holes considering a shorter run time and an element thickness to width ratio of 2.5 to perform HDMS analys is efficiently at multiple instants of loading . The chosen mesh density near and around the holes in the submodel is displayed in Figure 7 . It should be noted that p rocessors with better performance can be used to employ smaller mesh size and increase the accuracy of HDMS while keeping the runtime short. Moreo ver, the purpose of this study is to demonstrate the methodology of HDMS primarily . A preprint - July 1 4, 2020 11 Table 1 Comparison of strain at SG1 and CPU time to run the plate with holes submodel for different mesh sizes Mesh size at holes, mm. (in.) Run time, (sec) Run time, (% change) 𝜺𝒚 at SG1, (𝝁strain) 𝜺𝒚 at SG1, (% change) Element thickness/ element width 2.54 (0.1) 7.40 -- 569.46 -- 3.175/1.905 = 0.833 1.905 ( 0.075 ) 9.70 31.08 585.01 2.74 3.175/1.905 = 1.67 1.5875 ( 0.0625 ) 16.20 67.01 593.04 1.37 3.175/1.5875 = 2 1.27 ( 0.05) 26.20 61.73 602.30 1.39 3.175/1.27 = 2.5 0.635 ( 0.025) 1520 480.15 618.75 2.90 3.175/0.635 = 5 Figure 7 Mesh density distribution of submodel of the plate with holes with 1.27 mm (0.05 in.) mesh size around the holes Figure 8 (a) Comparison of input file runtime , (b) 𝜀𝑦 at SG1 for different mesh sizes around the holes A preprint - July 1 4, 2020 12 The top and bottom edges of the submodel are the only cut -boundaries of th e submodel , and they pass through the locations of SG3 and SG 4, respectively , as shown in Figure 4 (c). The values of traction normal to the boundary 𝜀𝑦 and shear 𝜀𝑥𝑦 measured on these boundaries through DIC , and the strain gauges are interpolated with cubic spline interpolation to each element on the submodel cut-boundaries . As the cut-boundaries lie in the elastic region throughout the loading , the interpolated strains are converted to stress using linear constitut ive relation, Equation 3. The strain gauge recordings at SG3, SG5 on the top boundary , and SG4, SG6 on the bottom boundary are compared , as shown in Figure 9 . An average difference of 15 𝜇strain is observed during maximum loading between SG3 and SG5 readings on top submodel boundary , which was not observed between SG4 and SG6 on the bottom boundary. This can be because of slight bending in the plate experience by the plate during loading. As the diff erence is as small as 15 𝜇strain we neglected the slight bending in the plate while loading. In addition to th is, as the boundaries are away from the stress concentrations , we can a ssum e absence out of plane stresses at these locations and that the strain s measured on the surface of the plate with DIC at the boundaries do not change through the thickness of the plate . The s ame stress is applied to boundary elements through the plate thickness . The HDMS analysis is performed for load stages at every 5 seconds interval to estimate the strain distribution in the submodel. A nonlinear global FE model with the same mesh density in the submodel region is created at the corresponding load stages to compare wit h the performance of the HDMS method. Figure 9 Comparison of st rain gauge recordings between front and back of the plate (a) Comparison strain gauges SG3 and SG5 on the top submodel boundary, (b) Comparison strain gauges SG4 and SG6 on the bottom submodel boundary 4.3 Results of Case Study 1 Figure 10 show s the strain distribution s 𝜀𝑦 near the critical locations of the plate obtained from DIC measurement and NL submodel estimation at a load stage of 204.72 kN. A similar pattern of strain distribution is observed in the DIC measurement and the nonlinear HDMS submodel. The NL submodel 𝜀𝑦 distribution points are same as the NL submodel grid. It can be observed that the spatial refinement of strains estimated through the HDMS method depends on the mesh refinement of the submodel . The white strips in the color ma p are locations of the plate where speckle patter ns are not visible to the DIC cameras due to the holes and co nnecting slit, wires connecting strain g auges, and masking tape covering the strain gauges at critical locations . A preprint - July 1 4, 2020 13 Figure 10 Distributions of strain obtained from DIC measurements and NL hybrid submodel For a better comparison , the 𝜀𝑦measured and estimated 𝜀𝑦 values are compared at SG1 and SG2 (Figure 4 (c)) i.e. the critical locations. The strain s at these locations from the linear submodel (EL submodel), nonlinear submodel (NL submodel), nonlinear global model (NL global model) and the measured strains of SG1 ( left hole) and SG2 (right hole) are compared at a 5 seconds interval, as shown in Figure 11 . The 𝜀𝑦 values measured at SG1 are lower than at SG2 , which is attributed to the asymmetry in the actual loading on the structure from left to right . This phenomenon is not captured by the load cell in the test and hence could not be replicated in the global FE model , which is driven with the load cell readings. As a result of this, no such asymmetric behavior is observed between 𝜀𝑦 values at SG1 and SG2 in the NL global FE model . This can be observed from the FE global model 𝜀𝑦 values during the minimum and the maximum loa d plateaus , shown in Figure 11 . The asymmetr ic behavior at the critical locations is captured by both EL and NL HDMS submodels. In the case of HDMS submodels , the measured boundary conditions that drive the HDMS models capture the asymmetry , which is propagated to the estimated response. Although both EL submodel and NL submodel provide a similar response , EL submo del adhering closely to response measu red by the strain gauges, the NL submode l gives a more accurate estimation of the behavior at critical regions a s the location SG2 goes into nonlinear strain rage for the A572 steel at 𝜀𝑦~2215 𝜇strain in the maximum load plateau. It should also be noted the strain gauges give the average of strain values over 1cm length at the cr itical locations where the spatial variation of strains is very high. Whereas, the NL submodel provides strain values averaged over integration points of neighboring elements of 1.27 mm ( 0.05 in ), which can provide a better estimate at the critical locations. A preprint - July 1 4, 2020 14 Figure 11 Comparison of strains from EL submodel, NL submodel, NL global model and at SG1 and SG2 Figure 12 (a) Absolute difference between strain estimated by NL global model and measured 𝜀𝑦 at SG1 and SG2 , (b) Absolute difference between strain measured by SG1 and SG2 and estimated with HDMS NL submodel The average difference between the 𝜀𝑦 measured with strain gauges and estimated and the NL HDMS at the critical locations of SG1 and SG2 is 39.8 𝜇strain and 56.7 𝜇strain, respectively, as shown in Figure 12. We also observe that this difference is higher during the loading and unloading phases of the test compared to the constant loading phases . This is due to the higher floor noise in the DIC measurements that were used to drive the HDMS submodel in the loading and unloa ding phase s compared to the constant A preprint - July 1 4, 2020 15 loading phase. Figure 13 shows the load vs. strain behavior at SG1 and SG2 form the strain gauge measurements and the NL submodel. The strain 𝜀𝑦 at SG2 shows permanent set at the critical location upon loading indicating the onset of nonlinear behavior. However, SG1 does not demonstrate any permanent deformation. NL submodel also successfully demonstrates the nonlinear response, which can be corroborated by the NL submodel response curve in Figure 13 . The widening of load vs. strain curves corresponding to the NL submodel can be attributed to the noise in the DIC measurements. Figure 13 Load Vs. strain at SG1 and SG2 locations from strain gauges and NL submodel 5 Case study 2: Plate with wide slits This case study demonstrates the implementation of HDMS of steel plate with 1.5 mm wide slits to estimate response at the critical regions of the plate , i.e., the slit tips . The following sections present a description of the test to acquire the cut -boundary response, building , and implementation of the HDMS model of the region containing the slit and a discussion on the response estimated by the HDMS method. 5.1 Test Description The steel plate specimen is 355.6 mm ( 14 in ) wide, 711.2 mm ( 28 in ) tall, and 6.35 mm ( 0.25 in) thick. Let the mid-point of the plate ’s front surface facin g DIC cameras be the origin of a 2D coordinate system . Slit1 of 51 mm length is located between coordinate s (63.5, 90) and ( 114.3 , 90) in the top right region and Slit2 of 25 mm length is located between coordinates ( -101.6 , -63.5) and ( -76.2, -63.5) in the bottom left region as shown in Figure 14 (a). The plate is loaded vertically in tension via channel sections attached to one side of the plate at the top and bottom. Figure 14 (b) shows the profile of applied load from 22.24 kN (5 kips) to 222.4 0 kN (50 kips) in increments of 22.24 kN (5 kips) with alternate stages of 22.24 kN (5 kips) linear increments and constant loading for 20 min utes. Two thousand images of the plate under the incremental loading are captured to estimate strain distribution at each of those time instants using the DIC method . For further details o f the test and image acquisition , please refer to Gulgec et al . (2020) 62. A preprint - July 1 4, 2020 16 Figure 14 (a) Test plate geometry, (b) Loading profile for plate with slits , (c) Sample vertical strain distribution from test measured with DIC. A sample of vertical strain ( 𝜀𝑦) distribution on the surface of the plate, processed from one of the images is obtained via DIC at 222.40 kN, is shown in Figure 14 (c). The discontinuit ies in the colormap the slits , where strains are not computed because of the lack of speckled pattern. The color map indicates the stress concentrations near the slit ends in hues of yellow . The o nly Slit2 on the lower left side of the plate is chosen as the critical region to demonstrate hybrid submodeling , as the other slit is too close to the stress concentrations from the channel sections . A 127mm ( 5 in) × 152 mm ( 6 in) rectangular region around the Slit2 is chosen as th e submodel region , with its top-right corner at (12.7mm, 12.7mm) or (0.5, 0.5) , as shown with a gray boundary in Figure 14 (a). This region is chosen as the submodel such that the cut - boundaries do not cross stress concentrations originating from load grips or the slits. One DIC measurement from each of the constant loading phases is selected to estimate the strain fields. A 3D FE submodel of the chosen region in the plate is built in ABAQUS 2017 with 8-node solid linear elements with reduce d integration . Grade 36 steel properties with Young’s modulus 𝐸=2×105 MPa (29000 ksi), and Poisson’s ratio 𝜈=0.3 are used to define the material p roperties of the elastic submodel. The mat erial properties of the nonlinear submodel are defined as elastic -plastic strain harde ning with a yield stress of 248 MPa ( 36 ksi ), according to Gu lgec et al . (2020)62. The submodel is created with differ ent mesh sizes at the curved tips of the Slit2, as shown in Table 2 to compare input file run time, strain at one of the ends of the slit (left) , and the element aspect ratio s. As the remaining parameters did not change significantly, mesh size of 0.3175mm ( 0.0125 in) is used at the ends of the slit , so that the element asp ect ratio is within the advisable limit (10). The mesh de nsity distribution of the submodel is shown in Figure 15 , with denser mesh close to the slit ends where the stress concentrations are expected under loading . A preprint - July 1 4, 2020 17 Table 2 Comparison of strain at critical locations and runtime for the plate slit submodel for different mesh sizes Mesh size at slit ends, (mm (in) ) Run time (sec) 𝜺𝒚 at highest stress location (𝝁𝒔𝒕𝒓𝒂𝒊𝒏) Element thickness/ element width 0.1127 (0.05) 2.40 1265.9 1 2.5 0.0635 (0.025) 2.50 1260. 74 5 0.03175 (0.0125) 2.60 1275.0 1 10 0.0254 (0.01) 2.80 1272. 19 12.5 0.015875 (0.00625) 2.80 1288. 17 20 Figure 15 FE mesh design of submodel for slit 2 5.2 Implementation of HDMS and Results HDMS is performed with elastic and nonlinear submodels in ABAQUS 2017 at each constant load stage. The DIC strain measurements on the four cut-boundaries are interpolated using a cubic spline to each of the element faces on the FE submodel boundaries. The interpolated normal and shear strains on the respective boundaries that is , 𝜀𝑥 and 𝜀𝑥𝑦 on the right and left boundaries and 𝜀𝑦 and 𝜀𝑥𝑦 on the top and bottom boundaries , are converted to stresses (Equation 3) and used to drive the submodels. Figure 16 (a) shows the DIC 𝜀𝑦 distribution and Figure 16 (b) shows the 𝜀𝑦 distribution estimated with the nonlinear submodel at 222.4 0 kN. The DIC 𝜀𝑦 distribution is sparse compared to the refined submodel and the discontinuity near the slit omits the values near the stress concentration locations at the slit ends , so it is difficult to quantify the accuracy at all point s in the submodel . However, f rom observation, both have similar 𝜀𝑦 distribution pattern including high strain in the top right region of the submodel . The pronounced contrast in the 𝜀𝑦 distribution around the slit in in the estimation by the submodel is not observed in the DIC measurement . This is the result of the larger facets used in post-processing DIC measurements to limit the spatial noise while comput ing strain , which did not allow t he capture of sharp changes in st rain near the critical regions. For a better understanding of the performance of HDMS, the strains at eight individual locations around the Slit2 in the submodel region are compared to interpolated DIC measurements that are available closest to these locations . Figure 17 shows the locations of these individual points and comparison of strain from HDMS and DIC strains at each location . The behavior of both linear and nonlinear submodels coincides A preprint - July 1 4, 2020 18 at all stages of loading and the selected 8 locations depicted in Figure 17 as the submodel behavior remains linear at all these locations, even at the maximum load stage. Figure 18 and Figure 19 show 𝜀𝑦 difference and percentage error of estimated 𝜀𝑦values from the DIC measurements at the selected locations for each loading stage . The discrepancy is observed higher for the location s closest to the stress concentrations , such as location 1. The strain estimated at location 4 shows a high discrepancy compared to the DIC values as the DIC post-processing could not sufficiently capture the high contrast near stress concentrations . At the remaining locati ons, the strain difference is within 50 𝜇strain or 20% error at higher load phases . The high % error at initial loading phases is attributed to the low signal to noise ratio at smaller loads. Figure 16 (a) Strain distribution form DIC measurement in the submodel region at 224 kN (b) Strain distribution of nonlinear submodel at 224 kN Figure 17 Comparison of strain measured with DIC and estimated with HDMS submodeling for plate with slits A preprint - July 1 4, 2020 19 Figure 18 Difference between strain m easured with DIC and estimated with HDMS submodeling for plate with slits Figure 19 Percentage error of strain estimated with HDMS submodeling for plate with slits A preprint - July 1 4, 2020 20 6 CONCLUSIONS This study proposes a hybrid data + model -based submodeling (HDMS) method for a refined response estimation near the critical locations in a structure. To imp lement HDMS, a chosen submodel region around the critical locations is driven by sensor measurements interpolated from cut -boundaries on the actual structure to submodel boundary elements. The implementation of the HDMS method is demonstrated through two c ase studies of plate structures in tension with different types of critical locations, i.e., Case study 1: holes connected with a horizontal slit as the critical feature and case study 2: wide slits in the plate, perpendicular to the direction of loading. The boundary strain measurements were obtained with DIC and strain gauges for the implementation of HDMS with force -based submodeling, which is more feasible in real structures. The HDMS method successfully estimated the refined response measurements with a strain error of <= 50 𝜇strain at critical locations in the first case study and at locations closest to the critical locations in the second case study. The error is within 20 % higher loads compared to lower loads due to the low signal to noise ratio of DIC in the lower load cases. The HDMS method could predict the effects of unsymmetrical loading mimicking different strain values at the critical locations on the outer edge of the holes in case study 1. This could not be captured by a global model driven by the load cell readings. The HDMS method could also estimate the nonlinear response of the material at the critical locations mimicking the permanent set experienced by just o ne hole because of the unsymmetrical loading. Hence, the proposed HDMS is a reliable and a versatile, refined full -field response estimation method for fatigue damage -prone critical locations in real structures at a lower cost. Given the geometry of the region containing critical locations and enough sensors on the submodel boundary to capture any boundary stress variations, the HDMS method is easily extendable to response estimation at critical features of even three - dimensional structural components. ACKNOWLEDGMENTS Research funding is partially provided by the National Science Foundation through Grant No. CMMI - 1351537 by Hazard Mitigation and Structural Engineering program, and by a grant from the Commonwealth of Pennsylvania, Department of Community and Economic Development, through the Pennsylvania Infrastructure Technology Alliance (PITA). BV would like to thank N S Gulgec for providing DIC data for the case study of the plate with slits that are used in t his article. BV would like to thank J Saunders and Y Chen for their help with data collection in the case study of the plate with holes. ORCID Bhavana Valeti https://orcid.org/0000 -0003 -2189 -6339