Time delays in numerical modeling of frontal thoracic blast pressure wave responses

In this work a lumped mass mechanical model of a thorax subject to a blast pressure wave is taken into consideration. A thorax spring-dashpot model developed by Lobdell is implemented in numerical modeling of dynamics of the multibody system. The five degrees of freedom mechanical model of a chest adjacent to the elastic backrest is subject to an impulse loading generated by the blast pressure wave released by an explosion. The so-called coupling of the pressure wave to the thorax is reconsidered. With respect to the evident existence of inherent time delays of displacements, the system of coupled bodies is described by a time delay differential equations that are derived from the large-scale systems approach. Numerical solutions present interesting dynamical behavior of the bio-inspired system resulting from inherent time delays and a time of arrival of the blast pressure wave. There is pointed out that the inherent state time delays change dynamical response of the multibody system. Proper time of deployment of the foam-based armor plate reduces relative compression of the thorax, which is to be protected by a bullet-proof waistcoat.


Introduction
Intrinsic delays in states of physical quantities characterize many dynamical systems in physics, material engineering, ballistics, biology and chemistry [1][2][3][4][5]. Natural or artificial control systems have delays from the sensing of a variable and the initiation of appropriate response. Mathematics of systems with time delays pose basic mathematical challenges. They can be described mathematically by delay differential equations, which belong to the class of functional differential equations [6]. The effect of a delay influences dynamical responses of investigated systems, and is strongly visible in behavior of multibody systems. The time delay terms of differential equations produce an infinite number of roots of the characteristic equation, making the corresponding dynamical behavior difficult to analyze. One solves such problems indirectly by applying some approximations, but a limitation in accuracy that leads to the instability of systems can occur [7]. More effective methods based on an analytic approach to obtain the complete solution of systems represented by the delay differential equations based on the concept of Lambert W function was developed in [8].
From the other side, numerical methods can be used to aid any mathematical modeling.
Numerical methods are today popular in biodynamics and, in particular, in dynamical analysis of primary blast lung injuries caused by the pressure wave released by an explosion [9]. Rapid motion of the human's chest wall creates a pressure wave in the lung material [10]. A concept to protect the chest is to use a layer foam material behind a massive armor plate worn over the chest. Coupling of the blast wave to the thorax causes that soft tissues like lungs placed in the contents of the thorax can sustain large stress and strain rate [11]. To investigate the mechanical responses of the internal organs, a complex modeling is required. Homogeneous and linear elasticity material properties are assigned to each part of the model, whereas the human cartilages and bones may have different material properties. Such conditions are taken into consideration by Lobdell's model which has been reconsidered in this work to perform numerical solutions of a large-scale time-delay system.
In this paper the problem of modeling of a multibody biomechanical system of thorax is considered. On the basis of the theory of large-scale continuous-time systems an uncertain model of thorax has been written in the representation allowing to define its parametric uncertainties and complex interactions between its subsystems. The parametric uncertainties make the discontinuous system more difficult to solve, but it is compensated by more accurate numerical solution and evident possibility of inclusion of time delays in the impulse responses.
The fact is that, if we consider an exponential decaying response that takes about 0.2 ms, then the time delay superposed on each body of the system play an important role.
A large-scale dynamical system can be characterized by a large number of state variables, system parametric uncertainties, and a complex interaction between sub-systems [12,13]. In view of reliability and practical implementation, time delays have to be incorporated into the numerical modeling of the large-scale physical systems due to the real transport of mass, propagation of vibrations and computation times.

