Analysis of liquid spring damper for vertical landing reusable launch vehicle with network-based methodology

This paper presents the network-based modeling, validation and analysis of the nonlinear liquid spring damper model under vertical landing conditions of reusable launch vehicle. The impedance function of damper model is derived first. Then, its mechanical and hydraulic networks are newly established based on the hydro-mechanical analogy and network-based analysis. By comparing the networks between the corresponding symmetric and asymmetric structures, the meaning of each branch in the network is elucidated. After that, the validity of the network-based model for the liquid spring damper is confirmed by comparison against the experimentally verified nonlinear model in both frequency and time domain. The force and energy absorption characteristics of the damper model are further decomposed, and, specifically, the influence of the orifice area and orifice length on the attenuation performance is studied. The results show that the network-based model provides predictions consistent with those generated by the nonlinear model. The main discrepancy is attributed to the inaccuracy caused by the equivalent fluid bulk modulus. The network-based analysis indicates that the orifice area mainly influences the damping force in the network, which further affects the loads and efficiency of the damper. The orifice length mainly influences the inertia force in the network, which should be limited to a small value. The proposed novel interpretation of the damper models and responses under impact conditions constitutes a framework suitable for systematic design of typically highly nonlinear landing systems in reusable launch vehicles.


Introduction
Due to the anticipated flight cost reduction and an increase in launch flexibility to be brought by the vertically landing reusable launch vehicle (RLV) technologies, many organizations consider them to be a promising mode of space transportation [1]. For these RLVs, the touchdown event is the final step for their recovery which closely influences the success of the whole mission [2]. The liquid spring dampers are commonly used in RLV landing systems to attenuate the touchdown impact due to their compact structure, light weight, high reliability, and reusability [3]. Consequently, it is necessary to study the buffering characteristics and influence of the key parameters on the resultant forces and energy absorption of liquid spring dampers and other related devices.
Many researchers have conducted theoretical and experimental analysis of liquid-based dampers used in landers or other vehicles. Welsh [4] analyzed the dynamic characteristics of a helicopter hydraulic damper with the gas chamber and spring-assisted valve under low and high frequency excitation using simulations and tests. Raja et al. [5] developed a liquid spring damper with high spring rate for a vehicle suspension system. In that work, a fluid-based numerical model was proposed to predict performance of the system at different conditions. Yue et al. [6] designed a novel landing system for the landers using doublechamber and single-chamber hydraulic dampers in the primary and auxiliary struts, respectively. The dynamic models were established, and their buffering performance was analyzed under the critical landing conditions numerically and experimentally. Choi et al. [7,8] developed the design-focused analysis and control of adaptive magnetorheological hydraulic dampers to enable adaptive shock mitigation in a lightweight helicopter. A mathematical model was established to obtain the time-domain landing responses under impact conditions. Wang et al. [9] proposed an adaptive landing system with liquidbased dampers for an RLV, and its model was also implemented. Both, the influence of the structural flexibility and friction on the landing performance was analyzed. Lei et al. [10] established the dynamic model of a single chamber hydraulic damper used in the RLV landing legs. The design parameters were then optimized to reduce the maximum acceleration and strut force based on the dynamic model. Gan also performed the analysis and optimization design of an oleo-pneumatic landing gear to tune the resulting damping force [11]. The above research mainly focused on building nonlinear dynamic models of the landing suspension systems to analyze their dynamic performance in the time domain. More focused analysis and interpretation of the force generation and energy absorption mechanisms, as well as their links to the device architecture, have been typically omitted.
To enable analysis of the liquid spring damper attenuation characteristics in the frequency and time domains, the network analysis concepts are used in this paper. This approach was originally applied to design the electrical circuits with the required impedance characteristics by constructing and then modifying their network diagrams. In recent years, with the introduction of the inerter component leading to the full analogy between the electrical and other physical domains, including a hydraulic one, the passive network synthesis concepts were made applicable for the design and analysis of liquid-based dampers [12]. Some researchers studied the dampers based on this approach and furthered the idea of the equivalent mechanical networks. For instance, Hu et al. [13] used six selected passive mechanical networks to conduct analysis of the vehicle suspension system by balancing the ride comfort, suspension stroke and tire grip ability with multi-objective optimization method [14,15]. Shen et al. [16] designed a vibration mitigation approach for vehicles based on a selected spring-damping-inerter (ISD) network configuration. The damping characteristics of the network were then analyzed in the frequency domain, and validity of the model was verified through a 1/4 vehicle vibration experiment. Giaralis [17] applied the mechanical network approach to develop the wind-resistant vibration mitigation in buildings. The analytical solution of the dynamic responses was deduced, and vibration performance under the different frequencies and peak wind loads was studied. The accuracy of the network representation was then verified experimentally. Li et al. [18] studied and optimized the aircraft landing performance using a range of candidate shock absorber designs represented by their corresponding mechanical networks.
Some researchers used the domain analogies to establish the equivalent mechanical networks from the underlying hydraulic damper networks aiming to study their predictive potential and basic damping characteristics. For example, Swift et al. [19] built the dynamic model of the liquid damper with a helical tube and then proposed its corresponding hydraulic and equivalent mechanical networks. The validity of these models was verified using the harmonic vibration tests. The energy absorption characteristics and their frequency dependence under the different fluid bulk modulus values were also analyzed. Dario et al. [20] established the mechanical and hydraulic networks of a liquid-based damper design used in buildings and optimized the buffering parameters present in the network. Liu et al. [21] built the dynamic model of a symmetric double-rod hydraulic damper while taking into consideration the liquid elastic, inertial and nonlinear frictional effects. The equivalent mechanical and underlying hydraulic damper networks were then deduced. The device network topology was investigated highlighting the relationship between the damper parameters and characteristics. Traditionally, the verification of the established network was completed through a range of quasi-static and dynamic tests. Liu et al. [22] completed investigations of the mechanical and hydraulic networks for multiple liquid damper designs. In particular, the correlation between the elements from the equivalent mechanical and underlying hydraulic network was elaborated. The network parameters were identified from the dynamic tests under the range of conditions, and predictive potential of the proposed network was experimentally shown by studying the magnitudefrequency and phase-frequency transfer characteristics. It can be concluded from the above research that the key in establishing a reliable network representation is to determine the topology and then coefficients of the equivalent spring, damping and inerter elements. Once the equivalent mechanical network is established, the next typical step in network analysis is optimal vibration damping or absorption system redesign with the prescribed degree of increased system complexity, e.g., in terms of the number of additional components. However, the equivalent network models can also be used in other way. This work introduces an alternative and novel point of view where the equivalent network model is used to focus the investigation on the performance characteristics of either realized or realizable damping systems, their interpretation and analysis of the contributing factors or constituent elements. This is done to improve understanding of the device when in operation. Here, the construction of the network model is an initial enabling step which supports further analysis. Having established this model opens new routes to the analysis of the factors contributing to the produced shock absorbing forces and temporal evolution of the energy absorption in the system. The main advantage of this approach, compared to the classical strategy which uses the overall and integral measures such as the total forces or energy, is the ability to interpret the origins of the responses at their fundamental source level. The proposed process is illustrated through the study of the previously validated nonlinear landing system model with the added challenge arising from the presence of an asymmetric piston. The broad contributions of this work therefore cover the methodology to investigate the shock absorbing systems in highly transient conditions and the development and analysis of the equivalent mechanical network model of the previously unexplored liquid spring system with an asymmetric piston.
The network-based methodology is used to derive the impedance function of the linear liquid spring damper based on the previously developed and experimentally verified nonlinear model [23]. The mechanical and hydraulic network representations of the damper are then established using the hydromechanical analogy, with the physical meaning of each branch in the network explained. The validity of the equivalent network-based model is further confirmed by its comparison with the original nonlinear model in the frequency and time domains. After that, the transient force and energy absorption characteristics of the damper are decomposed in the time domain. This leads to the novel insights and a way of analyzing the spring, damping and inertial forces separately under highly transient conditions. The influence of the orifice area and orifice length on the attenuation performance is studied with the aim to further optimize the damping and inertial factors in the future damping device architectures.
In summary, to rationalize the design and analysis of nonlinear liquid-based damping devices which operate in highly transient conditions, this paper proposes a new methodology to investigate their shock-absorbing characteristics. Within this context, the contributions of this work are summarized as follows: (1) the introduction of the new methodology where a piecewise linear equivalent mechanical network model with the time-invariant topology is used to perform the component-based response force and energy absorption analysis; (2) the development and validation of an equivalent network model of the liquid spring damper with asymmetric piston followed by its functional and performance analysis based on the proposed component-based analysis; (3) the application of the methodology in the analysis of the influence of the primary flow restrictor design parameters on the shock absorbing qualities of the investigated liquid spring damper.
The paper is organized as follows: Sect. 1 presents an overview of the relevant research and the main contributions of this work. Section 2 provides the deduction of the equivalent mechanical and hydraulic network for the liquid spring damper. Section 3 contains the parameter identification and validation of the damper network model and landing impact model. Section 4 investigates the composition of the damper force and absorbed energy under the landing conditions. Section 5 concludes this research with analysis of the selected network-based model parameters.
2 Equivalent liquid spring damper network

