Theory manual

Newmark time integration

In Newmark time integration (Newmark, 1959), the acceleration and velocity at are written in terms of the displacement (), velocity () and acceleration () at time and the displacement at .

(2)

In the above equations, and are Newmark time integration parameters. Substituting the above two equations into the equation of motion will result in a linear system of equations () from which can be estimated.

For and , the Newmark time integration scheme is the same as the trapezoidal rule. The trapezoidal rule is an unconditionally stable integration scheme, i.e., the solution does not diverge for any choice of , and the solution obtained from this scheme is second order accurate. One disadvantage with using trapezoidal rule is the absence of numerical damping to damp out any high frequency numerical noise that is generated due to the discretization of the equation of motion in time.

The Newmark time integration scheme is unconditionally stable for and . For , high frequency oscillations are damped out, but the solution accuracy decreases to first order.

Rayleigh damping

Rayleigh damping is the most common form of classical damping used in modeling structural dynamic problems. The more generalized form of classical damping, Caughey Damping (Caughey, 1960), is currently not implemented in MASTODON. Rayleigh damping is a specific form of Caughey damping that uses only the first two terms of the series. In this method, the viscous damping is proportional to the inertial contribution and contribution from the stiffness. This implies that in the matrix form of the governing equation, the damping matrix () is assumed to be a linear combination of the mass () and stiffness () matrices, i.e., . Here, and are the mass and stiffness dependent Rayleigh damping parameters, respectively.

The equation of motion (in the matrix form) in the presence of Rayleigh damping becomes:

The same equation of motion at any point in space and time (in the non-matrix form) is given by:

The degree of damping in the system depends on the coefficients and as follows: (3)

where, is the damping ratio of the system as a function of frequency . The damping ratio as a function of frequency for and is presented in Figure 1.

Figure 1: Damping ratio as a function of frequency.

Simulation of a constant damping ratio

For the constant damping ratio scenario, the aim is to find and such that the is close to the target damping ratio , which is a constant value, between the frequency range . This can be achieved by minimizing the difference between and for all the frequencies between and , i.e., if

Then, and results in two equations that are linear in and . Solving these two linear equations simultaneously gives:

Rayleigh damping for soils

Small-strain material damping of soils is independent of loading frequency in frequency band of 0.01 Hz - 10 Hz (Menq (2003), Shibuya et al. (1995),Presti et al. (1997), and Marmureanu et al. (2000)). The two mode Rayleigh damping is frequency dependent and can only achieve the specified damping at two frequencies while underestimating within and overestimating outside of these frequencies. The parameters and for a given damping ratio can be calculated as follows:

In case of two mode Rayleigh damping, Kwok et al. (2007) suggests to use natural frequency and five times of it for the soil column of interest. In addition, selecting first mode frequency of soil column and higher frequency that corresponds to predominant period of the input ground motion is a common practice.

Heterogeneities of the wave travel path may introduce scattering effect which leads to frequency dependent damping (Campbell (2009)). This type of damping is of the form (Withers et al. (2015)):

(4)

where, is the frequency above which the damping ratio starts to deviate from the constant target value of , and is the exponent which lies between 0 and 1. Minimizing the difference between Eq. (4) and Eq. (3) with respect to and for all frequencies between and gives:

where, the functions and are given by:

Also, for soils is inversely proportional to the shear wave velocity () and a commonly used expression for of soil is:

where, is in m/s.

Frequency-independent damping

As seen in the previous sub-section, the damping ratio using Rayleigh damping varies with frequency. Although the parameters and can be tuned to arrive at a constant damping ratio for a short frequency range, as the frequency range increases, the damping ratio no longer remains constant. For scenarios like these, where a constant damping ratio is required over a large frequency range, frequency independent damping formulations work better. This formulations is under consideration for adding to MASTODON.

Soil layers and meshing

Small strain properties (shear wave velocity, small strain modulus etc.) as well as mobilized shear strength of soils change with depth. Thus, in numerical models, soil profile (layers) is constructed to incorporate the depth dependent properties. The ground surface as well as layers that define the soil domain can be horizontal or non-horizontal. For the horizontal ground surface and layering scenario, the location of the interfaces can be provided as input and MASTODON will use that information to generate a set of soil layers, each with a unique identification number. These layer ids are later used to assign material properties to each soil layer. The same procedure can also be used for non-horizontal but planar soil layers by specifying the normal to the plane and the interface locations measured along the normal direction.

