Surrogate models for dynamic load factor and ductility ratio in blast response determination

In this article, the blast response of structures by the Single Degree of Freedom (SDOF) method is revisited. The existing literature in which the Biggs’ chart is used to determine the ductility ratio is examined in detail. The numerical method determines the SDOF response by considering the elastic perfectly plastic behaviour of the structure. The numerical values of the ductility ratio in various plastic Dynamic Load Factors (DLF) in each td/T ratio are collected. Such verifiable data can represent the whole chart in the elastic and plastic regions. The only available implicit formula suitable for sharp pulses is tested. It is shown that substantial inaccuracy exists in the formula. Since the error in the existing formula is up to 100%, an alternative formula is suggested to reduce the error. The new nonlinear surrogate model describes the chart by introducing 3 regions. In each region, an approximate formula is developed. The model computes the ductility ratio with less than 2% error. The newly suggested functions are nonlinear quadratic types and have been developed by using high-order polynomial optimization. In the numerical example, the result of the new surrogate model is commented on in comparison with SDOF and FEM. It is shown that it can be used in the computational design of protective structures without using Biggs’ chart. It is concluded that similar surrogate models can be developed for unsymmetrical blast pulses.


Introduction
Modern technologies for condition monitoring of assets are growing fast.They are based on a new concept known as digital twins [1] of real physical assets.Each component of the asset requires a dynamic model that is included in the digital twin.
Improvement in the twin capability, for fast prediction of the outputs, is an active field of research.The tool is based on expressing the results of the simulation of underlying dynamical systems, via the tailored surrogate models [2].Such surrogate models are under investigation in the field of structural dynamics [3].
The blast-resistive protective structures are included in all offshore platforms.They are faced with possible hydrocarbon explosions.In the design of protective structures, a mixture of analytical, numerical, and graphical methods is used [4].The blended approaches in the design can use an accurate surrogate model.
The initial design of any protective structure is based on SDOF analysis, in which one can calculate the ductility ratio from Biggs' chart [5].The chart is still a tool for designers, in blast wall design [6].It is desirable to replace the chart with a surrogate function, such that it can be used not only for designers but also for the digital twin of the protective structure.
There is an alternative simple method to SDOF which is known as the Pressure-Impulse (P-I) diagram method and is described in [4].The method is graphical, it is appropriate when a linearly decreasing blast pulse is applied [7][8][9].
The P-I method is based on curves that are enveloped in horizontal asymptotes (long pulses) and vertical asymptotes (short pulses).Therefore, they are accurate for short and long pulses only.In order to make the P-I diagram accurate, it can be developed individually for particular structures by using high-fidelity models [10,11].
The P-I model can also be modified to include different material models [12].However, in the blast pulses where pressure rising time is significant (like gas explosions), the corresponding P-I diagram should be modified and redeveloped using the SDOF method [13].
The SDOF method is not graphical and can be used for a variety of blast pulses.It is described in [4,5].SDOF method requires the equivalent mass and stiffness, that can be found via simple methods [14][15][16].
High-fidelity models may also be required to provide the appropriate equivalent mass and stiffness [6,17].For the nontriangular pulse shape and nonlinear material model, the SDOF can also be implemented successfully [18].The validity of the method is also tested, see for example [19].
In this article initially, the elastic perfectly plastic model in SDOF analysis is revisited and the ductility ratio is determined and stored as verifiable data that can accurately express the Biggs' chart.Thereafter it is shown that the existing surrogate function for the chart in [4] is substantially inaccurate and therefore it should be updated.
The new nonlinear quadratic surrogate function is developed by using a high-order polynomial optimization technique.Thereby 3 surrogate functions are found for 3 regions of the chart.In each region, the computed ductility ratio indicates less than a 2% error.
Region III which covers the elastic region, provides a surrogate expression for the elastic dynamic factor and was previously reported by the author in [20].In this article, regions I and II are also developed that provide the surrogate expressions for ductility ratio in plastic regions.
It is discussed that the existing Biggs' chart assumes the blast pulse is a symmetrical isosceles triangular type, and this is a hypothetical assumption.The approach in this article is based on practically accepted pulses with unsymmetrical shapes.It is concluded that similar surrogate functions can be developed for unsymmetrical pulses via the approach in this article.