Working principle and nonlinear model
The overall configuration of the chosen RLV landing system is shown in Fig. 1. The structure of this retractable four-legged landing system is inspired by the existing design [24]. Each landing leg includes the main strut and auxiliary strut. The main strut consists of the four deployable cylinders and one liquid spring damper for energy absorption. The footpad is mounted at the lower end of the auxiliary strut and is in direct contact with the ground during the touchdown.
Before this, four deployable cylinders are unfolded using the pneumatic control system until they reach the locking position. The liquid spring dampers, located at the lower end of the main struts, are responsible for dissipation of the landing energy. The composition of the considered damper is shown in Fig. 1. The damper mainly includes the outer cylinder, piston, piston rod and sealing components. The piston rod is connected to the deployable pneumatically controlled strut cylinders, and the main (outer) damper cylinder is attached to the auxiliary strut through a spherical joint. When the piston rod compresses, the fluid flows from the lower chamber to the upper chamber, which generates the pressure difference between the two chambers [25,26]. This process converts a significant proportion of the landing energy into heat. Besides this, with increasing compression stroke, the piston rod causes reduction of the total volume of the working chambers and further compresses the fluid. Therefore, the fluid, which is a type of highly compressible dimethyl silicone oil, generates the elastic effects as well. The friction effects due to the sealing components should also be considered during a full working cycle of the damper.
Based on the authors' previous work [23], the threestate lumped parameter damper model, which was verified experimentally, can be represented by the following set of three first-order nonlinear ordinary differential equations: where, respectively, P 1 and P 2 are the pressures of the upper and lower damper chamber; V 1 and V 2 are the volumes of the upper and lower chamber; B 1 and B 2 are the fluid bulk moduli of the upper and lower chamber; A P,1 and A P,2 are the wetted piston areas in Fig. 1 Overall configuration of the RLV landing system the upper and lower chamber; x h and v h are the piston stroke and velocity; q 1 and q 2 are the densities of the upper and lower chamber, q avg is the average density of the two chambers; Q o is the volumetric fluid flow rate through the orifice; F f is the damper friction force, F h is the damper force obtained from the piston rod force analysis. It should be noted that the explicit form of the liquid flow rate model is discussed in detail in Ref. [23].
To enable initial linear analysis, the above model is considered with the constant coefficients relative to the initial working conditions of the chamber pressures and fluid flow rate. Thus, the equivalent mechanical, as well as underlying hydraulic networks, can be formed and the corresponding solution procedure developed. This approach sets the framework for the subsequent force decomposition analysis and parameter influence analysis.

