Towards a Unified Single Analysis Framework Embedded with Multiple Spatial and Time Discretized Methods for Linear Structural Dynamics

David Tae and Kumar K.Tamma

Department of Mechanical Engineering,University of Minnesota,Twin Cities,Minneapolis,55455,USA

ABSTRACT We propose a novel computational framework that is capable of employing different time integration algorithms and different space discretized methods such as the Finite Element Method,particle methods,and other spatial methods on a single body sub-divided into multiple subdomains.This is in conjunction with implementing the well known Generalized Single Step Single Solve(GS4)family of algorithms which encompass the entire scope of Linear Multistep algorithms that have been developed over the past 50 years or so and are second order accurate into the Differential Algebraic Equation framework.In the current state of technology,the coupling of altogether different time integration algorithms has been limited to the same family of algorithms such as the Newmark methods and the coupling of different algorithms usually has resulted in reduced accuracy in one or more variables including the Lagrange multiplier.However,the robustness and versatility of the GS4 with its ability to accurately account for the numerical shifts in various time schemes it encompasses,overcomes such barriers and allows a wide variety of arbitrary implicit-implicit,implicit-explicit,and explicit-explicit pairing of the various time schemes while maintaining the second order accuracy in time for not only all primary variables such as displacement,velocity and acceleration but also the Lagrange multipliers used for coupling the subdomains.By selecting an appropriate spatial method and time scheme on the area with localized phenomena contrary to utilizing a single process on the entire body,the proposed work has the potential to better capture the physics of a given simulation.The method is validated by solving 2D problems for the linear second order systems with various combination of spatial methods and time schemes with great flexibility.The accuracy and efficacy of the present work have not yet been seen in the current field,and it has shown significant promise in its capabilities and effectiveness for general linear dynamics through numerical examples.

KEYWORDS Time integration; structural dynamics; multiple scale and multiple methods; ordinary differential equations;differential algebraic equations

1 Introduction

As numerical analysis is utilized in the various fields and applications,the desire for analyzing complex problems more efficiently has grown.The Finite Element Method(FEM)is one of the most widely used numerical analysis methods.While it features simplicity and versatility,it is susceptible to numerical errors such as shear locking and complexity in adopting it for large deformation problems and the like.

Concurrently,the strong form based particle methods has been gaining interest in the field of numerical simulations as the particle methods have the capability and flexibility for handling large distortion,crack propagation and free surface detecting,etc.Some of the notable strong form particle methods,Smoothed Particle Hydrodynamics(SPH)[1–4]and Moving Particle Semi-implicit Method(MPS) [5,6],have been used to solve partial differential equations.Later,Silling et al.proposed a nonlocal extension of classical continuum mechanics called peridynamics(PD)[7,8],which represents spatial derivatives of field quantities in terms of internal interactions between points/particles via integrals involving differences.The peridynamics model shows a comprehensive ability to solve linear and nonlinear second order system problems and simulate structural dynamics and molecular dynamics as well as multi-scale problems[9].More recently,a mixed form representation of the strong form particle method,where multiple primary variables are solved at once,is presented in[10].Other related works appear in references[11–14].

As problems become more complex,controlling the spatial discretization locally is a desirable feature.While an area where fast dynamics are presented is discretized with high spatial resolutions,a coarser mesh may be sufficient in the majority of the structure.In addition,the selection of spatial methods in different parts of a body can also be an effective strategy such as by utilizing the FEM for the majority of a body,and applying particle based methods only on an area where localized phenomena such as crack or fracture may occur.

One of the ways to locally control spatial discretization is dividing a body into regions of subdomains in which different methods can be applied in each subdomain.The challenge in using subdomains is handling the interface condition between the subdomains,and over the years,numerous coupling methods have been developed.Interface coupling is largely divided into two categories:overlapping and non-overlapping domain coupling.In the case of the overlapping domain decomposition,a portion of the subdomains are overlapped and the shape function similar to the framework of the FEM is used to constrain the common Degrees of Freedom (DOFs).However,the large difference in each of the subdomain’s stiffness may lead to the global stiffness matrix to be ill-conditioned,as noted by Orsini et al.in[15].Moreover,Metsis et al.[16]implied that the weak form based domain decomposition methods may be more prone to this issue due to the implementation of the Gaussian integral.

On the other hand,the non-overlapping technique is more straight forward where the subdomains are only connected by an edge or a surface in 2D or 3D,respectably,and are constrained by the Lagrange multiplier at the interface.At the interface,usually the nodes/particles are overlapped;however,Park et al.[17]and Brezzi et al.[18]introduced an approach where the subdomains are not directly interacting but through an intermediate interface.This approach allows more flexibility in meshing as one-to-one matching of the nodes is not required.

Implementing the different spatial discretization methods for localized phenomena has been a study which has been explored for many years.The coupling of the FEM with molecular dynamics via a bridging domain which enforces compatibility by utilizing the Lagrange multiplier is proposed by Xiao et al.[19].Alternately,the Discrete Element Method (DEM) [20] is a widely used meshfree method in the rock engineering simulation field and various FEM-DEM coupling methods have been developed [21,22].Utilizing the peridynamics method for localized phenomena such as crack propagation is also a topic with extensive studies including coupling the FEM and the PD[23–25].In this work,one to one node/particle matching is used to simply demonstrate the concept of coupling different spatial methods and different time integration schemes unlike state of art which is limited to just the spatial methods.We couple the FEM and particle methods with different time schemes for each subdomain and such features are not feasible in the state of art.

