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

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 mesh-free 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 inter-field 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 displacement uh is approximated within each element by the basis functions. The most common basis functions have the Kronecker delta δij property so that the nodal values u^ can be accurately represented. The integration of the finite dimension is performed by summing up the integrals over each element domain.

∫Ωg(u)dΩ=∑j=1Nh∫Ωjg(uh)dΩj=∑j=1Nelem∫Ωjg(Niu^i)dΩj(1)

where Ni are 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:

u=Niu^i=[N1(ξ,η)N2(ξ,η)N3(ξ,η)N4(ξ,η)][u1u2u3u4]T(2)

where

N1(ξ,η)=14(1−ξ)(1−η)N2(ξ,η)=14(1+ξ)(1−η)N3(ξ,η)=14(1+ξ)(1+η)N4(ξ,η)=14(1−ξ)(1+η)(3)

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)

∂σ∂x+bx=ρ∂2u∂t2(4)

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

g(w,u,σx)=w(x)(ρ∂2u∂t2−∂σ∂x−bx)=0(5)

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

G(w,u,σx)=∫Ωw(x)(ρ∂2u∂t2−∂σ∂x−bx)dx=0(6)

The stress term is then integrated by parts as

−∫Ωw(x)∂σx∂xdx=∫Ω∂w∂xσxdx−w(x)nxσx|Γ(7)

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

G(w,u,σx)=∫Ωw(x)(ρ∂2u∂t2−bx)dx+∫Ω∂w∂xσxdx−wtx|Γu−wt¯x|Γt=0(8)

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

G(w,u)=∫Ωw(x)(ρ∂2u∂t2−bx)dx+∫Ω∂w∂xE∂u∂xdx−wt¯x|Γt=0(9)

It is worth noting that we removed tx by setting w = 0 where the displacements are specified.

The resulting weak formulation in space is given as

ρ∫Ωw∂2u∂t2dΩ+∫Ω∇wE∇udΩ=∫ΩwbdΩ+∫∂Ω∖Γuwt¯dΓ(10)

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

ρ∫Ωwh∂2uh∂t2dΩ+∫Ω∇whE∇uhdΩ=∫ΩwhbdΩ+∫∂Ω∖Γuwht¯dΓ(11)

Setting uh=Niu^hi and wh=Niw^hi, for the Ritz-Galerkin approximation, we have

ρw^hiT∫ΩNTNu^¨hidΩ+w^hiT∫Ω∇NTE∇Nu^hidΩ=w^hiT∫ΩNTbdΩ+w^hiT∫ΓNTt¯dΓ(12)

The above equation is valid for any w^h. Therefore we may write

ρ∫ΩNTNu^¨hidΩ+E∫Ω∇NT∇Nu^hidΩ=∫ΩNTbdΩ+∫ΓNTt¯dΓ(13)

leading to

Mu^¨hi+Ku^hi=F^(14)

M=ρ∫ΩNTNdΩ(15)

K=E∫Ω∇NT∇NdΩ(16)

F^=∫ΩNTbdΩ+∫ΓNTt¯dΓ(17)

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 X_∈Ω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⊂Ω0 as a domain where the center point position vector X is in, then we have

φ(X,t)=limVΩX→0∫VΩXφ(X_,t)dVVΩX(18)