Linear damper model and its equivalent mechanical network
To facilitate the application of the network as well as classical frequency domain analysis methods during the transient landing conditions, the nonlinear model in Eq. (1) is initially simplified by the following assumption: The relations between the flow rates and pressure drops in the damper are assumed to be linear, and the bulk modulus and the volume of each chamber are assumed to be constant in this section. Further, to provide a complete description of this model, a particular pressure-flow model is also assumed here where the liquid flow inertia I h and resistance R h are included in their commonly accepted linear form [23]. This model can be written as follows: where B L , V L , R hL , q 1L , q 2L and q avgL are taken as constants in the linear model; x hL is the damper stroke of the linear model; I h and R h are the inertance and resistance of fluid in the orifice, respectively; Q fL is the volumetric fluid flow rate through the orifice in the linear model. The inertia of the piston is neglected because of its comparatively low value. Using the Laplace transformation with the initial pressure and flow rate set as the atmospheric pressure P atm and zero in the simulation, respectively, the liquid spring damper model can be written as: where P 1L0 and P 2L0 are the reference or initial pressures of the lower and upper chamber, respectively; Q fL0 is the initial flow rate, which is 0 before the damper compression; s is a complex variable. c ðÞ represents the Laplace transformation of the respective states and of the piston motion function. By joint solution of the linear system from Eq. (3) and the damper force model F h from Eq. (1), the following relationship between the overall damper force and the relative velocity between the damper attachment ends (terminals) is derived: where F h0 ¼ A 1 ðP 1L0 À P atm Þ À A 2 ðP 2L0 À P atm Þ, which represents the reference or initial damper force due to the pressurized fluid. c F f is the Laplace transformation of the friction force. From the above equations, the parameters k 1 and k 2 , whose physical units are N/m, represent the damper stiffness characteristics. The parameters b 1 and c 1 , whose physical units are kg and Ns/m, represent the liquid flow inertia and damping characteristics, respectively.
To identify the underlying structure of the equivalent mechanical network, based on Eq. (4), the impedance function [12] of the linear liquid spring damper is expressed further as follows: On the basis of Eq. (5), the corresponding equivalent mechanical network of this liquid spring damper can be obtained as shown in Fig. 2. It can be seen that there are four main branches in the network which constitute the load paths for the forces F h1 , F h2 , F h0 and F f . In the following, these labels denote both the forces and the branches through which they flow. The reference or initial force branch (F h0 ) is connected in parallel with the other branches. The F h1 branch contains a spring element with the spring coefficient k 1 . In the F h2 branch, a damping element c 1 is in parallel with an inerter element b 1 , both of which are in series with a spring element k 2 . Since the friction force is mainly caused by the contact between the piston rod and cylinder seals, the friction branch (F f ) is considered in parallel with all other branches.
The v k2 , v c and v b in Fig. 2a represent the relative terminal velocities of the spring, damper and inerter in the F h2 branch, respectively. The relationship between v k2 , v c , v b and v h can be written as: To clarify the physical significance of the F h1 branch, the coefficient of the spring element k 1 in this branch is simplified with the assumption that the fluid bulk modulus and density of both working chambers are equal (that is B 1L = B 2L = B, q L = q 2L = q). The simplified coefficient k 1 can be written as: Because the expression ðA P;1 À A P;2 Þ=ðV 1L þ V 2L Þ represents the volumetric strain of the liquid spring damper with a unit stroke, based on the definition of the fluid bulk modulus [27], BðA P;1 À A P;2 Þ=ðV 1L þ V 2L Þ is the chamber pressure change with a unit stroke under the quasi-static loading conditions. Therefore, the coefficient k 1 can be interpreted as the stiffness of the quasi-static spring force which is experienced when engaging the asymmetric liquid spring damper quasi-statically. This can also be confirmed by observing that k 1 = 0 when the liquid spring damper has a symmetric double rod cylinder (A P,1 = A P,2-= A), which can be taken as an extreme design condition.
To better interpret the other components within the configuration, the network of the corresponding symmetric damper design is established and shown in Fig. 2b. The coefficients of the elements of the symmetric structure can also be obtained based on Eq. (4), and the assumption of the equivalent bulk modulus and density for both working chambers is also applied here: where A is the wetted piston area of the symmetric double rod cylinder; the subscript s represents the symmetric structure. Based on Eq. (9), k 2s represents the sum of the effective stiffness due to the upper and lower chambers in parallel generated by the fluid during the transient conditions. The damping force F cs and inertial force F bs in the F h2s branch are determined by the flow restriction and inertance characteristics of the orifice, respectively. Through similar reasoning, the spring k 2 in the network of the asymmetric liquid spring damper represents the fluid transient compression, and the damping force F c and inertial force F b are related to the orifice properties.
The hydro-mechanical analogy [28] can now be used to develop the corresponding hydraulic network which links every mechanical component directly to its hydraulic counterpart. It should be noted that the F h0 and F f branches are omitted from this process. This is because the F h0 is generated by the initial chamber pressure and the F f is generated through the mechanical interface interactions, both of which have no relationship with the liquid flows in the hydraulic network. According to the analogy, the damper force (F h ) corresponds to the pressure difference (DP) and the relative terminal velocity (v h ) is equivalent to the total flow rate (Q). Due to the asymmetric configuration of the liquid spring damper, the average area of the upper and lower chamber (A P,1 ? A P,2 )/2 is chosen to establish the relationship between the two domains. The relations between the two sets variables are obtained by introducing the following definitions: Substituting Eqs. (10) and (11) in Eq. (5) and excluding F h0 and F f, , the relationship between Q and DP is obtained as: where DP and Q are the fluid-domain and network- Fig. 3 The hydraulic network of the liquid spring damper specific pressure difference and total flow rate between the chambers, respectively, which will be more closely interpreted in subsequent discussion; C n1 and C n2 are the hydraulic compliances, R nL and I n are the hydraulic resistance and inertance, respectively. The hydraulic network of the liquid spring damper with asymmetric structure can be thus attained, which is shown in Fig. 3a. The relationship between the flow rate, pressure difference of each element in the hydraulic network and the corresponding force, velocity in the mechanical network can be written as: where Q 1 and Q 2 are the flow rates of the two branches in Fig. 3a, respectively; DP c1 and DP c2 are the pressures of C n1 and C n2 , respectively. Similarly, the hydraulic network for the liquid spring damper with symmetric structure can be obtained, which is shown in Fig. 3b. In this figure, the hydraulic compliance C n1 is removed from the hydraulic network due to the equivalent area of both working chambers. Then, the relationship between the hydraulic network elements and equivalent mechanical network elements can be simplified as: When A P,1 = A P,2 (symmetric structure), it can be seen that the pressure and flow rate in Fig. 3b have clear physical meaning based on Eqs. (15) and (16). Q s represents the total flow rate between the upper and lower chambers. It is separated into the transientcompression-related flow (Q 1s ) and orifice-related flow (Q 2s ) components. DP s or DP c2s is the pressure difference between the chambers. DP s (or DP c2s ) consists of DP rs and DP is which represent the pressure difference through the orifice caused by the flow resistance associated with the fluid viscosity and inertance associated with the fluid mass, respectively. When A P,1 6 ¼ A P,2 (asymmetric structure), the additional component of the hydraulic compliance (C n1 ) is present in Fig. 3a due to the asymmetric design. DP and Q for this asymmetric structure do not have the same physical meaning as they have for the symmetric structure, but their hydraulic network is evolved from the symmetric structure. There, it is obtained that the hydraulic resistance is in series with the hydraulic inertance, both of which are in parallel with the hydraulic compliance. Considering only the approximate nature of the linear damper model when representing the nonlinear conditions with the large piston strokes and velocities, a more detailed piecewise linear model and its network is further proposed in this section.