There also has been high interest in coupling subdomains for transient analysis.Gravouil et al.presented the GC method in[26],employing different time increments in subdomains using the limited Newmark family of algorithms.However,it was later discovered that artificial numerical dissipation arises at the interface resulting in the loss of energy and loss of second order accuracy in the GC method.

Similarly,Pegon extended the GC method and proposed the PM method which employs an interfield parallel solution procedure to couple different subdomains[27,28].Prakash et al.[29]developed the PH method in order to alleviate the GC method’s energy issue by solving the interface at a large time scale.The drawback of the PH method was described by Karimi et al.who stated in[30],that the PH method can only apply to two subdomains.Karimi et al.[30] developed the subdomain DAE framework for the Newmark family of algorithms.A noteworthy aspect of this work is that the constrained continuous system in time is in fact a system of differential algebraic equations rather than ordinary differential equations in time which most coupling methods have been based upon.Nakshatrala et al.later applied the DAE approach to the Finite Element Tearing and Interconneting(FETI) method based domain decomposition to solve time dependent first order systems with the fourth order Nakshatrala et al.[31].More recently,Shimada [32] integrated the Generalized Single Step Single Solve(GS4)family of algorithms with the subdomain DAE framework using the FEM.And,Maxam et al.[33]proposed an unified approach of the DAE-GS4 framework on thermoelastic problems using the FEM.There exists interest in coupling different subdomains for transient analysis.However,the current algorithm technology is very limited and does not permit features to interface a wide variety of time algorithms and spatial methods.And,this is the focus of the advances and contributions of this paper.

In this work,the subdomain DAE framework is particularly advanced to the Generalized Single Step Single Solve(GSSSS or GS4)framework and family of algorithms developed by Tamma et al.[34–37].The versatility of the GSSSS family of algorithms in the treatment of numerical dissipation,overshoot and ease of adapting to explicit methods,proves extremely useful to the subdomain framework.In addition,the coupling of the FEM based on the weak form,and the generalized particle method (GPS) based on the strong form is also presented in this work to demonstrate integrating altogether different spatial methods as well.The coupling of different spatial and time discretization methods may lead to more accurate computation as it allows the selection of the appropriate method for a specific situation.The novel GS4 framework and its family of algorithms is a general purpose structure.It encompasses the entire class of Linear Multistep (LMS) methods including those that have been developed over the past 50 years or so and in addition,also includes new and advanced optimal designs of algorithms,in 1st/2nd order systems that are second order time accurate and used for commercial software to enable practical real word applications to be conducted.In addition,under one umbrella,the frameworks inherit and have new and optimal algorithms and designs in addition to existing methods such as Crank-Nicolson,Gear’s BDF,etc.for first order systems; and new and optimal choices in second order systems as well in contrast to existing ones such as the Newmark family,HHT,WBZ and the like.A particularly unique feature is that the framework permits interfacing a wide variety of implicit-implicit,implicit-explicit,and explicit-explicit algorithm designs in a single analysis unlike the existing methods which are limited and not capable of doing so;in particular coupling of different implicit-implicit algorithms with/without controllable numerical dissipation which is not feasible in the current state-of-the-art.This is a significant breakthrough and is capitalized herein.The generality,flexibility and robustness of the present computational framework which embraces subdomain and time integration algorithms and differing spatial methods fused into one is not possible with today’s traditional methods and technology that shows limitations such as reduced order accuracy,consistency problems and the like in various field variables such as displacement,velocity,acceleration and Lagrange multiplier.The present developments provide a wide variety of choices to the analyst.In this work,we consistently preserve and maintain second order time accuracy in all variables such as displacement,velocity and acceleration including the Lagrange multipliers.In this paper,the FEM and GPS are first reviewed to explain the basics in integrating different spatial methods as well as different time integration methods in a single analysis.Then in the following section,the incorporation of the GS4 family of algorithms with the DAE framework is investigated and validated through various numerical examples.

2 Finite Element Method

A brief overview of the finite element method is reviewed in this section for better understanding of the various numerical examples and its methodology that are presented in a later section.In addition,the objective is to get clarity of the computations involved in the calculations of the FEM integrals and particle method integrals.

2.1 Finite Element Approximations

In the finite element method,a body of domainΩis divided into numerous regions called elements formed by connecting points known as nodes.The finite dimensional field variable such as the displacementuhis approximated within each element by the basis functions.The most common basis functions have the Kronecker deltaδijproperty so that the nodal values ?ucan be accurately represented.The integration of the finite dimension is performed by summing up the integrals over each element domain.

whereNiare the elemental basis functions which varies for different types of element.The basis functions for a commonly used 2D bilinear quadrilateral are given as follows:

where

Quadrilateral elements with two Gauss points in each x and y direction will be used throughout in the numerical example section.