SDOF Response to an Unsymmetrical Pulse Force
When an unsymmetrical triangular pulse is applied to a mechanical system in Fig. 1 with mass M and the stiffness k the equations of motion in the SDOF approach are: Equations (1) are valid when the system is in elastic status i.e. ( x ≤ x el ).However, when the system faces plastic defor- mation, where the material model is Elastic Perfectly Plastic E-P-P material model, Eqs.(1) will change to: (1) The maximum resistance R m in Eq. ( 2) depends on maximum elastic deflection x el and the stiffness k , via this formula: The plastic dynamic load factor depends on maximum resistance R m and is defined by this equation: The natural period of structure T is given by:  1) and ( 2) can be converted into the following dimensionless form: In the elastic region, i.e. for x < 1 In perfectly plastic region i.e. for x > 1 (2) From the numerical solution of the equations in the elastic region (6) and plastic region (7) the history of the ductility ratio x( ) can be found and the ductility ratio in the Biggs' chart is: The full derivation of Eqs.(6,7)

Revisiting Existing Empirical Formula
The Biggs' chart was first developed by the US Army in [21] but appeared to the public in a book written by Biggs [5].It is known by some authors as design chart [22].It is developed by analogue computers and the results are graphical.In Fig. 2 the chart shows the ductility versus other parameters [4,5].The force pulse for this chart is isosceles triangular type i.e. = 0.5 .It can be reproduced by Eq. ( 8) in this article.
The only empirical formula that can express the chart is given by Newmark [23].It is used as an alternative to pressure impulse diagram [24] and is written as: This formula ( 9) is good for α = 0, but it was also used for any types of pulses for few decades [4,25], to estimate ductility without using Biggs' chart.Herein we have tested Eq. ( 9) to see if the results are matching with (8) and the comparison is shown in Fig. 3.
Figure 3 shows that the approximate formula (9) overestimates the values of the ductility ratios substantially.The accuracy is acceptable only if (DLF) R < 0.5 i.e. for strong blast forces whereas when the loading approaches the elastic limit i.e. (DLF) R = 1 the overestimated error sometimes is above 100% so ( 9) is not appropriate, since it is valid for α = 0 only.

New Surrogate Expressions
By using accurate surrogate models, engineers and researchers in field of blast response, do not need to learn sophisticated tools like digital spreadsheet analysis SBEDS.This tool is developed by the US Army and employs a numerical time-stepping procedure to solve the equation of motion of a mass-spring system.Some predictions of SBEDS are shown in [26].
In order to find an accurate surrogate expression, that can predict ductility ratio, a new polynomial optimization technique is used that is valid for 0 < α < 1.To present the whole chart, a high-order nonlinear quadratic form is required.However, for lowering the order, it is required to divide that chart into 3 regions and find the surrogate function in each region.The method herein is new in this field but used before at [27] to predict the surface of the satellite mirrors for online position control.It has been successful and has been updated recently in [28].The 3 regions are identified by examining the data that is generated from (8) which picks up the maximum from the numerical solution of ( 6) and (7).These 3 regions represent all the data in Fig. 2. Looking at Fig. 2, we set the vertical axis y = DLF R and horizontal axis T , then in the following range (I): The ductility ratio can be approximated with the following polynomial: (10) y = DLF R Fig. 3 Ductility ratios comparison using Eq. ( 8) and (9) However, in this range (II): The ductility ratio can be approximated with the following polynomial: In the region III (elastic) with y ≥ 1 , the polynomial is in terms of d only, such that: The coefficients of the polynomials (11-13) can be found from multiple regression analysis with arbitrary terms based on least square method.The procedure is explained in many references including [27][28][29] and is based on minimising the following error: Then the parameters p a,b in ( 14) can be found from the following system of equations: In (14) the N is the number of data point that is generated by (8).In (15) n is the number of parameters.The Eqs. (11)(12)(13) are rewritten into the following form: For region I the formulae includes 11 terms as follows: (11) Region II 0.75 ≤ (DLF) R ≤ 1: For region II the formulae include 14 terms as follows: For region III (elastic region) the formulae include 18 terms as follows In order to test the expressions ( 16) and ( 17), the Fig. 4 is produced.It shows that the results are comparable with Eq. ( 8).In both regions I and II there is perfect match between two sets of the results and the maximum error does not exceed 2%.It is obvious that, Eqs.(16,17) can replace the Biggs' chart in Fig. 2. Such formulas can be stored in calculators or inside the other computational software to determine the ductility ratio.Ultimately there is no need to refer to Fig. 2, and the whole design procedure can be computerised and inserted into a digital twin for online prediction of the damage in a protective structure.Moreover, formulas are suitable for design engineers and removes the errors that can be produced by visual examination of the chart in Fig. 2.
The author claims that Eq. ( 8) can reproduces the Biggs' chart by using digital computers, whereas the Biggs' chart in Fig. 2, is produced by analogue computers in 1957 [21].Therefore, we assume that Eq. ( 8) provides an accurate result for the Biggs' chart with zero error.