Variation of the air-blast over-pressure wave
Explosions in air create intense shock waves capable of transferring large transient pressures and impulses to the objects they intercept [14]. The free-field pressure time response from an explosion in air is described by known Friedlander's waveform where p(x, t) is the pressure at a point x and time t, p 0 is the blast over-pressure (maximum amplitude), t i is the wave decaying time and v 0 is the sound speed in air. A blast wave propagates outward from an explosion. It consists of a shock front, which precedes a phase of positive pressure and can be followed by a negative pressure phase, which is not taken into the analysis. For the numerical experiment carried out in this work one assumes that the wave has arrived after the explosion in time t arr at the buffer mass m 1 placed at x = 0 in front of the investigated multibody system (see Fig. 2). Therefore, the simplified time-dependent characteristics of the pressure wave [9] is given It will be more useful in numerical computations to represent the pressure wave as the effective impact force u b (t) per an effective area a of the chest subject to the air shock. Temporal variation of the blast wave that reaches the armor plate is estimated using the following conditions: (C1) The most harmful positive phase of the blast wave profile is assumed [15].  Table 1). (C10) Buffer deployment time for an active mitigation concept depends on σ u (ε) as well. (C11) Maximum incident over-pressure cannot exceed 13 kN, which is determined from the lung damage threshold on Bowen curve [16], and for duration of positive incident over-pressure equal to 1 ms.
If a sensor capable of detecting the electromagnetic emission [17] created at the instant of detonation affords a time delay, t a , between detonation and the arrival of the blast wave then a buffer deploys by using a high-speed actuator, such  as a propellant. The force exerted on the buffer as it deploys must assure that the reaction pressure exerted on the protected chest does not exceed 0.3 MPa. For detonations of high explosives the peak over-pressure can cause injury to the thorax. In the study, the foam deploys with a lower peak overpressure determined by σ u (ε) (see Table 1), and the impact time-response dynamics is assessed for a model problem consisting of 10 kg of a high explosion TNT at 2.2 m standoff. The relevant ConWep [16] computations of peak pressure p 0 over the atmospheric pressure, impulse I and arrival time t arr as a function of range are found in [14]. Finally, we can write Corresponding time history of positive phase of the blast wave is depicted in Fig. 1.

The foam-based armor with buffer plate
Foam materials have the ability to deform at low stress level while absorbing mechanical energy. Foam is used in applications of impact protection, to absorb the kinetic energy of an impact, and reduce the maximum stress on the protected object [9]. The foam material reduces peak acceleration while increasing the duration of the impact. Pressure waves released by explosions cause the socalled primary blast injuries. They significantly affect the air-containing organs of human body [18], and in particular, lungs. A lung injury caused by the pressure wave occurs after rapid motion of the chest wall, which creates a following pressure wave in the lung structure. This is referred to as coupling of the blast wave to the thorax. A concept to protect lungs against primary blast injury is to use a layer of foam material behind a massive armor plate worn over the chest [9].
Capability of energy absorbing is the basic feature of foams. They are deformed and absorb the impact energy [19] while keeping the stress acting on the armor plate loaded with the blast wave. One of the common features of energy absorption foam materials is that there is a discernible plateau in their compression stress-strain curves. It means, that the foam materials can absorb energy by deformation, but keep the stress almost constant [19,20].
The typical shape of the stress-strain curve for solid foams has been presented in [21], and including hysteresis, in [22]. The model proposed by Goog encompasses few parameters dependent on relative density, because Young modulus and plateau stress usually increase and densification strain decreases with increasing foam density.
In this work the assumed foam model is a little modification of Goog's model that is completed with a hysteresis. It consist of three systems in parallel, (see in Fig. 2): (G1) The Maxwell arm, containing a spring (linear stiffness k 16 ) and dashpot (viscosity c 62 ) in series. First component describes the first and second region of foam compression and deformation, if the plateau stress is constant. (G2) Linear spring k p . The second stiffness represents a shape of the plateau. It is integrated to describe the increasing or decreasing part of the plateau. (G3) Nonlinear spring k D . The third stiffness is responsible for densification part. It is given by the formula where γ and n are the model parameters, ε is the strain and h f is a thickness of the foam.
Each element of the spring-dashpot model will produce the following components of reaction force: where k p , γ, c 62 , k 16