2.2 Weak Form Formulation for FEM

The equilibrium equation expressed in terms of the displacement can be written as(shown for 1D as a simple illustration)

Next we introduce an arbitrary functionw(x) that is defined in the domainΩand multiply the equilibrium equation.

Then we integrate over the domain to obtain an integral form that is set to zero.

The stress term is then integrated by parts as

whereΓis the boundary ofΩandnxis the outward pointing normal to the boundary.By substitution of Eq.(7)into Eq.(6),we get the weak form of the equilibrium as

Introduce the constitutive equation to obtain governing equation based on displacement,and we get

It is worth noting that we removedtxby settingw=0 where the displacements are specified.

The resulting weak formulation in space is given as

and limiting the trial and test spaces to high fidelity finite element spaces,we get

Settinguh=Ni?uhiandwh=Ni?whi,for the Ritz-Galerkin approximation,we have

The above equation is valid for any ?wh.Therefore we may write

leading to

3 Generalized Particle Method

3.1 Gradient and Divergence of Classical Physical Field

A physical field can be thought of as the assignment of a physical quantity at each point in space and time.Consider a domain in Euclidean space,E.Let X ∈Ω0?E denotes the geometric position vector of a certain position in the domain and∈Ω0?E denotes another point which is close to the position X as shown in Fig.1a.Therefore,any physical quantity at each point in space and time can be defined as?(X,t):Ω0×I→E.DefineΩX?Ω0as a domain where the center point position vector X is in,then we have

Figure 1:(a)Continuous domain with position vector X and ,(b)Discretized domain with position vectors Xi and Xj

The weighted residual method has been widely used in solving different kinds of ordinary/partial differential equations,in particular,using the Galerkin method,which is the fundamental theory of the finite element method.In this work,we exploit the weighted residual method directly to the Taylor series expansion,such that a generalized approach of developing particle based methods can be explored.

In mathematics,the Taylor series expansion is a representation of a function as an infinite sum of terms that are calculated from the values of the function’s derivatives at a single point.The generalized Taylor series of a general physical quantityabout X at a fixed timetyields

where ??(X,t)is a vector when?(X,t)is a scalar value and is a second order tensor when?(X,t)is a vector value;andΛi∈R(i=1,2,3....)is the algorithmic parameters.

Omitting the higher order terms from Eq.(19),we have the first order Taylor series approximation

whereΛ1=1.Define a residual with regards to Eq.(20),R,as

According to the standard weighted residual method,we introduce a weight function,C1In this work,the weight function used in the weighted residual method is defined as C,and the weight function for describing the influence generated from neighboring particle is defined as W.,and minimize the residual.Eq.(21)is multiplied by the weight functionCand integrated over the domainΩX,which then yields

It is worth noting that the domain is the cut-off influence domain of the positionX,and the residual is a local value.The weight functionCprovides a wide scope of approximations for constructing the derivative definitions used in various particle-based methods and new approach of particle methods as well.

GradientRecall the residualRin Eq.(21)and the weighted residual,RC,with respect to a vector property?is

where C is a vector with the size of geometric dimension of a problem;for example,Cis a two by one vector in 2D problems.We postulate that Eq.(23)is the exact relationship between the point X and one of its neighboring pointlocated in the domainΩXin terms of the weighted residual,.Therefore,by integrating Eq.(23)over the domainΩX,the following formulations can be obtained as

Hence,the gradient of a vector can be obtained as

The proposed gradient operator has the ability to recover the′?′in the Moving Particle Semiimplicit (MPS) method,the deformation gradient in the state-based peridynamics as well as in the corrected Smoothed Particle Hydrodynamics(SPH)method,by selecting the appropriate correspondingC.

Remark 3.1.By selecting the arbitraryC=nij,the Gradient in terms of the original MPS method can be derived from the proposed approach as well.SubstitutingC=nijin the discretized system yields,

Discretizing into particles,and applyingnij=(Xj?Xi)/it becomes

Thus,Eq.(27)can be written as

Gradient Operator:

Introducing the weighting function Wijand supposingVi=Vj,VΩXi=VΩXjyields

It is worth noting that for a general particle within a uniformly distributed system,the weighted average matrixAjkis a diagonal matrix.Consequently,Ajk=dI,where I is identity matrix and d equals to the dimensional number.Therefore,Eq.(31)recovers the gradient calculation in the MPS method,and can be rewritten as

where n0=

DivergenceBy representing the divergence model as gradient operator,?,dot product with a physical quantity,?,the divergence model can be written as following

The discretized formulations is listed as following

LaplacianThe definition of the Laplacian can be written as

In order to obtain the Laplacian relationship between the center point X and its neighboring point,we reconstruct the gradient between two points via taking advantage of Eq.(26).One may notice that it is convenient to find the gradient from Eq.(23)by setting theas 0,which yields

Therefore,an alternative gradient of a vector is given by taking advantage of the average property,Eq.(37).

Consequently,the Laplacian can be obtained via combining Gauss’s Theorem and the gradient formulations.

Therefore,the Laplacian of?(X,t)becomes

Hence the following Laplacian formulations can be obtained as