The introduction of the piecewise linear model
The standard relationship between the terminal force and velocity for the spring, damper and inerter elements, after their Laplace transformation, can be written as follows: With reference to Eq. (9), when the network element coefficients change their values during the simulations, the linear force-velocity equation for the spring elements in Eq. (17) is no longer representative due to changing piston displacement, whereas the force-velocity equations for the damping and inerter elements are not affected by this consideration. This conclusion can be further extended to the analysis of the damper network.
If the spring coefficients k 1 , k 2 in Eq. (4) are considered to be constants and c 1 , b 1 become variables in simulation, where c 0 1 and b 0 1 are introduced to represent the variable damping and inertial coefficients, then the following equations are deduced for the general network topology in Fig. 2a: From Eq. (18), the impedance function under this scenario is obtained: It is seen that Eq. (19) and Eq. (5) have the same form, which suggests the possibility of the same underlying network topology. However, if the spring coefficients k 1 and k 2 are considered to be variables k 0 1 and k 0 2 , respectively, then owing to their nonlinear dependency on the total piston displacement, it can be It is therefore important to note that the above analysis is enabled by the linearity assumption of the underlying network element models. However, this assumption is significantly limiting for the liquid spring dampers during the normal RLV landing conditions. Therefore, to extend this approach to the full operational range of a working device, a new linearization strategy is developed next. With such approach, it will become feasible to conduct the energy analysis and research participation of each network branch and its constituent elements during the damper operation. The force and energy distribution relationships for all elements during the different landing phases can then be further obtained and analyzed, which is useful for detailed inspection of the buffering mechanism.
A piecewise linear damper modeling and analysis approach underpinned by the network model with the time-invariant topology is therefore introduced here to complement the original fully nonlinear model. In this new model, the overall response is split into the N p calculation steps. The spring, damper and inertance coefficients remain constant in each step, and they are varied between the steps. The number N p of the operational points can be chosen such that the piecewise linearized model yields responses which are representative of the full nonlinear behavior. To derive the equivalent network model, the parameters which change with the damper stroke and velocity are divided into two categories. The first category includes the chamber volume, fluid bulk modulus, density and viscosity, while the second category contains the orifice flow resistance parameter R h . It can be seen from Eq. (4) that the coefficients k 1 , k 2 , b 1 and c 1 change with the parameters in the first category, while only the coefficient c 1 changes with the parameter R h in the second category. Whereas the variation of b 1 and c 1 parameters do not affect the network configuration, the potential use of the displacement-sensitive coefficients k 1 and k 2 would lead to the network topology different from the one introduced in the previous sections. To retain the network topology consistent with the one based on the linear arguments, and thus enable the proposed decomposition analysis, only the parameter R h is taken as a variable during the linearized simulation, while all parameters in the first category are assumed to have their equivalent and constant values. The chamber volumes undergo changes with the varying damper stroke during the touchdown simulations. The chamber volumes are approximated as constant quantities instead of employing the real chamber volumes at each step. As indicated above, the reason for this approximation is to maintain the identical (invariant) network configuration (topology). This is because the change of the chamber volume would cause the change of k 1 and k 2 based on Eq. (4), which would further result in Eqs. (18) and (19) no longer holding and thus the change of the network topology. This could bring the additional difficulties in the following force and energy decomposition analysis. Three different constant chamber volumes are set to verify this approximation, which corresponds to the damper in its full extension state (Condition 1), half extension state (Condition 2) and fully deployed state (Condition 3), respectively. Because the constant volumes in the piecewise linear model mainly affect the spring coefficients k 1 and k 2 , the spring coefficients for these three conditions are calculated and listed in Table 1.
From Table 1, the value of k 1 experiences little change (1.642 9 10 5 N/m-1.669 9 10 5 N/m) for the three different conditions due to the small damper stroke and this trend has slight influence on the landing responses. For k 2 , even though it changes from 1.698 9 10 9 to 1.    Fig. 5, the damper forces for the cases with the constant chamber volumes (network model) and variable chamber volumes (nonlinear model) under impact conditions display only minor differences. The trends followed by these damper forces are consistent, and the numerical discrepancies between the network and nonlinear models are within acceptable range. This confirms the feasibility of the constant chamber volume assumption. The damper force predicted under condition 1 is the closest to that of the nonlinear model. Consequently, the constant chamber volumes are set to the conditions of the damper's full extension state. Based on this analysis, the equivalent chamber volumes are denoted as: where V equ1 and V equ2 are the equivalent volumes of the upper and lower chamber, respectively; V 01 and V 02 are the initial volumes of the upper and lower chamber, respectively. The maximum damper stroke x hmax is used to calculate the equivalent chamber volume.
In addition, the equivalent fluid bulk moduli B equ for both chambers are assumed to be equal. B equ and x hmax are determined with the help of the following procedure: B equ is first determined by means of a simulation test comparison for the full range of the x hmax values. Then, the relationship between B equ and x hmax is established through the polynomial fitting. Finally, under the specific linearization case, B equ and x hmax are found simultaneously by an iterative calculation which is elaborated in Sect. 2.4.

Mechanical and hydraulic network
This piecewise linear modeling approach is introduced to complement the original fully nonlinear model. The overall response is split into different calculation steps. The piecewise linear model can be described by the following set of differential equations: where index i represents the ith calculation step; P 1i , P 2i and Q fi are the state variables of the piecewise linear model; B equ is the equivalent fluid bulk modulus; P 1i and P 2i are the upper and lower chamber Fig.6 The mechanical and hydraulic network of the piecewise linear model Fig. 7 The piecewise linear model evaluation flow chart pressures in ith step, respectively. Q fi is the orifice flow in ith step. The initial value of the state variables is P 1i0 , P 2i0 and Q fi0 . x hi and x hi Á are the damper stroke in the global coordinate system and velocity, respectively. R hi , which is a constant in the ith step, is defined as follows: where l h is the length of the orifice; u avg is the average viscosity of the two chambers used due to the viscosity discontinuity; A o is the orifice area; C d is the discharge coefficient; Q fi0 is the initial orifice flow for the ith step.
On the basis of Eq. (21), and in analogy with the development presented in the previous sections, the mechanical and hydraulic networks of the piecewise linear model are shown in Fig. 6a and b, respectively. The symbol '' in the figure represents the possible nonlinearity in the damping element.