are stresses in MPa.
Polymeric open-cell foams exhibit complex nonlinear behavior [22]. As it was mentioned, the stress-strain curve for uniaxial compression shows three well distinguishable regimes (G1-G3). The same three regimes are present during unloading of the foam but the response exhibits a hysteresis loop. For the numerical experiment, loading σ (ε)|ε >0 and unloading σ (ε)|ε ≤0 curves are estimated, respectively, for the model parameters of foam densities ρ l = 50 kg/m 3 and ρ u = 40 kg/m 3 (divided by two) listed in Table 1 (compare with [21]).

Mechanical model of the idealized thorax supported from behind
The problem of time-varying position-and velocitydependent system parameters is reflected in the modeling by many factors, i.e. Maxwell elements. Skeletal deflection of the thorax was determined by the difference between displacements of its anterior and posterior walls. The parameter discontinuity observed in connections between the directly coupled lumped masses of the mechanical idealization appear in four cases: (U3) If a foam material is compressed (in loading state) or relaxed (in unloading state), then two different phenomenological solid foam models, which are nonlinear according to (G1-G3) are assumed (see σ l (ε) and σ u (ε) in Table 1). (U4) If the foam is fully deployed and there is not compressing force acting on it, then its modeling is changed from the full nonlinear viscoelastic foam model (G1-G3) to a reduced one, which is modeled by a spring connection of elasticity k p between the proof mass m 1 and the front wall's mass m 2 of the thorax (see σ d (ε) in Table  1). (U5) If the foam achieves a state of full deployment, then the proof mass m 1 (the armor plate) starts to pull the posterior wall of the chest via the waistcoat. In consequence, the spring-dashpot model is being activated, so k 13 and c 13 become different from zero. It is because a deployable foam-based armor plate is assumed in the experiment to be integrated on an outer side of the bullet-proof waistcoat, which is worn over the chest.
In general, discontinuity sources (U1-U5) in the system's stiffness and damping produce coexisting time-varying uncertainties. Modeling of the investigated system dynamics is quite complicated, but the large-scale system's representation does make it easier. An idealized biomechanical model of a thorax supported from behind has been depicted in Fig. 2. Extended models of a sitting human are analyzed in [23]. The concept of the thorax model bases on Lobdell's approach that was developed by General Motors to study the response of the human thorax in automobile crashes [24]. An application of the model has been presented in [9], where the model was consisted of a configuration of springs and dashpot elements. Injury which results from the pressure wave released by an explosion is referred to as primary blast injury [14]. Primary blast injury most significantly affects the air-containing organs of the body [18].
The Lobdell's model was developed through measuring the thoracic response of a human subject to an impact load. Use of the Lobdell model has been extended by researchers to the field of protection against blast, to predict the thoracic response to a blast wave [25].
In this paper, the model and some associated concepts of limiting injuries caused by air blasts mentioned above are reconsidered in one system and supplemented too. The following ones are incorporated: (a) a support attached to the thorax from behind, (b) basic approximation of air pressure wave given in Friedlander's form [26], (c) a phenomenological model for solid foams introduced in [21], An idealized air-blast pressure wave reaches an effective area of the foam-based armor plate, which is represented by the body of mass m 1 , which is by the assumed foam model attached elastically to the second body of mass m 2 (front wall of the thorax, see Fig. 3).