Numerical Example and FEM Verification
A triangular pressure pulse with peak of p max = 1.933 bar and duration of t d = 58 ms is applied to a steel blast wall.The pitch is p = 1.2 meter and the other dimensions are shown in Fig. 5.It is one of the existing profiles of blast wall that is described in [6] The second moment of cross section I = 8.767 × 10 −5 m 4 , the section modulus W pl,y = 4.37 × 10 −4 m 3 , mass per pitch M = 410 kg, Young modulus E = 210 GPa, and yield stress f * y = 400 MPa, the length L = 3 m.According to instruc- tions in [6], the natural period for the blast wall is calculated T = 16.1 ms.The SDOF model for this blast wall results the deflection y max = x el ≅ 58.4 mm and μ is found from Eq. ( 8).This is shown in Table 1.
In order to verify the SDOF result a comparison is done by using FEM technique via ABAQUS modelling [30] for the blast wall in this example.The meshing is shown by a snapshot in Fig. 6.In this model there are 6500 shell type S4R elements, each with nine internal integration points.Moreover, there are substantial FEM outputs, including the Fig. 4 Ductility ratios comparison using Eq. ( 8) and (16,17) Fig. 5 A typical cross section (one pitch) of a blast wall [6] local buckling details in bottom flanges as shown by red colour in Fig. 6.
The only result that can be compared with SDOF analysis is the displacement of the node in middle of the top flange.This has been extracted from output files and is y max ≅ 52.99 mm .Since ductility ratio is not defined the in ABAQUS, the Table 1 is provided to compare the (maximum deflection) in each approach.
The 2 nd row of the Table 1, is found from history of the displacement of the middle of the top flange of the blast wall.In the 3 rd row of the table, the displacement is found by using y max = x el ≅ 62.4 mm and is found from surrogate Eq. ( 16).This shows that Eq. ( 16) provides a conservative result as well.
The structure remains in elastic status if the following inequality holds: The inequality (A-6) can be expanded as follows: Substituting (A-5) into (A-4) and the result into (A-7), also using (A-3) and (A-1) in right side of (A-7) yields to: By numerical simulation we can find the maximum deflection and we can check if (A-6) or (A-8) holds, then we can find if the structure is in elastic or plastic status.
When an unsymmetrical triangular pulse is applied to a mechanical system with mass M and the stiffness k the equations of motion in SDOF approach are: The equations (A-9) are valid when they system is in elastic status i.e. x ≤ x el .However, when x > x el the system facing plastic deformation.For the Elastic Perfectly Plastic E-P-P material model, the equations (A-9) will change to: The maximum resistance R m is defined by equation (A-5).The equations in (A-9) and (A-10) can be changed to: However, M k can be expressed in terms of natural period of structure T as follows: Substituting (A-13) into (A-11) and (A-12) results: Substituting (A-5) into (A-15) then dividing the result by x el results (A-17).Similarly dividing (A-14) by x el results (A-16).
Considering dimensionless parameter, and together with substitution of equations (A-2), (A-4) and (A-5) changes the equations (A-16) and (A-17) into this form: T are introduced which yields to: Considering (A-20) the equations (A-18) will change in elastic region, i.e. for to: Also equations (A-19) will change in perfectly plastic region i.e. for x > 1 to:

Fig. 1
Fig. 1 Mass spring system with unsymmetrical pulse force