Procedure for the piecewise linear model solution
By solving the piecewise linear model, the participation of the individual spring, damping and inerter elements in the network on the device response can be further analyzed. The procedure for the piecewise linear damper model evaluation is described in Fig. 7.
Step 1: The initial value of the damper maximum stroke is first estimated as x hmax0 , and from there, the equivalent fluid bulk modulus and volume of the upper and lower chamber are calculated.
Step 2: In the first and the ith calculation step: a. based on Eq. (22), the parameter R hi is first calculated. b. R hi is substituted in Eq. (4) to obtain the damping coefficient c 1vi , and then, the damper network model is substituted in a vehicle landing model (Eq. (26)). The damper stroke x h(i?1)0 and velocity x hðiþ1Þ0 Á at time t i?1 are calculated based on Eq. (26) and used as the initial values for the next step. c. Based on Eq. (21), the orifice flow Q f(i?1)0 and the upper and lower chamber pressures P 1(i?1)0 , P 2(i?1)0 at time t i?1 are obtained. Q f(i?1)0 is used for substitution in Eq. (22) to attain R h(i?1) of the (i ? 1)th step, and all three parameters Q f(i?1)0 , P 1(i?1)0 , P 2(i?1)0 are used as the initials for the calculation in the (i ? 1)th step.
Step 3: If t i \t end (t end is the simulation end time), then i = i ? 1 which is followed by the return to step 2. Otherwise, output the new x hmax and proceed to step 4.
Step 4: If x hmax0 À x hmax j j ! e (e is set to 0.01 mm in this paper), then substitute x hmax0 by x hmax , and return to step 1. Otherwise, the whole process ends.
The outer loop is used to come up with an improved estimate of the maximum stroke under different landing conditions. Thus, the selection of the initial stroke x hmax0 does not influence the solution process and determination of the final value. In addition, note that the parameter choices and calculation of the Fig. 8 Landing responses with the different calculation step lengths damper friction force are same to those introduced in Ref. [23].
It can also be seen from Fig. 8 that the quality of the piecewise linear network-based model depends on the calculation step length. Different step lengths (1 9 10 -4 s, 5 9 10 -4 s and 1 9 10 -3 s) are set to obtain the landing responses, which are used for the determination of the calculation step size. The results, represented by the main and auxiliary force, are shown in Fig. 8.
From Fig. 8, the different calculation step sizes between 1 9 10 -3 s and 1 9 10 -4 s slightly influence the predicted landing responses during the initial touchdown instant and damper compression. The change of the step size causes the change of the damping and inertia coefficients. This results in the change of the main and auxiliary forces from -2.68 KN and -0.3 KN to -1.83 KN and -0.01 KN during the initial touchdown, respectively. Also, the minimum damper force and maximum auxiliary force undergo change from -0.82 and 1.32 KN to -0.72 and 1.35 KN at the maximum compression state. In general, the discrepancies between the network-based model and nonlinear model are reduced with the decrease in the step size from 1 9 10 -3 to 1 9 10 -4 s. To maintain predictive quality, the step size 1 9 10 -4 s is adopted for the following analysis.
For the computations, the network-based models and nonlinear models are implemented in MATLAB [29]. Its ordinary differential equation (ODE) solvers are employed to compute the required transient responses. The damper model parameters are adopted from the authors' previous work [23]. The test structure parameters are directly obtained from the test prototype. The network-based model parameters, such as the spring, damping and inertance coefficients, are derived from the nonlinear model. The fluid bulk modulus is identified via the model-experiment comparison.
The established methodology enables the topologydriven response force decomposition which, in turn, allows targeted analysis of all network branches individually during the highly transient damper operation. The calculated force and energy temporal relationships can be thus obtained and used when investigating the benefits of the obtained buffering characteristics.

Parameter identification and model validation
In this section, initially, the equivalent fluid bulk modulus at different maximum strokes is identified. Then, the frequency domain characteristics of the piecewise linear network-based damper model and the original nonlinear model are compared. After that, the time domain response comparison for these two models under impact conditions is conducted.

Identification of fluid bulk modulus
The fluid bulk modulus is identified under quasi-static conditions, where the damper is compressed and stretched at a low velocity v h . The damper is assumed to be in the quasi-static loading conditions, with the aim to determine its steady-state response. Therefore, the stiffness effects are separated from the inertial and damping effects under such conditions. In this section, the axial force and corresponding mechanical network of the damper under the chosen quasi-static conditions are first discussed. Then, the equivalent fluid bulk modulus under the varying damper stroke is obtained through minimization of the discrepancy between the network-based model and nonlinear model. Finally, the identified equivalent fluid bulk modulus and the corresponding damper stroke data are fitted by an approximating function.
According to Fig. 6a, the force of F h1 branch, which includes the spring element k 1 , solely depends on the damper stroke, while the initial force F h0 and friction force F f remain constant during the quasistatic conditions. Therefore, compared with Fig. 6a, Fig. 9 The mechanical damper network under quasi-static operational conditions the only modification of the damper mechanical network under this condition is the F h2 branch. Assuming a constant damper velocity and applying the final value theorem [21], the steady-state response of F h2 can be described as: Note that Eq. (23) shows that the stiffness and inertance items will not affect the steady-state response of the branch. Then, the total axial force can be denoted as: Based on Eq. (24), the mechanical damper network under quasi-static conditions is derived and shown in Fig. 9a. Since v h is a small value, the force F h2 is relatively small compared with other forces. Neglecting the F h2 branch, the damper network under the stated quasi-static conditions can be further simplified as shown in Fig. 9b. There are three branches in this simplified network which represent the spring force F h1 , initial damper force F h0 and friction force F f , respectively. Based on this simplified damper network, the axial forces under the different stroke intervals are obtained, which are later compared with the nonlinear damper forces using the relative root mean square (R RMS ) measure [30]. The nonlinear damper forces of the four different liquid spring dampers, which belong to the individual landing legs, were compared with experiments previously [23].
The stroke-dependent equivalent fluid bulk modulus is identified next. The identified equivalent fluid bulk moduli and corresponding minimum R RMS values are shown in Table 2. The identified bulk moduli under the varying strokes are close to each other, demonstrating thus a low level of discrepancies in machining as well as the assembly and damper filling procedures. The force R RMS values between the two compared models are below 12.5%. The equivalent bulk modulus of the network model remains constant within each simulation, causing a nearly linear relationship between the axial force and stroke. The bulk modulus of the nonlinear model is mainly influenced by the entrapped air during the initial small stroke intervals and by the high chamber pressures under large stroke intervals causing a distinctly nonlinear force-stroke relationship. This factor causes a varying level of the R RMS discrepancy across the considered stroke intervals.
By fitting the polynomial function across the mean equivalent bulk moduli for all four dampers (B equ ) at their corresponding stroke intervals, the following relationship is obtained:  network model under different stroke intervals. These intervals then correspond to the different identified B equ . In Fig. 10e, the red points represent the original B equ and the blue curve represents Eq. (25). The R RMS value between the original and fitting data in Fig. 10b is 1.34%, which demonstrates a good fitting process. Generally, the network-based model with the stroke-dependent equivalent bulk modulus reflects the damper characteristics under the quasi-static conditions with an acceptable level of discrepancies between this and the reference nonlinear model.