Formulation of the large-scale problem
Let us take into consideration a class of uncertain continuous-time system composed of N -coupled sub-systems as follows respectively: vectors of system states, control inputs, and system outputs, t j = t − τ j . The dynamical system (6) of i coupled sub-systems is described by the internal behavior time-independent state matrices A n i ×n i i , while the control inputs matrix B n i ×m i i , the system output matrix C l i ×n i i and the control inputs transition matrix D l i ×m i represent connections between the external world and the system, τ j is the time delay of jth coupled system. In current investigations, control inputs do not directly influence the system outputs, therefore matrix D l i ×m i i is zero. For the purpose of solution of the proposed problem there are introduced in Eq. (6) time-dependent matrices ΔA n i ×n i ii ΔA n i ×n i ji and ΔB n i ×m i i that define, respectively: the system's state and control input uncertainties. It allows for a more or less precise inclusion of the system's parameters disturbances (U1-U5) given in a form of known time-dependent function or in a quite different form of a dynamically dependent function on the internal system state (time-or state-varying properties are allowed to be modeled as well). Matrices ΔA n i ×n i ji represent all possibilities of connections between interconnected sub-systems of the entire system. Particular view of the biomechanical model of a thorax subject to frontal blast pressure wave has been shown in Fig.  3. From the left side: effective foam armor's mass m 1 under the pressure impact, front body of mass m 2 in anterior surface of the chest wall, rear body of mass m 3 in posterior surface of the chest wall, m 4 -mass of the support.
Equations (6) can be expanded to the following forṁ impact force of blast pressure wave u b (t) is given by Eq. (3), and a reaction force of the support Note, displacements x 5 and x 6 marked in Fig. 3 of massless points in Maxwell elements are expressed in Eq. (7e) by x 51 and in Eq. (7f) by x 61 , while corresponding velocities are obtained directly from a two-point method for approximating the derivative of displacementṡ Non-zero state-space matrices are as follows: In (7g), a difference x 21 (t) − x 31 (t) in displacements of bodies denoted by m 2 and m 3 will be the observed system output. In Eq. (7) zero matrices are neglected.

Uncertainties and the switching matrices
The problem definition uncovers some new features of the investigated bio-inspired system shown in Fig. 3. The problem of time-varying parameters that introduce to the model some uncertainties have been solved numerically by definition of multivalued matrices switched in accordance to cases (U1-U5). Stiffness of chest cave interior with organs increases at condition (U1) about two times to appropriately approximate the real chest's compression. It obviously means, that stiffness of the rheological coupling increases discontinuously with regard to a greater than d compression of the thorax x r , and that damping ability of the coupling varies in time as the thorax is subject to a suitable compression or relaxation. Such discontinuity in the subsystem's stiffness and damping produces parameter uncertainties dependent on state variable.

Uncertainties (U1-U2)
One can encounter in Eqs. (7b) and (7c) four ΔA i j (t) parameter uncertainties of the analyzed dynamical system: where i, j = 2, 3. Distribution of entries in ΔA i j (t) does not change over their all possibilities. Therefore, the unknown time-varying real-valued matrix F(t) is assumed to be defined in the same way, but F T (t)F(t) ≤ I must hold for t ∈ R + , I is the identity matrix. The introduced decomposition DFE includes some known constant real-valued matrices D i and E i j of appropriate dimensions that need to be estimated to use them, for instance, in a solution of LMI problems [12,27].
Dependently on x r (t) and v r (t) the following forms of the uncertainty matrix ΔA 22 (t) = ΔA (k) 22 in Eq. (7b) are possible: where the remaining ΔA i j (t), for i, j = 2, 3 are to be defined in a similar way, defines rheological properties of the biomechanical system. Switching conditions s i will select only one of k possibilities ΔA (k) i j for k = 1 . . . 4 dependently on values of pairs (x r (t), v r (t)) creating the discontinuous time history of uncertainties ΔA i j (t) of the state matrices of coupled two degrees-of-freedom neighboring subsystems.
To find better description of the existing switching nature it is now required to choice D i and E i j matrices of a decomposition. For instance, an exemplary decomposition of ΔA (k) 22 , for k = 3 could be made accordingly to the scheme presented in [28].
Perturbation matrix F(t) will depend on k cases that have been delivered in the case statement (10). Therefore, with regard to F T (t)F(t) ≤ I and using Eq. (9) the following formula reads where i, j = 2, 3, σ(i, j) is given in Eq. (9b), and F(t) will switch between The case statement (12) captures switching properties of ΔA i j (t) described in comments to (10). As it was expected, F(t) is defined in the same way for four uncertainties of the model, matrices D i and E i j are constant and their entries will depend on the sub-system that they are related to.
To check correctness of the previously made estimations let the following parameters be assigned for the uncertainties:

Uncertainties (U3-U4)
While a foam material is compressed (in loading state) or relaxed (in unloading state), two different phenomenologi-cal solid foam models σ l (ε) and σ u (ε), which are nonlinear according to (G1-G3) are assumed (see in Table 1). If the foam is fully deployed and there is not compressing force acting on it, then its modeling is changed from full nonlinear viscoelastic foam model (G1-G3) to a reduced one σ d (ε), which states a spring connection of elasticity k p between the proof mass m 1 and the front wall's mass m 2 of the chest (see in Table 1).

Uncertainty (U5)
While the foam achieves a state of full deployment, the proof mass m 1 (the armor plate) starts to pull the posterior part of the thorax via the waistcoat. The spring-dashpot model is being activated, k 13 and c 13 become different from zero. It is because in the experiment a deployable foam-based armor plate is integrated on an outer side of the bullet-proof waistcoat, which is worn over the chest.
Numerical integration of 12 differential equations of first order is a bit conditioned, but it better approximates complex dynamics of the biomechanical model of the chest subject to an impulsive loading. It is equipped with additional bodies like armor plate and backrest of a seat which have an influence on behavior of the investigated multibody system.

Numerical experiments
Two sources of variability are applied in the numerical simulation to investigate responses of the analyzed multibody system: (1) Various arrival times of impact force u b (t) at foam-based armor plate. (2) Various inherent time delays τ i of each subsystem of the thorax.
Stress F ch in the thorax's model is calculated according to the formula where k 23 ,x 2 ,x 3 are time-dependent. Deployment of the foam initiates all the numerical experiments presented below.

Influence of blast wave time arrivals
Figures 4 and 5 exhibit significant differences in the system's dynamical behavior. The time histories were computed for a few time delays of arrival (ranging from 2.5 ms to 0) of the air-blast pressure wave that reaches the foam-based armor plate. It is seen that for t a equal about 0.5 ms the relative deformation is about 4 cm. If the foam's deployment delay is greater, then the strain in the thorax rises comparing with the non-delayed counterpart. The results prove that even small time delay affects the dynamics, but in this particular case, both the deformation (Fig. 7) and maximum stress (Fig. 8) decrease.

Influence of inherent time delays
Figures 6 and 7 exhibit significant differences in behavior of the dynamical system.
The time histories were computed at very small time delays (ranging from 5h to h = 2 × 10 −6 ) and compared with the non-delayed counterpart. The results prove that even small time delay affects the dynamics, but in this particular case, both the deformation (Fig. 7) and maximum stress (Fig. 8) decrease, as well as a side effect in a form of small amplitude vibrations of higher frequency appears. Figure 8 presents a comparison of stress-strain characteristics of the chest's model deformation estimated for a non-delayed (for τ 6 ) and the delayed displacements of bodies m i , i = 1, 2, 3. One could observe, that even small time delays of about a few time steps of numerical integration play significant role in the blast pressure wave response. Inherent time delays in modeling of the chest are important. Condition (C11) will be satisfied if the time delay τ = 5h (see Fig. 8).

Animation of motion
The six degrees of freedom dynamical system has been numerically solved using dedicated procedures coded in Python. Numpy and matplotlib libraries served as the tools for manipulation on multidimensional tables, plotting of frames and saving them in png files. A set of 1000 frames has been generated and joined to make an animation. Exemplary frame is shown in Fig. 9.

Conclusions
Importance of time delays in numerical solution of the analyzed multibody system is confirmed. Modeling of response Fig. 9 A frame taken from animation of motion of the investigated multibody system of the thorax when energized by air-blast over-pressure has gained new quality. There were mentioned in the literature overview some attempts regarded to the study, but significance of time delays have not been sufficiently emphasized. Therefore, this work sheds new light on mathematical description and modeling of the phenomenon. The system under investigation has received a new useful representation by application of the large-scale systems approach. It allowed to include many uncertainties of parameters and time-delay dependencies of state variables making the modeling more flexible and ready for future improvements.
Interesting dynamical behavior of the bio-inspired system is solved numerically. There is pointed out that the inherent state time delays change dynamical response of the multibody system. Proper time of deployment (initiated about 0.6 ms before an impact of the blast wave) of the foam-based armor plate reduces (at some conditions of the experiment) relative compression of the thorax.