The discretized formulations is listed as following:

3.2 Stiffness Matrix for Particle Method

We will explore how the stiffness matrix is formed in 2D using the GPS method.We start with the same governing equation and constitutive law,

whereDis elastic matrix as following

The first step to get the stiffness matrix is to describe ?uwith the gradient operator from the GPS method as described in Eq.(32).

where[A]is a two by two matrix that is computed first.The summation is applied for every particle within particlei’s influence circle.But for simplicity,we will consider there are two particles for the system,particles 1 and 2.The Eq.(47)will be written as following for particle 1,

Thus,

Next,we need to convert the 2 by 2 matrix on the left hand side of Eq.(50)to a 4 by 1 vector form.

We can also get the engineering strain as follows:

Eq.(52)can be simply written as

Stressσcan be computed by multiplying elastic matrix,[D],from Eqs.(45)and(46).

It is worth noting that[D][?][u]includes all the components in the particlei’s influence domain.

where[?1]and[?2]represent the gradient operator obtained from the first and the second particle as theith particle,respectively.This suggests that,in general,the gradient operator for every particle in the system needs to be computed before applying divergence operator.

Next,we will show how to construct the divergence operator into a matrix form similar to how the gradient operator matrix was formed as shown above.

whereB1=andB2=for simplicity.Here we further isolateσtogether.

or for the engineering stress,

Thus,we can get the stiffness matrix for the particle 1,by combining Eqs.(55)and(60),[?·][D][?][u]=F

This process is repeated for each particle in the system.After generating [Ki] matrices,they are assembled together to form the system stiffness matrix[K].

4 Equations of Motion for a System of Multiple Subdomains

Differential Algebraic Equations (DAE) are differential equations with algebraic constraints.Consider the equation of motion of a body that is subject to constraintsg(u)=0.The action integral is given as

and we can take the variation of the action integral,

whereG(u)=?g(u)/?u.For the standard DAE equations of second order transient systems,the above equation can be written as

where,Qiner,Qintand Qapplare the inertial,internal and applied forces,respectively.The Lagrange multiplier variable,λ?(t),is introduced in order to accomplish the constraint equation.

Now,consider a domain of interestΩthat is divided intondomnumber of subdomains,Ωi(fori=1,2,...,ndom),with each subdomain joined by a non-overlapping interface in between.An example with three subdomains is shown in Fig.2.In the case of linear semidiscrete second order transient systems,the equation of motion for subdomainΩiis,(fori=1,2,...,ndom),

where ui(t),(t):=dui/dt,and(t):=d2ui/dt2are the displacement,velocity,and acceleration vectors of theith-subdomain,respectively;Miis the mass matrix of theith-subdomain;and Qiint(ui)andare the internal and applied force of theith-subdomain,respectively.

Figure 2:A domain Ω divided into three subdomains with non-overlapping interface

The matrix Giis a boolean matrix with entries either?1,0,or+1,that is defined by the subdomain interface conditions.The nodes or particles on one side of the interface will have?1 while the corresponding nodes or particles on the other side of the interface will have value of +1 and all other nodes/particles will have 0.Eq.(67b)indicates that the velocity continuity along the interface is enforced.

4.1 DAE Framework with GS4 Family of Algorithms for a System with Multiple Subdomains and Multiple Time Step Sizes

LetΔtandΔtibe the system time step size and the subdomain time step size forΩi,respectively,and define the ratio asτi=Δt/Δti≥1 fori=1,2,...,ndom.The system and subdomain time step sizes are defined as

wherendenotes the system time step,andmrepresents subdomain time step.The visual representation of time subcylcing is shown in Fig.3.

Figure 3:Time step subcycling for different subdomains

Now applying the GS4 family of algorithms,the equation of motion of each subdomain can be written as following:

where

It is worth noting that the Lagrange multipliers are approximated in the system time level rather than in the subdomain time level.This factor suggests that the solution must be in the system time level even though the displacement,velocity,and acceleration are approximated in the subdomain time step.

The updates of each of the unknown variables with the GS4 framework of algorithms are

WΛandλare the GS4 algorithmic parameters and can be found in [37].The advantage of the GS4 algorithm is that it allows one to easily switch between various time integration algorithms as well as explicit and implicit families based on the GS4 algorithmic parameters.A single analysis code readily enables current and future researchers and analysts to investigate a multitude of algorithms and combinations in the academic and commercial software world;thereby making it a much desired computational simulation toolkit.

In linear dynamical systems,Eq.(67)may be written as

whereCiand Kiare the damping and stiffness matrices forΩi,respectively;and Qappli (t)is the external force vector forΩi.Since the stiffness matrix is different for each subdomain,various methods can be coupled together with each subdomain having the stiffness matrix derived from the corresponding methods.For example,subdomain 1 and 2 in Fig.2 can use the stiffness matrix derived from the FEM while subdomain 3 adopts the stiffness matrix formulated using the GPS method as outlined in the previous section.

A set of all kinematic unknowns forΩiin its subdomain time step is defined as

Eq.(71a)can then be discretized and written as follows:

where

with

and