For scenarios where the soil layers are non-horizontal and non-planar, images (.jpg, .png, etc.) of the soil profile can be provided as input. The different soil layers are distinguished from the image by reading either the red, green or blue color value (as per user's directions) at each pixel. Gray scale images in which the red, green and blue values are all the same also work well for this purpose. For creating 3D soil layers, multiple 2D images with soil profiles at different 2-D cross-sections of the soil domain can be provided as input.

Once the soil layers have been distinguished, it is necessary to ensure that the different soil layers are meshed such that they can accurately transmit waves of the required frequency. The optimum element size for each soil layer depends on the type of element used for meshing, cut-off frequency (f) of the wave and the shear wave velocity () of the soil layer. A minimum of 10 points is required per wavelength of the wave to accurately represent the wave in space (Coleman et al., 2016). The minimum wavelength () is calculated as:

If linear elements such as QUAD4 or HEX8 are used, then the optimum mesh size is . If quadratic elements such as QUAD9 or HEX27 are used, then the optimum mesh size is . Using the minimum element size information, MASTODON refines the mesh such that the element size criterion is met and at the same time the layers separations are visible. An example of this meshing scheme is presented in Figure 2 where a 2D soil domain is divided into 3 soil layers and these soil layers are meshed such that the element size criterion is satisfied. A denser mesh is created at the interface between different soil layers.

Figure 2: Auto-generated mesh for a soil domain with three non-horizontal non-planar soil layers.

Linear elastic constitutive model

In scenarios where the material exhibits a linear relation between stress and strain, and does not retain any residual strain after unloading, is called a linear elastic material. In linear elasticity, the stress tensor () is calculated as , where is the elasticity tensor, and is the strain tensor. This material model is currently used for numerically modeling the behavior of concrete and other materials used for designing a structure in MASTODON.

Nonlinear hysteretic constitutive model for soils (I-soil)

I-soil (Numanoglu (2017)) is a three dimensional, physically motivated, piecewise linearized nonlinear hysteretic material model for soils. The model can be represented by shear type parallel-series distributed nested components (springs and sliders) in one dimensional shear stress space and its framework is analogous to the distributed element modeling concept developed by Iwan (1967). The model behavior is obtained by superimposing the stress-strain response of nested components. Three dimensional generalization follows Chiang and Beck (1994) and uses von Mises (independent of effective mean stress) and/or Drucker-Prager (effective mean stress dependent) type shear yield surfaces depending on user's choice. The yield surfaces are invariant in the stress space Figure 3 . Thus, the model does not require kinematic hardening rule to model un/reloading stress-strain response and preserves mathematical simplicity.

Figure 3: Invariant yield surfaces of the individual elastic-perfectly curves (after Chiang and Beck, 1994).

The current version of I-soil implemented in MASTODON utilizes Masing type un/reloading behavior and is analogous to MAT79 (LS-DYNA) material model but does not exhibit numerical instability observed in MAT79 (Numanoglu et al. (2017)). Masing type un/reloading is inherently achieved by the model because upon un/reloading the yielded nested components regain stiffness and strength. The cyclic response obtained from current version of the model is presented in Figure 4. Reduction factor type modification on un/reloading behavior (Phillips and Hashash (2009); Numanoglu et al. (2017)) is an ongoing study within MASTODON framework.

Figure 4: I-soil model details: (a) 1D representation by springs; (b) example monotonic and cyclic behavior of four nested component model (reprinted from Baltaji et al., 2017).

Main input for the current version of I-soil in MASTODON is a backbone curve at a given reference pressure. MASTODON provides variety of options to built backbone curve for a given soil type and reference pressure using following methods :

1. User-defined backbone curve (soil\_type = 'user_defined'): The backbone curve can be provided in a .csv file where the first column is shear strain points and the second column is shear stress points. The number of nested components that will be generated from this backbone curve depends on the number of discrete shear strain - shear stress pairs defined in the .csv file. When layered soil profile is present, .csv file for each reference pressure can be provided to the corresponding elements in the mesh.

2. Darendeli backbone curve (soil\_type = 'darandeli'): The backbone curve can be auto-generated based on empirical relations obtained from laboratory tests. Darendeli (2001) presents a functional form for normalized modulus reduction curves obtained from resonant column - torsional shear test for variety of soils. MASTODON utilizes this study and auto-generates the backbone shear stress - strain curves. The inputs for this option are (1) small strain shear modulus, (2) bulk modulus, (3) plasticity index, (4) overconsolidation ratio, (5) reference effective mean stress () at which the backbone is constructed, and (6) number of shear stress - shear strain points preferred by user to construct piecewise linear backbone curve. All the other parameters except the number of shear stress - shear strain points can be provided as a vector for each soil layer.

   Darendeli (2001) study extrapolates the normalized modulus reduction curves after 0.1% shear strains. This extrapolation causes significant over/under estimation of the shear strength implied at large strains for different type of soils at different reference effective mean stresses (Hashash et al., 2010). Thus user should be cautious about implied shear strength when utilizing this option.

3. General Quadratic/Hyperbolic (GQ/H) backbone curve (soil\_type = 'gqh'): Darendeli (2001) study presented in previous item constructs the shear stress - shear strain curves based on experimentally obtained data. At small strains, the data is obtained using resonant column test, and towards the moderate shear strain levels, torsional shear test results are used. Large strain data are extrapolation of the small to medium shear strain data. This extrapolation underestimates or overestimates the shear strength mobilized at large shear strains. Therefore, implied shear strength correction is necessary to account for the correct shear strength at large strains. GQ/H constitutive model proposed by Groholski et al. (2016) has a unique curve fitting scheme embedded into the constitutive model that accounts for mobilized strength at large shear strains by controlling the shear strength. This model uses taumax, theta\_1 through 5, small strain shear modulus, bulk modulus and number of shear stress - shear strain points preferred by user to construct piecewise linear backbone curve. The parameter taumax is the maximum shear strength that can be mobilized by the soil at large strains. The parameters theta\_1 through 5 are the curve fitting parameters and can be obtained using DEEPSOIL (Hashash et al., 2016). Other than the number of points, all the other parameters can be given as a vector for the different soil layers. The number of points, which determines the number of elastic-perfectly plastic curves to be generated, is constant for all soil layers.

Once the backbone curve is provided to I-soil, the model determines the properties for nested components presented in Figure 4. The stress integration for each nested component follows classic elastic predictor - plastic corrector type radial return algorithm (Simo and Hughes (2006)) and model stress is obtained by summing the stresses from each nested component:

where is the total shear stress, is the shear modulus of the nested component, is the shear strain, is the yield stress of the nested component, and represents the number of components that have not yet yielded out of the total nested components.

Small strain shear modulus can be varied with effective mean stress via:

where, is the initial shear modulus, is the tension pressure cut off and is mean effective stress dependency parameter obtained from experiments. The shear modulus reduces to zero for any mean effective stress lower than to model the failure of soil in tension. Note that the mean effective stress is positive for compressive loading. Thus, should be inputted as negative.

Yield criteria of the material can also be varied with effective mean stress dependent behavior as:

where, , and are parameters that define how the yield stress varies with pressure.

Lead-Rubber Isolator

Lead-rubber (LR) seismic isolators are composed of alternating layers of natural rubber and steel shims with steel flange plates at top and bottom, and a central cylindrical lead core. Continuum modeling of LR bearings is computationally expensive and impractical for applications involving large structures with tens to hundreds of isolators. Hence, the two-node model developed by Kumar et al. (2015) is adopted here. The isolator model has six degrees of freedom (3 translations and 3 rotations) at each node and is capable of simulating nonlinear behavior in both axial and shear directions. The physical model of the two-node isolator element is shown in Figure 5 below. For a detailed description on the implementation of numerical model, refer to (Kumar et al., 2015).

Figure 5: Discrete, two-noded model of an LR bearing: (a) degrees of freedom, and (b) discrete spring representation (Kumar et al., 2014).

This model captures the following physical behaviors:

• Axial

  • Buckling in compression

  • Buckling load variation with shear deformation

  • Cavitation and post-cavitation behavior in tension

• Shear

  • Viscoelastic behavior of rubber

  • Hysteretic behavior of the lead core

  • Strength degradation due to the heating of the lead core

• Coupling in axial and shear

  • Shear stiffness variation with axial load

  • Axial stiffness variation with shear deformation

Axial (local X direction)

Compression

The composite action of the rubber and the steel shims result in high axial stiffness of a bearing. The expression for the critical buckling load in compression is derived from the two-spring model proposed by Koh and Kelly (1987) as shown in Eq. (5).

(5)

where,

Eq. (5) provides the critical buckling load at zero shear displacement. Experimental investigations have shown that the critical buckling load is a function of the shear deformation. Warn and Whittaker (2006) proposed a simplified expression for the critical buckling load as a function of the overlap bonded rubber area. The critical buckling load is updated at every analysis step based on the current shear deformation in the bearing.

(6)

where,
is the rotational modulus
is the height of isolator excluding the end plates
is the total rubber thickness
is the shear modulus of the rubber
is the initial bonded rubber area
is the current value of the critical buckling load
is the critical buckling load at zero shear displacement
is the current overlap area of bonded rubber

Once the bearing buckles, the axial stiffness is reduced to zero. However, to avoid numerical problems, a very small value of axial stiffness is assigned to the post-buckling region.

Tension

Beyond a certain tensile stress, rubber undergoes cavitation, which is an expansion of existing voids or defects in the material. Cavitation results in an inelastic deformation due to the onset of permanent damage. The LR isolator exhibits linear elastic behavior, till the point of cavitation (,), where, and are the deformation and force, respectively, at the onset of cavitation. Due to the breakage of rubber-filler bonds, the post-cavitation stiffness is reduced. Constantinou et al. (2007) suggested an expression for post-cavitation stiffness as:

(7)

Integrating the Eq. (7), the expression for the post-cavitation tensile force is obtained as:

(8)

where,
is the initial cavitation force
is a cavitation parameter
is the initial cavitation deformation
is the tensile deformation
is Young's modulus of the rubber

When the bearing is loaded beyond the point of cavitation and unloaded, it does not unload with the initial elastic stiffness. The unloading stiffness is smaller than the elastic stiffness because of formation of cavities in the rubber. The cavitation strength is a function of a damage parameter which varies from zero to one. At maximum damage, the tensile strength of rubber layers is not completely zero. Therefore an upper bound of 0.75 value is assumed for the damage parameter.

(9)

where,
is the damage constant
is the current value of the damage parameter
is the upper bound on the damage parameter

For (namely, at the cavitation point), is equal to 0 implying no damage. As the tensile deformation increases, the damage parameter converges to . The reduced cavitation strength as a function of damage parameter is calculated as

(10)

Figure 6 shows the assumed response of a LR isolator to cyclic axial loading (both compression and tension).

Figure 6: Axial response of an elastomeric isolator (Kumar et al., 2014).

where,
is the axial stiffness of the isolator
is the current shear deformation
is the radius of gyration of the isolator

Shear behavior (local Y and Z directions)

The smoothed hysteric model proposed by Park et al. (1986) and extended by Nagarajaiah et al. (1989) is used to model the behavior of a lead-rubber bearing in shear. The shear resistance of the LR bearing is a combination of the viscoelastic behavior of the rubber layers and the hysteretic behavior of the lead core, as shown in Figure 7. The behavior of rubber is characterized by elastic stiffness, viscous damping terms and the nonlinear hysteretic contribution of the lead core is modeled using a Bouc-Wen formulation.

Figure 7: Shear behavior of the LR isolator (Kumar et al., 2015).

(11)

where,
subscripts and correspond to shear 1 (local Y) and shear 2 (local Z) directions respectively, as shown in Figure 5
is the viscous damping co-efficient of the rubber</br> is the shear stiffness of the rubber material </br> is the dynamic yield strength of the lead core </br> is the cross sectional area of the lead core

The Bouc-Wen formulation introduces two hysteresis evolution parameters, and , which are functions of the shear displacements and , in the local Y and Z directions respectively.

(12)

where,

(13)

Eq. (12) is iteratively solved using the Newton-Raphson method to calculate the values of and at every analysis step in MASTODON.

Friction Pendulum^TM Isolator

Friction Pendulum^TM (FP) bearings are used to seismically isolate structures, including critical facilities such as nuclear power plants. The single concave FP bearing consists of a spherical sliding surface, a slider coated with PTFE-type composite material, and a housing plate. Figure 15 and Figure 16 shows components and part sectional view of a single concave Friction Pendulum^TM bearing, respectively.

Figure 15: Components of a single Friction Pendulum^TM bearing ((EPS), 2012)

Figure 16: Sectional view of a single concave Friction Pendulum^TM bearing (Kumar et al., 2015b)

Continuum modeling of FP bearings is computationally expensive and impractical for applications involving structures with tens to hundreds of bearings. Hence a two-node model developed by Kumar et al. (2015b) is adopted here. The bottom node represents the concave sliding surface, and the top node represents the articulated slider. The two-node model has six degrees of freedom (3 translations and 3 rotations), at each node and is capable of simulating nonlinear hysteric behavior in the FP bearing. The physical model and spring representation of the two-node FP isolator element is similar to that of the lead-rubber isolator as shown in fig.

Force-displacement relationship

Friction Pendulum^TM bearings dissipate the energy induced during an earthquake through nonlinear hysteretic behavior in shear. The lateral force-displacement relationship of a single concave FP bearing is characterized by two parameters: yield strength and post elastic stiffness, and can be represented by the bi-linear curve of Figure 8.

Figure 8: Lateral force-displacement relationship of a single concave FP bearing

The yield strength and the post-elastic stiffness can be calculated as,

(17)

(18)

where,
is the yield strength
is the post elastic stiffness
is the instantaneous axial load, and
is the effective radius of curvature of the sliding surface.

Parameters and depends on the coefficient of sliding friction and the axial load on the bearing. For detailed information on the numerical implementation of the nonlinear hysteretic behavior of the FP bearing in shear, refer to Kumar et al. (2015) (page:106-107). Additionally, the FP bearing exhibits very large stiffness in axial, torsional and rotational directions and is therefore modeled as a linear elastic material with large stiffness in these degrees-of-freedom.

Coefficient of sliding friction

The lateral force-displacement relationship of the FP bearing is governed by the coefficient of sliding friction. At a given instant of time the coefficient of sliding friction depends on the instantaneous sliding velocity, instantaneous axial pressure, and the temperature of the sliding interface.

(19)

where,
is the coefficient of sliding friction at given time
is the instantaneous sliding velocity
is the instantaneous axial pressure, and
is the instantaneous temperature of the sliding surface.

Kumar et al. (2015b) developed a framework to account for the variation in due to the combined effect of velocity, pressure, and temperature and it is adopted here. The sliding velocity and the axial pressure on the bearing are governed by the response of the superstructure. The temperature at the sliding surface is a function of the coefficient of friction, sliding velocity, instantaneous axial pressure and, thermal and heat flux parameters of the material (typically stainless steel):

(20)

where,
k is the thermal conductivity of the sliding interface material, and
D is the thermal diffusivity of the sliding interface material.

It is evident from Eq. (19) and Eq. (20), that and are interdependent. This complex interdependence of the coefficient of sliding friction, sliding velocity, instantaneous axial pressure, and the temperature at the sliding surface is described using the flowchart of Figure 9.

Figure 9: Interdependence of parameters governing the lateral force-displacement relationship of a FP bearing (Kumar et al., 2015b).

Dependence of the coefficient of friction on the sliding velocity

Mokha et al. (1988) proposed an expression to account for the dependence of the coefficient of friction on the instantaneous sliding velocity as

(21)

where,
is the friction coefficient at a very small sliding velocity
is the friction coefficient at a very high sliding velocity, and
is a rate parameter.

Based on the experimental data, the ratio is assumed to be a constant value of 0.5 (Kumar et al., 2015b). The modified expression for coefficient of friction as a function of instantaneous sliding velocity, at constant axial load and temperature, can be written as:

(22)

Dependence of the coefficient of friction on the axial pressure

Kumar et al. (2015b) proposed an expression for the coefficient of sliding friction as a function of axial pressure at very high sliding velocity, and is given by

(23)

where,
is the reference axial pressure
is the co-efficent of sliding friction when the axial pressure on the bearing () is equal to the reference pressure (), and
and are constants.

The assumption of for all values of axial pressure allows the coefficient of sliding friction to be expressed as the product of the coefficient of friction measured at a reference pressure and a term that depends on the instantaneous axial pressure. Based on the experiments, and can be set equal to 0.7 and 0.02 respectively (Kumar et al., 2015b). Rewriting the expression for as a function of :

(24)

Dependence of the coefficient of friction on the temperature at sliding interface

The coefficient of sliding friction decreases as the number of cycles of loading increase. Kumar et al. (2015b) proposed an expression to define the coefficient of friction as a function of the instantaneous temperature, at a reference axial pressure and very high sliding velocity as

(25)

where,
is the coefficient of friction at a temperature
is the coefficient of friction measured at the reference temperature , and
is the instantaneous temperature.

During the analysis, it is difficult to monitor the temperature over the entire sliding surface. To address the heating, Kumar et al. (2015b) proposed an alternate approach namely, to monitor the temperature at the center of sliding surface only. When the slider is above the center of sliding surface, the temperature increases, and drops as the slider moves away. This increase in temperature at the monitoring location is governed by the heat flux parameters and the path of the slider, and is given by Eq. (26)

(26)

where,
is the increase in temperature
is the depth measured from the sliding surface,
is the thermal diffusivity of the sliding interface material,
is the thermal conductivity of the sliding interface material, and
is the heat flux.

At a given instant of time, the instantaneous heat flux, at the monitoring location can be expressed as

(27)

Dependence of the coefficient of friction on the combined effects of temperature, axial pressure and sliding velocity

During an earthquake, the sliding velocity, axial pressure and the temperature at the sliding surface of a FP bearing changes continuously. The combined effect of all the three has to be considered to determine the coefficient of sliding friction. From Eq. (22), Eq. (24),and Eq. (25), the factors accounting for friction dependence on the velocity, pressure, and the temperature () can be written as

(28)

(29)

(30)

The coefficient of sliding friction as a function of the instantaneous sliding velocity, axial pressure on the bearing and the temperature at the center of the sliding surface is updated as:

(31)

where,
is the factor accounting for friction dependence on the sliding velocity
is the factor accounting for friction dependence on the axial pressure on the bearing
is the factor accounting for friction dependence on the temperature at center of the sliding surface, and
is the coefficient of sliding friction at refernce pressure , measured at high velocity of sliding with the temperature at the center of the sliding surface being .

Nonlinear Fluid Viscous Damper

ComputeFVDamperElasticity material is used to simulate the hysteretic response of a nonlinear FVD element in MASTODON. Nonlinear fluid viscous dampers (FVDs) are seismic protective devices used to mitigate the effects of intense earthquake shaking on engineered structures such as nuclear power plants. Fluid viscous dampers operate on the principle of fluid flow through the orifice creating a differential pressure across the piston head and develops an internal resisting force in the damper. Reinhorn et al. (1995) proposed a simplified expression for the force in the damper as a function of the fractional power of relative velocity as shown in Eq. (32)

(32)

where, is the damping coefficient, is the velocity exponent, and is the instantaneous relative velocity in the damper.

The response of the FVD is characterized by two parameters, , and . For = 1, Eq. (32) models a linear damper and the typical range for is between 0.3 and 1.0 for seismic applications. Fluid viscous dampers are typically installed in the structure with an in-line brace, which provides stiffness to the damper assembly. Reinhorn et al. (1995) observed that such stiffness had to be addressed for response calculations and therefore modeled explicitly. The brace can be modeled as a spring in series with the nonlinear dashpot, and the combined behavior of spring and dashpot is best represented using the Maxwell formulation, as shown in Figure 10.

Figure 10: Mathematical model for a nonlinear FVD element (Akcelyan et al., 2018).

The Maxwell model of Figure 10 also includes the flexibility of the gussets, brackets and clevises. The overall stiffness of the damper assembly is calculated as

(33)

In the Maxwell model, the nonlinear dashpot and the spring are in series and therefore the force in the dashpot and spring are equal. This force is given by,

(34)

and the expressions for the total displacement and the total velocity of the damper are

(35)

(36)

where, and are the relative displacement and the relative velocity in the spring, and and are the relative displacement and the relative velocity in the dashpot.

From Eq. (34), the rate of change of force can be written as

(37)

Using Eq. (35) & Eq. (36), and rewriting the expression for the rate of change of force

(38)

Eq. (38) is a first-order differential equation in the form of an initial value problem (IVP)

(39)

Dormand and Prince (1980) proposed an iterative integration scheme known as Dormand Prince (DP54), to solve a generalized initial value problem of the form shown in Eq. (39) above. Iterative algorithms such as DP54 are computationally more efficient than traditional integration schemes for numerical solution of initial value problems. Akcelyan et al. (2018) proposed a framework to obtain the numerical solution of nonlinear FVD based on DP54 algorithm and that has been adopted here. For detailed information on the numerical implementation, specific to the FVD element, refer to Akcelyan et al. (2018).