Figure 1: (a) Continuous domain with position vector X and X, (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 quantity φt(X_) about X at a fixed time t yields

φ(X_,t)=φ(X,t)+Λ1∇φ(X,t)(X_−X)+Λ2H(X,t)(X_−X)(X_−X)+…+O(‖X_−X‖n)(19)

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

φ(X_,t)−φ(X,t)≅Λ1∇φ(X,t)(X_−X)(20)

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

R=φ(X_,t)−φ(X,t)−∇φ(X,t)(X_−X)(21)

According to the standard weighted residual method, we introduce a weight function, C1, and minimize the residual. Eq. (21) is multiplied by the weight function C and integrated over the domain ΩX, which then yields

∫ΩXRCdV=0(22)

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

Gradient Recall the residual R in Eq. (21) and the weighted residual, RC, with respect to a vector property φ is

(RC)XX_=(φ(X_,t)−φ(X,t))⊗C−∇φ(X,t)XX_(X_−X)⊗C(23)

where C is a vector with the size of geometric dimension of a problem; for example, C is 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 point X_ located in the domain ΩX in terms of the weighted residual, (RC)XX_. Therefore, by integrating Eq. (23) over the domain ΩX, the following formulations can be obtained as

∫VΩXRCdV=∫VΩX[(φ(X_,t)−φ(X,t))⊗C]−∇φ(X,t)XX_[(X_−X)⊗C]dV=0(24)

By treating (φ(X_,t)−φ(X,t))⊗C and ∇φ(X,t)(X_−X)⊗C in the integral presented above as field properties in the domain ΩX, the following formulations can be obtained as

1VΩX∫VΩX[(φ(X_,t)−φ(X,t))⊗C]dV=1VΩX∫VΩX∇φ(X,t)XX_[(X_−X)⊗C]dV(25)

Hence, the gradient of a vector can be obtained as

∇φ(X,t)=limVΩX→0(∫VΩX(φ(X_,t)−φ(X,t))⊗CdV)limVΩX→0(∫VΩX(X_−X)⊗CdV)−1(26)

The proposed gradient operator has the ability to recover the ′∇′ in the Moving Particle Semi-implicit (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 corresponding C.

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

∇φi=[∑j∈ΩXi(Xj−Xi)⊗nijVj]−1[∑j∈ΩXi(φj−φi)nijVj](27)

Let A=(X_−X)⊗nij and it can be replaced by its normalized weighted average value within an infinitesimal domain via exploiting Eq. (18),

AXX_avg=1VΩXlimVΩ→0∫ΩX[(X_−X)⊗C]dV(28)

Discretizing into particles, and applying nij=(Xj−Xi)/||Xj−Xi||, it becomes

Aij=1VΩXi∑j∈ΩXi(nij⊗nij)Vj(||Xj−Xi||)=∑j∈ΩXi(nij⊗nij)WijVjVΩXi∑j∈ΩXiWij⏟NormalizedWeightedAverageA(||Xj−Xi||)(29)

Thus, Eq. (27) can be written as

∇φij=[∑k∈ΩXj(njk⊗njk)WjkVkVΩXj∑k∈ΩXjWjk]−1(φj−φi)nij||Xj−Xi||(30)

Gradient Operator:

Introducing the weighting function Wij and supposing Vi=Vj,VΩXi=VΩXj yields

∇φi=1VΩXi∑j∈ΩXiWij∑j∈ΩXi[∑k∈ΩXj(njk⊗njk)WjkVkVΩXj∑k∈ΩXjWjk]⏟Ajk−1(φj−φi)nij||Xj−Xi||WijVj=1∑j∈ΩXiWij∑j∈ΩXi[∑k∈ΩXj(njk⊗njk)Wjk∑k∈ΩXjWjk]−1(φj−φi)nij||Xj−Xi||Wij(31)

It is worth noting that for a general particle within a uniformly distributed system, the weighted average matrix Ajk is 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

∇φi=dn0∑j∈ΩXi(φj−φi)||Xj−Xi||2(Xj−Xi)Wij(32)

where n0=∑jWij.

Divergence By representing the divergence model as gradient operator, ∇, dot product with a physical quantity, φ, the divergence model can be written as following

∇⋅φ(X,t)=∇Tφ=limVΩX→0(∫VΩX(φ(X_,t)−φ(X,t))⋅CdV)limVΩX→0(∫VΩX(X_−X)⊗CdV)−1(33)

The discretized formulations is listed as following

∇⋅φi=[∑j∈ΩXi(φj−φi)⋅CjWijVj][∑j∈ΩXi(Xj−Xi)⊗CjWijVj]−1(34)

Laplacian The definition of the Laplacian can be written as

∇2φ(X,t)=∇⋅∇φ(X,t)(35)

In order to obtain the Laplacian relationship between the center point X and its neighboring point X_, 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 the (RC)XX_ as 0, which yields

∇φ(X,t)XX_=[(φ(X_,t)−φ(X,t))⊗C][(X_−X)⊗C]−1(36)

However, the rank of matrix [(X_−X)⊗C] is 1, which is not available for the calculation of the inverse when the gradient is carried over from the vector gradient calculation. To circumvent the issue, we introduce A=(X_−X)⊗C. Then it can be replaced by its average value within an infinitesimal domain via exploiting Eq. (18),

AXX_avg=1VΩXlimVΩX→0∫ΩX[(X_−X)⊗C]dV(37)

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

∇φ(X,t)XX_=[(φ(X_,t)−φ(X,t))⊗C](AXX_avg)−1(38)

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

(∇⋅φ(X,t))XX_VX_=φ(X_,t)⋅nXX_SX_(39)

where S_=V_||X_−X||. As will be explained in the later section, we can also replace nXX_ with C. Substitution of S_ into Eq. (39) yields

(∇⋅φ(X,t))XX_||X_−X||=φ(X_,t)⋅C(40)

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

(∇⋅∇φ(X,t))XX_||X_−X||=[(φ(X_,t)−φ(X,t))⊗C](AXX_avg)−1⋅C(41)

Hence the following Laplacian formulations can be obtained as

∇2φ(X,t)=∇⋅∇φ(X,t)=limVΩX→01VΩX∫ΩX[(φ(X_,t)−φ(X,t))⊗C]AXX_avg−1⋅C||X_−X||dV(42)

The discretized formulations is listed as following:

∇2φi=1VΩXi∑j∈ΩXi[(φj−φi)⊗Cj]Aijavg−1⋅Cj||Xj−Xi||WijVj(43)

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,

∇⋅σ=Fσ=D∇u(44)

where D is elastic matrix as following

D=E1−ν2[1ν0ν10001−ν2]for plane stress(45)

D=E(1+ν)(1−2ν)[1−νν0ν1−ν00012ν]for plane strain(46)

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

∇ui=[∑j∈ΩXi(uj−ui)⊗CijWijVj][∑j∈ΩXi(Xj−Xi)⊗CijWijVj]−1⏟[A](47)

where [A] is a two by two matrix that is computed first. The summation is applied for every particle within particle i’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,

∇u1=[(u2x−u1x)(u2y−u1y)][(C12x)(C12y)][A11A12A21A22]W12V2(48)

=[C12x(u2x−u1x)C12y(u2x−u1x)C12x(u2y−u1y)C12y(u2y−u1y)][A11A12A21A22]W12V2(49)

Thus,

[∂u12x∂x∂u12x∂y∂u12y∂x∂u12y∂y]=[(A11C12x+A21C12y)(u2x−u1x)(A12C12x+A22C12y)(u2x−u1x)(A11C12x+A21C12y)(u2y−u1y)(A12C12x+A22C12y)(u2y−u1y)]W12V2(50)

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.

[∂u12x∂x∂u12x∂y∂u12y∂x∂u12y∂y]=W12V2[(A11C12x+A21C12y)(u2x−u1x)(A12C12x+A22C12y)(u2x−u1x)(A11C12x+A21C12y)(u2y−u1y)(A12C12x+A22C12y)(u2y−u1y)]=W12V2[−(A11C12x+A21C12y)0(A11C12x+A21C12y)0−(A12C12x+A22C12y)0(A12C12x+A22C12y)00−(A11C12x+A21C12y)0(A11C12x+A21C12y)0−(A12C12x+A22C12y)0(A12C12x+A22C12y)]×[u1xu1yu2xu2y](51)

We can also get the engineering strain as follows:

[∂u12x∂x∂u12x∂y12(∂u12x∂y+∂u12y∂x)]=W12V2[−(A11C12x+A21C12y)0(A11C12x+A21C12y)00−(A12C12x+A22C12y)0(A12C12x+A22C12y)−(A12C12x+A22C12y)2−(A11C12x+A21C12y)2(A12C12x+A22C12y)2−(A11C12x+A12C12y)2]⏟[∇1]:gradient operator for particle 1[u1xu1yu2xu2y]⏟[u](52)

Eq. (52) can be simply written as

[ϵ1]=[∇1][u](53)

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

[σ1]=[D][∇1][u](54)

It is worth noting that [D][∇][u] includes all the components in the particle i’s influence domain.

[σ]=[D][∇][u]=[[D][∇1][D][∇2]⋮[D][∇n]][u](55)

where [∇1] and [∇2] represent the gradient operator obtained from the first and the second particle as the i th 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.

∇⋅σ1=[C12xC12y]⋅[A11A12A21A22][σ112−σ111σ122−σ121σ212−σ211σ222−σ221]W12V2(56)

=[A11C12x+A21C12yA12C12x+A22C12y][σ112−σ111σ122−σ121σ212−σ211σ222−σ221]W12V2(57)

=[B1(σ112−σ111)B2(σ122−σ121)B1(σ212−σ211)B2(σ222−σ221)](58)

where B1=(A11C12x+A21C12y) and B2=(A12C12x+A22C12y) for simplicity. Here we further isolate σ together.

[∇⋅σ1x∇⋅σ1y]=W12V2[−B1−B200B1B20000−B1−B200B1B2][σ111σ121σ211σ221σ112σ122σ212σ222](59)

or for the engineering stress,

[∇⋅σ1x∇⋅σ1y]=W12V2[−B10−B2B10B20−B2−B10B2B1]⏟[∇⋅]:divergence operator for particle 1[σxx1σyy1τxy1σxx2σyy2τxy2](60)

Thus, we can get the stiffness matrix for the particle 1, by combining Eqs. (55) and (60),

[∇⋅][D][∇][u]=FW12V2[−B10−B2B10B20−B2−B10B2B1][[D][∇1][D][∇2]]⏟[K1][[u]]=[F1xF1y](61)

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 constraints g(u) = 0. The action integral is given as

S(u,λ∗)=∫t0tfL∂t+λ∗Tg(u)(62)

and we can take the variation of the action integral,

δS(u,λ∗)=0→∂∂t(∂L(u,u.)∂u.)−∂L(u,u.)∂u+GT(u)λ∗=0(63)

g(u)=0(64)

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

Qiner(u,u.,u..)+Qint(u)=Qappl(u,u.,t)+GT(u)λ∗(65)

g(u)=0(66)

where, Qiner, Qint and Qappl are 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 into ndom number of subdomains, Ωi (for i=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 Ωi is, (for i=1,2,…,ndom),

Miu..i=Qiint(ui)+Qiappl(ui,u.i,t)+GiTλ∗(t)(67a)

0=∑i=1ndomGiu.i(67b)

where ui(t), u.i(t):=dui/dt, and u..i(t):=d2ui/dt2 are the displacement, velocity, and acceleration vectors of the ith-subdomain, respectively; Mi is the mass matrix of the ith-subdomain; and Qiint(ui) and Qiappl(ui,u.i,t) are the internal and applied force of the ith-subdomain, respectively.

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

The matrix Gi is 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 Δt and Δti be the system time step size and the subdomain time step size for Ωi, respectively, and define the ratio as τi=Δt/Δti≥1 for i=1,2,…,ndom. The system and subdomain time step sizes are defined as

Δt=tn+1−tn(68a)

Δti=tn+m+1τi−tn+mτi(68b)

where n denotes the system time step, and m represents 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:

Mia~i(n,m)=Qiint(u~i(n,m))+Qiappl(u~i(n,m),v~i(n,m),tn+m+W1τi)+GiTλ~∗(n,m)(69a)