Model validation in the frequency domain
The piecewise linear network-based model is compared here with the nonlinear model in the frequency domain to continue its validation [31]. The comparisons, which use the impedance function magnitude and phase responses of the damper under different excitation frequencies and amplitudes, are shown in Fig. 11 [21]. To further interpret the difference between the nonlinear model and piecewise linear network-based model, a linear network-based model without the friction force is introduced as a reference. It is noted that the equivalent damping coefficient of damper c 1 in the linear network model is taken as the average damping coefficient of the nonlinear model Three comparison indices are introduced to analyze the response magnitudes of the nonlinear and piecewise linear network-based models. These are the maximum and minimum value of the responses (shown in the middle column of Fig. 11), and the amplitude of the underlying first-order harmonics is obtained with the help of the harmonic balance method (HBM) approach (shown in the left column of Fig. 11) [32]. In addition, the phase-frequency curves of the HBM responses for nonlinear and piecewise linear network-based model are shown in the right column.
The following analysis describes the differences in the dynamic responses between the nonlinear model and network-based model across the frequency domain while also further highlighting the different dominant features during the excitation. The frequency response of the linear model is first analyzed and shown in the first and third column of Fig. 11. With the excitation frequency increasing from 0.1 Hz to approximately 1 Hz, the magnitude curves of the linear model first decrease with -20 dB/decade slope, demonstrating that the spring k 1 dominates the Fig. 11 The Bode plot under different excitation amplitudes low frequency dynamic response. Moreover, the phase values in this region show a gradual increase from approximately -90 degrees. This could be caused by the increasingly dominant damping element c 1 and decreasing influence of the spring element k 1 (since the phase for an ideal spring and damper are -90 and 0 degrees, respectively). In addition, for the linear model, it is noted that the phase angle growth toward the smaller phase lag angles becomes more pronounced with the increasing excitation amplitude. With the excitation frequency rising from 1 Hz to approximately 35 Hz, all curves of the linear model, except the d 1 curve under 0.003 m excitation amplitude, show a 0 dB/decade slope in the magnitude and they gradually approach a 0-degree phase angle. This is caused by the c 1 damping element's dominance in the system responses in this frequency region. It is noted that this phenomenon does not appear in the d 1 curve under 0.003 mm excitation amplitude due to the relatively low value of the damping coefficient. Here, with the excitation frequency increasing from approximately 35 Hz, the d 1 curve under 0.003 m excitation amplitude experiences an increase with ? 20 dB/ decade slope in the magnitude and it gradually approaches a ? 90-degree phase angle. This is due to the dominant influence of the inertial effect, since an ideal inerter element introduces both the ? 20 dB/ decade growth in the magnitude and the 90 degrees lead in the phase angle. When the excitation frequency is above approximately 60 Hz, both d 2 and d 3 curves under 0.015 m and 0.03 m excitation amplitudes show a decline with the -20 dB/decade slope, which is dominated by the elastic effect caused by the transient compression of the liquid associated with the spring element k 2 in the network model. This is confirmed by the phase decrease to -90 degrees simultaneously. Additionally, it is noted that with the increase in the excitation amplitude, the magnitude curves tend to shift left-ward to the low frequency regions. This is because the growth of the excitation amplitude at the same frequency would increase the corresponding excitation speeds, which further induces the shift of the regions dominated by the damping element c 1 and spring element k 2 , respectively.
Following the analysis of the linear model, the frequency responses of the nonlinear model and piecewise linear network-based model are studied. The curves of the nonlinear and piecewise linear model agree well with the d 1 curve in the low-frequency range [0.1 Hz, 1 Hz], indicating the influence of the spring element k 1 on the damper performance. However, the phase ranges from -63 degrees under the 0.003 m excitation amplitude to -80 degrees under the 0.03 m excitation amplitude instead of -90 degrees in the d 1 curve, which is mainly influenced by the damper friction force. This influence decreases with the growth of the excitation amplitude due to the lower relative presence or participation of the friction force in the total damper force. In the highfrequency range of the 0.015 m and 0.03 m excitation amplitude, the curves of the nonlinear and piecewise linear model display the similar trend with the d 3 curve, indicating the dominance of the spring element k 2 . In the mid-frequency range of 0.015 m and 0.03 m excitation amplitude and the high-frequency range of the 0.003 m excitation amplitude, the magnitude curves increase with the ? 20 dB/decade slope and the phases gradually increase from the negative values toward a 0-degree value. This is caused by the dominance of the nonlinear damping element c 1v . Based on Eqs. (11) and (22), the damping coefficient c 1v is linear with the velocity, leading to the damping force which is proportional to the square of the velocity across the element. The slope of the magnitude curve reaches ? 20 dB/decade, and the phase angle in this case is 0 degrees due to the absence of any phase lead or lag.
Through the comparison between the nonlinear and piecewise linear network-based models, it is observed that they have good consistency in the low-frequency range, while the discrepancies between these two models become larger over 100 Hz, 50 Hz and 30 Hz for 0.003 m, 0.015 m and 0.03 m excitation amplitudes, respectively. The main reasons for these discrepancies are: 1. The equivalent fluid bulk modulus of the piecewise linear network-based model remains constant within each simulation, whereas the fluid bulk modulus for the nonlinear model varies continually with the chamber pressures; 2. the equivalent fluid bulk modulus is identified from the quasistatic conditions, which would cause the difference under the landing impact condition due to its fast changes during the touchdown.
Overall, with the increasing excitation frequency, the responses of the nonlinear model and piecewise linear network model would undergo changes from the low-stiffness-spring-dominated region under the 0.015 m and 0.03 m excitation amplitude and low-stiffness-spring-friction-dominated region under the 0.003 m excitation amplitude through the damping dominated region and then toward the high-stiffnessspring-dominated region. The responses of the network-based model correlate well with the nonlinear model in the lower frequency ranges, which are [0,100 Hz], [0,50 Hz] and [0,30 Hz] under the 0.003 m, 0.015 m and 0.03 m excitation amplitudes, respectively. The discrepancy between the two models increases in the higher frequency range due to the simplified representation of the equivalent fluid bulk modulus, which is still acceptable for this application.