where O represent zero-matrix,ncandnirefer to the number of constraint degrees of freedom,and total number of degrees of freedom in one domain,respectively.Now,we define a set of all kinematic unknowns,i.e.,the acceleration,velocity,and displacement vectors,forΩiover its subdomain time step and a set of all kinematic unknowns of all subdomains in the system over a system time step as

The constraint equation at the velocity level,i.e.,Eq.(71b),can be written in the form

where

and

with

Eq.(82) shows the velocity constraint.The size of Giisncby 3niτisince the system solves for the displacement,velocity and acceleration at every subcycling step simultaneously.However,Eq.(81)suggests that the constraint is enforced only at then+1 time step.

The equation of motion for the entire system for a system time step level can be represented in a simple matrix form as following:

where

in which

in which

and

in which

4.2 Computational Procedure

The computation flowchart for the GS4-II DAE framework for the subdomain method employing different spatial methods,algorithms,and time steps is shown below in Fig.4 which here which provides a wide variety of algorithms and choices as options to the analyst.

?for dom:=1 to Ndom do Calculate M,C,K using corresponding method for the domain Calculate Li,and Ri using Eqs.(74)and(75a)Calculate ?Gi using Eq.(82)for τ:=1 to n do Assemble Ai using Eq.(85)Assemble Bi using Eq.(87)end for end for Assemble A using Eq.(84)Assemble B using Eq.(86)Transpose B to obtain G for(t:=1 to Timesteps)do for dom:=1 to Ndom do Calculate Fi using Eq.(89)end for Assemble F using Eq.(88)Assemble ∧Rglob and ?Jglob from Eq.(83)Solve(?JglobXn+1=∧Rglob)end for

Figure 4:Flowchart for constructing DAE second order system

4.3 Subdomain Coupling Condition and Time Level Consistency for Second Order Systems

The second order time accuracy of all the primary variables and the Lagrange multipliers can be achieved for interfacing different subdomains(this is in general not plausible with traditional methods)ifW1,a parameter in GS4,is the same for all subdomains.In order to guarantee the second order time accuracy in the acceleration,the algorithmic accelerationmust be an approximation at time leveltn:

whereφi:=W1Λ6?W1in whichW1andW1Λ6are the values used for subdomainΩi.

The GS4-II algorithms are written in the form of U0(ρmin,ρmax,ρs) for the U0 family and V0(ρmin,ρmax,ρs)for the V0 family.The relation of the variables follows 0≤ρs≤ρmin≤ρmax≤1.And depending on the value given to each variable,it recovers various existing time schemes including a wide variety of new and optimal methods and also provides new avenues to integrate multiple methods in a single analysis;for example,U0(1,1,0)recovers Newmark time scheme,and U0(0,0,0)is equal to the WBZ scheme.The V0 family represents an altogether new family of algorithms and is not well known to others and is relatively new and not well understood by the traditional research community at large;mathematically,the V0 family of algorithms are the mirror counterparts of the U0 family but with different overshoot behavior unlike several of the traditional methods.More recent efforts show the advances in the GS4 V0?family of algorithms that enable consistent second order time accurate features with and without physical damping[38].Certain algorithms within the U0 family are also new with certain exceptions which are existing methods.The GS4 encompasses the entire family of LMS methods that are second order time accurate.Only three parameters are needed to identify a particular algorithm;and no two algorithms can have the same DNA(Discrete Numerically Assigned markers).

With theW1parameter condition,the GS4-II family of time integration scheme combinations guarantee second order time accuracy of all variables and the Lagrange multipliers; this is unlike existing technology which does not permit such implicit interfacing of arbitrary methods whilst consistently preserving second order accuracy.Some of the possible GS4-II time scheme combinations are listed in Table 1.For example,one cannot simply or haphazardly couple the Midpoint rule with the Newmark method; or couple the Newmark,HHT,WBZ methods,etc.,arbitrarily in different subdomains in a single analysis without loss and reduced accuracy in the convergence rate of one or more variables.

Table 1: Coupling different time scheme algorithms

For the U0 family of GS4-II algorithms,theW1value is determined byρsparameter; whereas,for the V0 family,ρminandρmaxparameters control theW1value as shown below.It is also possible to couple U0 and V0 family if theW1condition is satisfied.For example,U0(1,1,0)scheme and V0(1/3,1/3,1/3)scheme both produceW1value of 1;thus,they can be coupled together.These criteria must be followed to preserve second-order time accuracy in all the variables and the Lagrange multipliers.

5 Numerical Examples

5.1 2D Wave Propagation Problem

Fig.5 shows the 2D geometry of a wave propagation problem.It is a rectangular bar of 10 m by 1 m with 1 m thickness,which is fixed on the left side wall as boundary condition.The semidiscrete equation of motion for the system is given as

Figure 5:Diagram of structure for wave propagation problem