Model validation in the time domain
The time domain landing impact responses of the nonlinear and piecewise linear network-based models of the liquid spring dampers are further compared. The minimal 3-DOF landing model, which represents a symmetric vertical landing regime, is utilized to study performance of the damper. Its governing equations can be written as [33,34]: where m u is the upper mass, m d is the lower mass, g is the acceleration due to gravity, v c is the velocity of upper main mass in the vertical direction in the global coordinate frame, v yf and v xf are the velocities of the lower (landing leg) mass in the vertical and horizontal directions, respectively, F a is the axial force in the auxiliary strut, F n and F t are the contact forces from the ground in the vertical and horizontal directions, respectively, and h and u are the rake angles of the main and auxiliary strut, respectively. Both the nonlinear and piecewise linear networkbased damper models are coupled with the landing model through the damper force F h , while the landing model is coupled with the damper models through the kinematic inputs, the stroke x h and velocity v h . The calculation scheme is summarized in Fig. 12.
Four landing conditions with the initial velocities of 0.626 m/s, 0.886 m/s, 1.4 m/s and 1.715 m/s are considered to evaluate the landing responses of the nonlinear and piecewise linear network-based damper models. The obtained vehicle acceleration, damper force, auxiliary force and the corresponding R RMS values are shown in Fig. 13.
From Fig. 13a, the vehicle acceleration reaches its peak at the touchdown instant, and then, it changes more gradually. The acceleration peak increases with the initial landing velocity from 0.626 to 1.715 m/s. . This discrepancy in the acceleration peaks comes from the difference in the maximum damper loads, which is shown in Fig. 13b and c. From Fig. 13b and c, the initial touchdown damper force peaks in the network-based model have the distinctly pointed and narrow shape, while the peaks in the nonlinear model have a relatively flatter and wider shape. This behavior can be explained with the help of Fig. 11. It is known that the damper undergoes highfrequency excitation during the touchdown instant and the typical range here is between 130 and 230 Hz. In this frequency range, the response magnitudes of the network-based model are larger than those of the nonlinear model, which can be seen in the middle plot of Fig. 11a. After the initial touchdown instant, the induced responses shift to the lower frequency bands, typically within the range of [3.5 Hz,6 Hz], where the two models have a good agreement and correlation with each other. From Fig. 13d, it can also be seen that the touchdown damper force discrepancies affect the agreement between the auxiliary forces predicted by the two models. After that, the trend and magnitude of  Figure 13e shows the R RMS value, the peak damper force and the stroke discrepancies between the two models. The peak damper force is separated into two parts (Max1 and Max2). They correspond to the maximum value during the touchdown instant and the maximum experienced during the subsequent times, respectively. From Fig. 13e, the differences in the Fig. 13 The summary of the landing impact responses in the time domain initial peak force between two models are the largest, between 25 and 35%, while the differences in the stroke are small, within 1.5%. Other discrepancies are generally less than 9%, demonstrating a good match between the two models.
Despite the discrepancy of 25-35% in the prediction of the initial peak touchdown damper force, the acceptable differences obtained in the subsequent stages (less than 9%) admit the force and energy decomposition analysis with the help of the piecewise linear network-based model performed in the next section.

Force and energy decomposition analysis
The resulting damper response is further analyzed based on the validated piecewise linear network-based model. Utilizing the network-based force decomposition, it is possible to perform the energy analysis and study selectively the participation of the individual network branches during the damper engagement. Specifically, the energy distribution relationships for the spring, damping and inerter elements during the different landing phases can be obtained and assessed. Analysis of the force and energy distribution across the equivalent network model can also be beneficial when interpreting the buffering mechanism and developing the guidance for rational and optimal design of dampers.
According to the network configuration shown in Fig. 6, the damper force can be separated into the spring force F k1 , damping force F c , inerter force F b and friction force F f , where these forces correspond to the spring k 1 , damping c 1v , inerter b 1 and friction element branches, respectively. The force of the element k 2 is in parallel with the sum of the damping force F c and inerter force F b . Since the initial pressure of the damper is assumed to be equal to the atmospheric pressure, the initial damper force F h0 in Eq. (24) is zero based on its definition F h0 ¼ A 1 ðP 1L0 À P atm Þ À A 2 ðP 2L0 À P atm Þ. Therefore, the total damper force can be written as: Based on the damper force decomposition, the equation of the energy distribution, which represents the energy stored or dissipated among the various network elements, can be further written as follows: where E in represents the total work done by the damper-induced force within the specified operational time interval and v h1 and v h2 represent the velocity between the two ends of k 2 element and b, c elements, respectively. The elastic energies stored in the spring element k 1 and k 2 are =2, respectively; the v h1 represent the relative velocity across the ends of the spring element k 2 , the s k1 , s k2 represent the relative displacements across the ends (terminals) of the spring elements; the energy dissipated in the damping the v h2 represent the relative velocity across the ends of the damping and inerter elements; the energy stored in the inerter element is =2 and the energy dissipated due to the friction branch is The landing condition with 0.886 m/s initial velocity is selected for the response analysis. The simulation time is set to capture the first compressionrebound cycle of the damper. Since the direction of the forces in each branch could change during landing, only the absolute values of the forces in each branch are shown for convenience of comparison. The results of the force and energy absorption decomposition are summarized in Fig. 14a and b, respectively.
As shown in Fig. 14, the compression stage of the damper located between 0 s and 0.128 s is followed by the rebound stage. The spring force of k 1 (F k1 ) increases from 0 N at the touchdown instant to 2906.3 N where the damper is in its maximum compression. After that, F k1 gradually decreases with the rebound of the damper. The energy storage by the k 1 component experiences an increase from 0 to 25.2 J initially and then a gradual decline. The damping force due to the c 1v (F c ) branch participates significantly on the total damper force during the initial touchdown but decreases gradually with the increasing compression. During the rebound stage, the direction of F c changes its sign and its participation on the total force stays within the relatively low levels due to the comparatively low rebound velocity. The contribution due to the friction force F f remains unchanged between 0 and 0.05 s, and then, it increases due to the surging chamber pressure caused by the damper compression [19]. The friction force component enters a viscous phase after 0.128 s and changes its direction after approximately 0.144 s. The orifice here is designed as having only 6 mm length, and the inertance coefficient b 1 is small. This causes the inerter forces to be very small when compared with the damper or spring forces. This small participation leads to the blue linelike regions, which correspond to the symbols ''E b '' and ''F b ''. The ratio of the spring coefficient k 1 and k 2 is 1.42 9 10 -3 , indicating thus that the stiffness of k 2 is much larger than k 1 . This condition causes only a small compression of the spring element k 2 , which further leads to its low energy storage participation. In addition, it is obtained from Fig. 14b that after the first compression-rebound cycle, the damping (c 1v ) and friction elements dissipate approximately 87% of the total landing energy, while the spring element k 1 stores approximately 13% of the total landing energy. The amount of energy absorbed by these elements keeps increasing monotonously because they dissipate the mechanical energy, while the energies absorbed by the springs first grow and then decline with the progressing compression-rebound cycle of the damper.

Parameter analysis
Based on the proposed methodology, the influence of various design parameters, e.g., the orifice area and orifice length, on the attenuation performance and landing response is studied in this section. The initial velocity of the landing condition is set as 0.886 m/s. It is noted that the analysis still focuses on the first compression-rebound cycle of the liquid spring damper.