The bar is compressed from the right side with p,the distributed load,of 10,000 N/m.This material has 72 GPa for Young’s modulus,0.3 for Poisson’s ratio,and 2712kg/m3for density.The bar is divided into two subdomains,and each subdomain is implemented with 17 nodes/particles in length,and 5 nodes/particles in height making it total of 170 nodes/particles with 5 nodes/particles overlap on the interface.Figs.6 and 7 show the FEM mesh and particle distribution for the two subdomains.The systemΔtfor time integration is set to be 0.8μs.The time history of the x-direction displacement,velocity and acceleration of a node/particle at the interface on the subdomain 1(SD1)side(point A in Fig.5)and the free end on the subdomain 2(SD2)side(point B in Fig.5)are evaluated.We compare the result in four sets:spatial method variance where we combine the different spatial methods;explicit and implicit variance where we pair explcit and implicit version of the same time scheme;time shceme variance where the different combination of the GS4-2 time schemes are tested;and time subcycling variance.

Fig.8 compares the three different cases of coupling different methods:FEM+FEM,GPS+GPS,and FEM+GPS; and shows the time history of the displacement,velocity,and acceleration at A in subdomain 1 and at point B in subdomain 2.It is clear that coupling the same methods or coupling different methods,be it weak form or strong form methods together,produce a comparable result with the result from commercial software such as ABAQUS.

Figure 6:FEM+FEM and FEM+GPS interface

Figure 7:FEM+GPS two domain mesh of the structure

Fig.9 shows the result for coupling explicit version and implicit version of the schemes.The evaluated coupling cases are explicit and explicit case,explicit and implicit case,and implicit and explicit case where the first condition refers to the time scheme algorithm for subdomain 1 and the latter for subdomain 2.This figure not only suggests that explicit and implicit time algorithms can be coupled together,it also shows that coupling is not limited to just Newmark family of algorithms as this example couples U0(1,1,0)(equivalent to Newmark,central difference for explicit,and average acceleration for implicit)and U0(0,0,0)(equivalent to WBZ).

To further present evidence of coupling various time scheme algorithms,five different cases listed in Table 1,are evaluated and illustrated in Fig.10.The time scheme algorithms are selected such that a combination has identicalW1to ensure the second order accuracy in time.It is evident that all the results align closely.

Figure 8: Comparison of transient response of displacement,velocity and acceleration (a) at point A on SD1 and (b) at point B on SD2,between three cases of different pairing of methods (Case 1:FEM+FEM,Case 2: FEM+GPS,Case 3: GPS+GPS) while using explicit version (with η3=0) of GS4-II U0(1,1,0)algorithm for subomain 1 and explicit version(with η3=0)of U0(0,0,0)algorithm for subdomain 2 with constant Δt in each case

Figure 9:Comparison of transient response of displacement,velocity and acceleration(a)at point A on SD1 and(b)at point B between different pairing of explicit(with η3=0)and implicit version(with η3=1)of GS4-II(Case 1:explicit+explicit,Case 2:explicit+implicit,Case 3:implicit+explicit)while using FEM with GS4-II U0(1,1,0)algorithm for subomain 1 and GPS with U0(0,0,0)algorithm for subdomain 2 at constant Δt in each case

Figure 10: Comparison of transient response of displacement,velocity and acceleration (a) at point A on SD1 and (b) at point B between different pairing of explicit version (with η3=0) of GS4-II algorithms listed in Table 1,while using FEM for subdomain 1 and GPS for subdomain 2 at constant Δt in each case

Lastly,the subcycling cases are also evaluated in Fig.11.A case where subdomain 2 has double the time steps than subdomain 1(time increment interval in subdomain 2 is half of that of subdomain 1)is compared with the no subcycling case and the results are very close.

Figure 11:Comparison of transient response of displacement,velocity and acceleration at both point A and point B between simulation with system Δt for both subdomains and with variable Δt=1/2 system Δt for subdomain 2 while using FEM with explicit version (with η3=0) of GS4-II U0(1,1,0) algorithm for subdomain 1 and GPS with explicit version (with η3=0) U0(0,0,0) algorithm for subdomain 2

Fig.12 demonstrates second order time accuracy in displacement,velocity,acceleration and the constraint variableλ?for different coupling conditions shown through Figs.8–10.

Figure 12: (Continued)

Figure 12: (Continued)

Figure 12:Convergence order for the solution with respect to the system time step for domain 1 and 2 in different cases:(a)FEM with explicit version of GS4-II U0(1,1,0)and FEM with U0(0,0,0),(b)GPS with explicit version of GS4-II U0(1,1,0)and GPS with U0(0,0,0),(c)FEM with explicit version of GS4-II U0(1,1,0)and GPS with U0(0,0,0)algorithm,(d)FEM with explicit version of GS4-II U0(1,1,0)and GPS with implicit version U0(0,0,0),(e)FEM with implicit version of GS4-II U0(1,1,0)and GPS with explicit version U0(0,0,0),(f)FEM with explicit version of GS4-II U0(1,1,0)and FEM with implicit version U0(0,1,0),(g)FEM with explicit version of GS4-II U0(0.5,1,0.5)and GPS with U0(0.5,0.5,0.5),(h)FEM with explicit version of GS4-II V0(1,1,0)and GPS with V0(1,1,1),(i)FEM with explicit version of GS4-II V0(0.7,0.7,0.7)and GPS with V0(0.7,0.7,0.21),(j)variable Δt=1/2 system Δt for subdomain 2 while using FEM with explicit version(with η3=0)of GS4-II U0(1,1,0)algorithm for subdomain 1 and GPS with explicit version(with η3=0)U0(0,0,0)algorithm for subdomain 2

Figure 13:Diagram of structure for bending problem

5.2 2D Vibration Problem

A beam with the same geometry as the one in the wave propagation problem is used to illustrate the application of the multi-domain method in a plane stress dynamic problem.A cantilever beam of 10 m by 1 m with 1 m thickness is fixed on the left side and is experiencing downward traction,p,of 100,000N/mas shown in Fig.13.The density,ρ,of the beam is set to be 271.2kg/m3,Young’s modulus,E,to be 72 GPa,and Poisson’s ratio to be 0.3 for the properties of the material.The beam is divided into two subdomains,and each subdomain is implemented with 17 nodes/particles in length,and 5 nodes/particles in height making it total of 170 nodes/particles with 5 nodes/particles overlap on interface.The semidiscrete equation of motion for the system is given as

The systemΔtfor time integration is set to be 0.8μs.The time history of y-direction displacement,velocity and acceleration of a node at the interface(point A)and the free end(point B)are evaluated.

Similar to the previous problem,a different combination of methods are compared in Fig.14.It is apparent that the different methods can be coupled regardless of whether one is weak form or strong form.

Figure 14: (Continued)

Figure 14: Comparison of transient response of displacement,velocity and acceleration (a) at point A on SD1 and (b) at point B between different pairing of methods (Case 1: FEM+FEM,Case 2:FEM+GPS,Case 3: GPS+GPS) while using implicit version (with η3=1) of GS4-II U0(1,1,0)algorithm for subdomain 1 and implicit version of U0(0,0,0)algorithm for subdomain 2 with constant Δt in each case

Fig.15 shows results for coupling the explicit version and implicit version of the time scheme algorithms.The evaluated coupling cases this time are the implicit and implicit case,explicit and implicit case,and implicit and explicit case.The response results for the different cases align closely.The same five cases of time scheme algorithm combinations used in the previous example are evaluated and illustrated in Fig.16.And the results for all cases match.The subcycling case with one domain having half theΔtas the other,is shown in Fig.17.It is apparent that subcycling still produces an accurate result.Lastly,all of the different cases and conditions are second order accurate in time as evident from Fig.18,illustrating slope of two in the log-log error convergence plot.

Figure 15:Comparison of transient response of displacement,velocity and acceleration(a)at point A on SD1 and(b)at point B between different pairing of explicit(with η3=0)and implicit version(with η3=1)of GS4-II(Case 1:implicit+implicit,Case 2:explicit+implicit,Case 3:implicit+implicit)while using FEM with GS4-II U0(1,1,0)algorithm for subdomain 1 and GPS with U0(0,0,0)algorithm for subdomain 2 at constant Δt in each case

Figure 16: Comparison of transient response of displacement,velocity and acceleration (a) at point A on SD1 and (b) at point B between different pairing of implicit version (with η3=1) of GS4-II algorithms listed in Table 1 while using FEM for subdomain 1 and GPS for subdomain 2 at constant Δt in each case

Figure 17:Comparison of transient response of displacement,velocity and acceleration at both point A and point B,between simulation with system Δt for both subdomains and with variable Δt=1/2 system Δt for subdomain 2 while using FEM with implicit version (with η3=1) of GS4-II U0(1,1,0)algorithm for subdomain 1 and GPS with implicit version(with η3=1)U0(0,0,0)algorithm for subdomain 2

Figure 18: (Continued)

Figure 18: (Continued)

Figure 18:Convergence order for the solution with respect to the system time step for domain 1 and 2 in different cases:(a)FEM with implicit version of GS4-II U0(1,1,0)and FEM with U0(0,0,0),(b)GPS with implicit version of GS4-II U0(1,1,0)and GPS with U0(0,0,0),(c)FEM with implicit version of GS4-II U0(1,1,0)and GPS with U0(0,0,0),(d)FEM with explicit version of GS4-II U0(1,1,0)and GPS with implicit version U0(0,0,0),(e)FEM with implicit version of GS4-II U0(1,1,0)and GPS with explicit version U0(0,0,0),(f) FEM with implicit version of GS4-II U0(1,1,0) and FEM with U0(0,1,0),(g)FEM with implicit version of GS4-II U0(0.5,1,0.5)and GPS with U0(0.5,0.5,0.5) algorithm for subdomain 2,(h) FEM with implicit version of GS4-II V0(1,1,0) and GPS with V0(1,1,1),(i)FEM with implicit version of GS4-II V0(0.7,0.7,0.7)and GPS with V0(0.7,0.7,0.21),(j)variable Δt2=1/2 system Δt for subdomain 2 while using FEM with implicit version(with η3=1)of GS4-II U0(1,1,0)algorithm for subdomain 1 and GPS with implicit U0(0,0,0)algorithm for subdomain 2

5.3 2D Vibration Problem-Three Subdomains