Orifice area
To analyze the influence of the orifice area on the damper force, energy absorption and vehicle landing responses, the orifice diameter is set as 1.6 mm, 3.2 mm and 4.8 mm, which is denoted further as the  landing condition 1, 2 and 3, respectively. Table 3 lists the relationship between the chosen landing response measures and the orifice area. Figure 15 shows the variation of the force and energy components for the different orifice areas. In Table 3, the descending and rebound distances represent the sinking and rebounding heights of the vehicle's main body in the vertical direction after the touchdown, respectively. The attenuation efficiency is denoted as: where E inc represents the absorbed energy of the damper during its compression, F max represents the maximum damper force during compression, and S max represents the maximum compression stroke of the damper. Based on Table 3 and Fig. 15, the damping force F c grows significantly with the decrease in the orifice diameter, especially in the initial landing period. With the decrease in the orifice diameter from 4.8 to 1.6 mm, the maximum damper force first decreases from 3.85 to 3.17 KN and then increases to 4.91 KN. This is because when the orifice diameter is 4.8 mm, the maximum force occurs at the maximum compression stroke which is dominated by the spring force, but when the orifice diameter is 1.6 mm, the maximum force occurs at the touchdown instant which is dominated by the damping force. The maximum compression distance and its corresponding force become smaller in response to the high damping force observed at the 1.6 mm orifice diameter condition. The attenuation efficiency increases from 56.74 to 74.47% and then decreases to 58.13% with the decreasing orifice diameter, showing that the optimal orifice diameter is located within the interval [1.6 mm, 4.8 mm]. Moreover, the maximum stroke is reduced significantly with the decrease in the orifice diameter, causing reduction in the participation of the spring force F k1 and friction force F f . In terms of energy conditions, the total energy absorbed (i.e., stored and dissipated) by the damper at the maximum stroke decreases significantly with reduction in the orifice diameter. This is because the overall gravitational potential energy absorbed by damper is reduced with the vehicle's decreasing descending distance from 89.62 to 38.98 mm. With the decreasing orifice diameter, combined with the significant reduction of Fig. 15 Damper responses with the different orifice areas the rebound distance from 77.84 to 0.72 mm, the energy fully dissipated by the damping element c 1v after completing a compression-rebounding cycle increases significantly. Further, it is noted that the inerter force and the energy absorbed by this element are significantly smaller than those of the damping or spring elements. This is because the inertance coefficient b is small for the cases where the orifice length is limited. Generally, the orifice area mainly influences the damping force F c in the network, and its suboptimal (e.g., too small, or large) value may worsen the induced loading conditions and efficiency of the damper.

Orifice length
To analyze the influence of the orifice length on the damper force, energy absorption and vehicle landing responses, the orifice length is set as 6 mm, 500 mm and 1000 mm, which is denoted here as the landing conditions 1, 2 and 3, respectively. It should be noted that the orifice length can be significantly adjusted by installing additional damper features, such as a helical tube or channel plate, outside the damper cylinder [35,36]. Table 4 lists the relationship between the chosen landing response measures and the orifice length. Figure 16 shows the variation of the force and energy components with the different orifice lengths. For better comparison, the vertical axis range in Fig. 16 is set to be the same under all three conditions. Based on Table 4 and Fig. 16, especially during the initial touchdown instants, the proportion of the inerter force F b gradually increases when the orifice length changes from 6 to 1000 mm. This causes a considerable surge in the maximum damper force from 3.17 to 6.45 KN as well as significant reduction in the attenuation efficiency from 74.47 to 35.15%. From Fig. 16b and c, the direction of F b changes three times during this simulation, which occurs once during the compression and twice during the rebound phase and is caused by the alteration of the local acceleration direction. The damping force F c at the touchdown instant is reduced due to the increase in F b . This interaction results from the parallel arrangement of the damping element c 1v and inerter element b 1 . Additionally, with increasing F b , the duration of the compression-rebound cycle also increases from 0.32 to 0.39 s. In terms of energy absorption, when increasing the orifice length, the total and elastic energies at the maximum stroke change from 41.1 and 25.1 J, respectively, to 34.5 J and 19.2 J, which is influenced by the reduction of the descending distance from 74.54 to 64.86 mm. This also leads to the reduction of the rebound distance.
With the introduced parameter variation, the energy absorption capacity of the inerter element b 1 generally increases, while the energy dissipation of the damping element c 1v and friction is slightly reduced. Therefore, the orifice length mainly influences the inerter force F b in the network, and it should be retained at a small value for the liquid spring damper discussed in this paper.

Conclusion
This research develops a novel methodology for analysis of damping devices applied in the highly transient conditions of the landing impact such as those experienced by the RLV landing systems. It consists of two complementary steps. The first step involves the development of an approximate piecewise linear model, along with its equivalent network realization, of the reference nonlinear model. To obtain a time-invariant network topology for the use in the proposed force decomposition analysis, the chamber volumes are assumed to be constant during the damping operation. This choice introduces an acceptable discrepancy between the two models. In the second step, the network model is used to interpret the device responses and parameter effects with the help of the force decomposition and energy absorption characteristics attributed to the individual network components and branches. The intended benefits of this approach stem from its inherent ability to reveal the role and significance of the individual design parameters. In that respect, this methodology can be used to rationalize the design and analysis of novel dampers used in landing systems. The methodology is developed and illustrated on a specific instance of the liquid spring damper which, owing to its simplicity and robustness, has been studied extensively both experimentally and analytically using detailed nonlinear dynamic models. Based on the network analysis approach and the hydromechanical analogy, the impedance function of the new and simplified linear model of the liquid spring damper is derived and its mechanical and hydraulic networks established. The physical significance of each branch in the network is then demonstrated by comparing the networks with the symmetric and asymmetric structures. After that, a new piecewise linear network-based model with the time-varying damping coefficient and equivalent fluid bulk modulus is proposed.
The magnitude and phase impedance characteristics of the original nonlinear model and network-based model for the liquid spring damper in the frequency domain are compared to validate the new model. It is shown that the damper's behavior undergoes specific changes that are attributed to a sequence of the regions characterized by the low-stiffness-spring, low-stiffness-spring and friction, damping, high-stiffnessspring dominated behavior. The landing impact responses of the two models are then further compared, and a good consistency is noted. A single area associated with relatively large, albeit conservative, discrepancy experienced during the initial touchdown instant is observed and linked with the inaccuracy in the equivalent fluid bulk modulus model. Despite this, the piecewise linear network-based model is deemed to reflect correctly the experimentally verified nonlinear model. The damper force and energy decomposition in the first compression-rebound cycle is then completed with the help of the network model. The effects of the orifice area and orifice length are selected for further studies. It is shown that the orifice area mainly influences the damping force component in the network, and that its suboptimal value may worsen the peak loads and efficiency of the damper. The orifice length is seen to primarily influence the inerter component in the network, and it should be restricted to within a small value for this specific design of the liquid spring damper.