Next,a body divided into three domains is examined next for the extension of the subdomain coupling concept.The same vibration problem above is evaluated again but this time,the beam is divided into three subdomains whereL2=L/10.The FEM is used on subdomains 1 and 3 with 22 nodes in length,and 7 nodes in height.GPS is implemented on subdomain 2,with 7 particles in length and 7 particles in height.The time scheme algorithm used in subdomain 1 and 3 is implicit U0(1,1,0);and in subdomain 2,U0(0,0,0)is used.The time history of y-direction displacement,velocity and acceleration of a node at the free end is evaluated.

Figs.20 and 21 clearly illustrate that the result at points A and B (shown in Fig.19a) from the three subdomain case matches closely with the result from the two subdomain case while achieving second order time convergence rate.The high frequency in the acceleration plot for the three domain free end node is caused by the use of non-dissipative U0(1,1,0)scheme.

Figure 19: (a) Structure divided into three domains (b) An example of mesh of the structure: FEM GPS FEM

Figure 20:Comparison of transient response of displacement,velocity and acceleration at point A and point B between three subdomain configuration and two subdomain configuration while using FEM with implicit version(with η3=1)of GS4-II U0(1,1,0)algorithm for subdomain 1 and subdomain 3 and GPS with implicit version of U0(0,0,0) algorithm for subdomain 2 with constant Δt in all subdomains

Figure 21:Convergence order for the solution with respect to the system time step for domain 1 and 2 and 3 in three domain configuration case:(a)subdomain 1 using FEM with implicit version of GS4-II U0(1,1,0)algorithm,(b)subdomain 2 using GPS with U0(0,0,0)algorithm,(c)subdomain 3 using FEM with GS4-II U0(1,1,0)algorithm

5.4 2D Notched Specimen Problem

Figs.22 and 23 show the geometry of a 2D plane stress problem and an example of the mesh,respectively.It consists of three subdomains where identical subdomain 1 and 3 are connected by a thin subdomain 2.The subdomain 1 is 3 m by 2 m rectangle which are discretized into 176 FEM nodes.On the other hand,subdomain 2 is 1 m by 2 m rectangle bar,with notches on top and bottom,discretized into 99 particles.And subdomain 3 is identical to the subdomain 1 geometrically,but discretized with 176 particles.The whole geometry has thickness of 1 m.The semidiscrete equation of motion for the system is given as

Figure 22:Diagram of structure divided into three subdomain.The red circles represent the locations where u,˙u,and ü are tracked in each subdomain

Figure 23:Mesh of all three subdomains

The bar is fixed at the left end and pulled down from the right end with a distributed load of 210 N/m.This material has 70 GPa for Young’s modulus,0.3 for Poisson’s ratio,and 1,000kg/m3for density.The simulation ran for 0.2 s with the systemΔtfor time integration set to 0.4 ms.To demonstrate the coupling of different spatial methods and different time schemes,four cases are studied as shown in Table 2.

Table 2: Coupling different time scheme algorithms

The time history of the y-direction displacement,velocity and acceleration of nodes/particles in each subdomain(points A,B,and C in Fig.22)are evaluated.

Figs.24–26 show the transient response of the displacement,velocity and acceleration for points A,B and C for all cases listed in the Table 2 compared with the reference results.All responses from the different spatial methods and different time schemes combinations match well.

Figure 24: Comparison of transient response of displacement,velocity and acceleration at point A between different cases with varying method and time integration schemes for each domain as listed in Table 2

Figure 25: Comparison of transient response of displacement,velocity and acceleration at point B among different cases with varying method and time integration schemes for each domain as listed in Table 2

Figure 26: Comparison of transient response of displacement,velocity and acceleration at point C among different cases with varying method and time integration schemes for each domain as listed in Table 2

6 Conclusions

In this paper,we have proposed a novel implementation of the Generalized Single Step Single Solve time integration framework for ODE’s into the DAE framework.This proposed method,unlike traditional approaches in current technology with LMS methods,allows the coupling of different numerical analysis methods such as FEM and particle methods to be used on a single body while ensuring the accuracy of the result with the use of a wide variety of time integration algorithms within the GS4 framework.For example,one cannot simply or haphazardly couple the Midpoint rule with the Newmark method;or couple a multitude of schemes such as the Newmark,Midpoint rule,HHT,WBZ methods,etc.,arbitrarily in different subdomains in a single analysis without loss and reduced accuracy in the convergence rate of one or more variables.However the GS4 framework circumvents various deficiencies,preserves accuracy,and permits a much broader design space for coupling algorithms in a single analysis and provides much desired robust features for applications to large scale industrial real world problems as well.It has the potential to increase the accuracy of the physics by selecting an appropriate method for only the area with a localized phenomena rather than utilizing only a single method that may not be suitable for certain applications on the whole body.The method is verified by applying it to 2D simple bar and beam problems.The results from various combination of methods and time scheme algorithms match closely with the reference result.With the basis of the proposed computational methodology established,the extension of the method for the future studies includes computational efficiency studies,extension to nonlinearity,and implementation of reduced order modeling.This work provides generality and versatility of the computational framework incorporating a wide variety of subdomain based spatial and time integration algorithms in a single analysis with great accuracy.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.