A Finite Element Flux-corrected Transport Method for Wave Propagation in Heterogeneous Solids

When moving discontinuities in solids need to be simulated, standard finite element (FE) procedures usually attain low accuracy because of spurious oscillations appearing behind the discontinuity fronts. To assure an accurate tracking of traveling stress waves in heterogeneous media, we propose here a flux-corrected transport (FCT) technique for structured as well as unstructured space discretizations. The FCT technique consists of post-processing the FE velocity field via diffusive/antidiffusive fluxes, which rely upon an algorithmic length-scale parameter. To study the behavior of heterogeneous bodies featuring compliant inter-phases of any shape, a general scheme for computing diffusive/antidiffusive fluxes close to phase boundaries is proposed too. The performance of the new FE-FCT method is assessed through one-dimensional and two-dimensional simulations of dilatational stress waves propagating along homogeneous and composite rods.


Introduction
The propagation of waves in elastic solids is governed by a second-order, hyperbolic differential equation. Waves therefore travel inside bodies with a finite speed, and cause an abrupt change of the local velocity and stress fields across their fronts.
Numerical methods have to accurately track such moving inner discontinuities. Since standard dis-placement-based finite element (FE) schemes adopt continuous interpolation fields to mimic the discontinuous ones inside the modeled domain, non-physical high-frequency oscillations show up around wave fronts; these oscillations are a numerical artifact and need to be filtered out of the solution. Algorithmic treatments of this issue have been proposed in the literature, and they typically consist in artificial viscosity or mesh adaption. Focusing on time integration, algorithms like the α−method [1] or the generalized α−method [2] were also devised to damp oscillations. All these methods can reduce in size such spurious effects, but sometimes entail energy dissipation which is, again, non-physical. An alternative treatment, whose roots can be traced back to the seminal work of Boris & Book [3,4], is the flux-corrected transport (FCT) method. The FCT algorithm consists of post-processing a standard FE solution: diffusive and antidiffusive fluxes, the latter being appropriately limited in size, are handled to improve the discrete velocity field around the discontinuities, and to filter out spurious oscillations. This method has been extensively used to simulate the propagation of shock waves in fluids [5,6]; recently, it has been adopted to simulate traveling stress waves in solids [7]. Noteworthy results have been obtained in [8] through a coupling of a displacement-based FE solution and the FCT algorithm; owing to the adopted structured meshes, results for bodies of arbitrary shapes were based on a partitionof-unity enrichment of the nodal shape functions.
In this work, to study the dynamics of heterogeneous bodies we propose two enhancements to the frame developed in [8,9]. First, to simulate the propagation of stress waves inside domains of arbitrary shape, an algorithmic length-scale ℓ FCT is introduced: this length-scale allows to define local supports of finite size, which are independent of the space discretization, wherein diffusive/antidiffusive fluxes are computed. Second, the rationale behind the computation of diffusive/antidiffusive fluxes close to body or phase boundaries is revisited, so as to permit the treatment of compliant interphases confined along loci of zero measure (i.e. surfaces in three-dimensional frames and lines in two-dimensional frames).
These two enhancements are of paramount importance when dynamic failure of quasi-brittle polycrystals (like, e.g., polysilicon) needs to be modeled, since damaging phenomena at the micro-scale are incepted as soon as the tensile strength is locally attained. Because of the polycrystal micro-structure, traveling waves are partially reflected by each grain boundary, and eventually lead to complex stress patterns. If the numerical solution does contain the aforementioned fictitious oscillations, the amplitude of the local stress field may be artificially increased and, therefore, damage wrongly started. Through simulations of stress waves traveling along homogeneous as well as heterogeneous (bimaterial) rods, we show that the proposed FE-FCT method can be used to accurately study the evolution of the stress field in quasi-brittle polycrystals.
As far as notation is concerned, a matrix one will be adopted throughout the whole paper with uppercase and lowercase bold symbols respectively denoting matrices and vectors, a superscript T standing for transpose, and a superposed dot representing time rates.

Governing relations
Let Ω be a heterogeneous three-dimensional body; its smooth outer boundary, with unit outward normal m, be constituted by the two disjoint sets Γ u and Γ τ , where displacements and tractions are respectively assigned. Without any loss of generality, let us assume that Ω is a continuum made of two phases (Ω + and Ω − , with Ω = Ω + ∪ Ω − ), tied together in the initial configuration along the flat interface Γ, see Figure 1. Γ actually represents the interphase between Ω + and Ω − ; since the thickness of this interphase is usually at least one order of magnitude smaller than the characteristic size of Ω, the interphase itself is modeled as a zero-thickness interface. Damaging processes along Γ may cause opening (mode I) and/or sliding (mode II and mode III) displacement discontinuities in the n direction and in the s 1 − s 2 plane, respectively.
The equilibrium of Ω at time t is governed by: where Γ + and Γ − are the sides of Γ respectively belonging to Ω + and Ω − , and according to Voigt's notation [10]: σ is the stress vector, which gathers the independent components of the stress tensor;b andτ are the assigned loads in the bulk Ω\Γ and along Γ τ ; ̺ is the mass density of the bulk material; u is the acceleration in Ω\Γ; C is the differential compatibility operator; M and N are the matrices collecting the components of the unit vectors m and n.
In the small strain regime, compatibility reads: where: ε is the strain vector; u is the displacement in Ω; [u] is the displacement discontinuity along Γ; u is the assigned displacement along Γ u .
The body is conceived to be initially at rest, i.e.: As far as constitutive modeling is concerned, the bulk is assumed to behave elastically according to: where E Ω is the bulk elasticity matrix. Along Γ, damage can be incepted once a local strength criterion is satisfied; to simplify matters we assume that the stress waves do not cause any dissipative phenomena. The interface thus behaves elastically too, according to: E Γ being the interface elasticity matrix.

Space-time discretized solution
The weak form of the equilibrium equations (1)-(4), allowing for bulk and interface constitutive laws (9)-(10), reads [11][12][13]: where: v is the test function, and ε v = Cv; U is the trial solution space, featuring displacement fields u continuous in Ω\Γ, possibly discontinuous along Γ, and fulfilling the boundary condition (7) on Γ u ; U 0 is the relevant variation space. In view of the assumed linearized kinematics, Γ ≡ Γ + ≡ Γ − holds for the interface. Now, let the finite element approximation of the displacement, displacement jump and deformation fields be: where: u h is the nodal displacement vector; Φ is the matrix of nodal shape functions; B Γ and B Ω are, respectively, the interface and bulk compatibility matrices. The semi-discretized equations of motion of the solid turn out to be: where the mass matrix M, the bulk and interface stiffness matrices K Ω and K Γ , and the external load vector F are: To advance the solution of (15) in time we employ an explicit Newmark method. Having partitioned the time interval of interest according to a predictor-integrator-corrector splitting is followed according to: • corrector: To resemble the average acceleration scheme, which is unconditionally stable and second-order accurate, β = 0.25 and γ = 0.5 are adopted in the preceding equations. The time step size ∆t has been always set so as to fulfill the Courant condition in the bulk Ω\Γ. Moreover, to speed up the explicit integrator phase (22), the mass matrix M has been diagonalized via row-sum lumping [14].
Customary FE methods finally assumeu i+1 =u i+1 , without any treatment to avoid the occurrence of spurious oscillations. With the FCT scheme, instead, a post-processing stage is added to deal with such an issue; details on this phase are furnished in Section 3.

Flux-corrected transport method
We already pointed out that the accuracy of customary FE simulations of moving discontinuities, as obtained via the procedure outlined in Section 2.2, is usually spoiled by oscillations first arising around the discontinuity fronts, where both the velocity and stress fields are discontinuous.
To avoid such kind of fictitious oscillations, the FCT algorithm has been adopted in [8,15,16]: it amounts to post-processing the local velocity field, which actually suffers discontinuities, via viscouslike diffusive and antidiffusive fluxes, the latter ones being appropriately limited in size. For ease of notation, we detail here the algorithm in the case of a one-dimensional domain Ω; pointing to the j−th node of the space discretization and assuming the mesh to be structured, the standard version of the FCT algorithm reads: • compute diffusive fluxes: • update the solution through diffusive fluxes: • compute anti-diffusive fluxes: • apply limitation to antidiffusive fluxes: Figure 3. Sketch of the shape of nodal supports, either when they are wholly included in Ω or when they intersect the boundary of any phase.
In the above equations: termsu j−1 andu j+1 respectively denote the velocities at the nodes on the left and right sides of the j−th one, see Figure 2(a); η D and η A are the diffusive and antidiffusive coefficients. In two-and three-dimensional settings the above procedure can be generalized by following the same path to correct the velocity component aligned with each reference axis.
As far as the computational burden of the proposed FE-FCT scheme is concerned, it can said that: the explicit Newmark time integration consists of the five loops (20)-(24) over the nodes; the FCT algorithm consists of the eight loops (25)-(32) over the nodes, and of some additional operations in (30)-(31) (i.e. min, max, positive part). Hence, in terms of computing time, simulations run with the present FE-FCT scheme turn out to be at least 13 5 = 2.6 times longer than the corresponding customary FE ones. When discontinuities travel inside irregular domains, ad-hoc procedures have been followed to fulfill boundary conditions or to model the response of multi-phase materials, see e.g. [8]. In this work we offer a simple procedure to allow also the adoption of unstructured meshes: in the above equations, instead of using nodal information to update the velocity field, we handle information on a local grid centered at the j−th node and characterized by a constant spacing among grid points, see Figure 2(b). This spacing between grid points, termed ℓ FCT , is not to be meant as the distance between neighboring nodes, but instead an algorithmic parameter to be tuned. The velocity at each point of this grid is obtained from the nodal ones through interpolation, using the shape functions gathered in Φ (see Equation 12).
A proof of the conservation of local momenta, or error estimators for the conservation laws are beyond the goals of this paper, and are left for future investigations.
To deal with compliant interfaces, whose behavior is described by constitutive laws like (10), a specific treatment of phase boundaries is needed. Here we propose the following rationale: during the FCT correction step, each phase is assumed to be an independent body. Therefore, if a local support intersects either the outer boundary Γ u ∪ Γ τ or the interface Γ, it is re-shaped according to what sketched in Figure 3. Since in Equations (25)-(32) the (j − 2)−th, (j − 1)−th, (j + 1)−th and (j + 2)−th terms are dropped when the relevant points (or nodes in the standard formulation) fall outside Ω, the algorithm is enhanced by adopting the same criterion if grid points fall outside the phase the j−th node belongs to. If this procedure is not followed, the FCT correction step can not feel interfaces endowed with their own constitutive laws. Forthcoming results will show that the two enhancements here proposed permit to increase much the accuracy of FE results, specially close to stress wave fronts, by preventing macroscopic fluctuations to show up in the velocity and stress fields, even when unstructured space discretizations are adopted.

Longitudinal waves in homogeneous and heterogeneous rods
To assess the capability of the newly proposed FE-FCT method, the propagation of (longitudinal) waves in homogeneous as well as heterogeneous slender bodies (rods) is simulated. To check the accuracy of the method, results are compared to an analytical solution derived next.

Analytical solution
Let Ω be a slender body, such that only the propagation of longitudinal waves is of interest. According to the schematic of Figure 4, the body is assumed to be loaded by a constant pressureτ , acting over the time interval ∆t τ . To get insights into the effect of inner boundaries on wave reflection and dispersion, a two-phase body in considered, with phases A and B (of length L in the direction of wave propagation) respectively featuring Young's moduli E A and E B , and mass densities ̺ A and ̺ B ; the relevant longitudinal wave speeds thus are c A = E A ̺ A and c B = E B ̺ B . In view of the elastic response of the bulk, a tensile wave propagating inside material A turns out to be characterized by the following relationship between the stress amplitude σ in and the particle velocity v in : Figure 7. Wave propagating along a homogeneous rod discretized with an unstructured mesh: stress field at t = 1.8 · 10 −2 s. Effect of parameters η D = η A on the capability of the newly proposed FE-FCT scheme to accurately track the stress discontinuity across the wave fronts.
(a) (b) After impinging upon the interface Γ at t =t, the incident wave induces reflected and transmitted waves, respectively characterized by relations [17]: where σ − and σ + are the amplitudes of the reflected and transmitted stress waves, and v − and v + are the relevant particle velocities. Equilibrium across Γ requires: while compatibility reads: [u] being the opening velocity on Γ. Account taken of the interface constitutive law (10), written for mode I loadings as τ ≡ σ + = E Γ [u], one arrives at the differential equation: where: t c = ̺ A c A 2E Γ is the time-scale of the opening process along interface Γ (sometimes called characteristic relaxation time [17]); γ = ̺ B c B ̺ A c A . Upon integration of (38), the solution in terms of time histories In case of perfect adhesion between the two phases, i.e. in case of E Γ → ∞, the exponential terms in (39) and (40) drop to zero. The time-dependent feature of the solution is thus lost in view of the non-evolving microstructure along the interface.

FE-FCT results
Similarly to what done in [9], all the results to follow have been obtained by adopting: L = 0.5 m; ̺ A = ̺ B = 10 2 kg/m 3 ;τ = 10 4 Pa, so as to initially cause the propagation of a tensile wave inside material A of magnitude σ in = 10 4 Pa; ∆t τ = 4 · 10 −3 s. As for the FCT parameters, if not otherwise stated we have adopted: η D = η A = 0.125 [9], and ℓ FCT = L 1000 . In a first example the rod has been assumed homogeneous, i.e. made of two phases featuring the   FCT scheme when either structured or unstructured space discretizations are adopted. In our analysis, unstructured one-dimensional meshes have been obtained from the relevant structured ones as follows.
Starting with a nodal spacing ∆x (here ∆x = 5 · 10 −4 m, meaning that 1000 quadratic elements, with characteristic size ℓ 1D e = 2∆x = L 500 , have been adopted to discretize the whole specimen), a random sequence ran[−0.5 0.5] is generated so as: where: x j struct and x j unstruct are the longitudinal coordinates of the j−th node in the structured and unstructured meshes, respectively; χ = 0.6 is a scaling factor, introduced to ensure positivity of the length of each element in the unstructured mesh. Figure 5 shows the stress distribution in the rod at time instants t = 0.5 · 10 −2 s and t = 1.8 · 10 −2 s; it can be seen that the tensile wave, after reflection at the rear free surface of the specimen, becomes a compressive wave of same amplitude. These results have been obtained with a structured mesh, hence the standard [8] and the newly proposed FE-FCT schemes furnish the same outcomes. Moreover, the analytical solution can not be distinguished from the FE-FCT ones; on the other hand, ripples show up in the FE solution, initially just behind the wave fronts (see Figure 5(a), where the discontinuities are propagating from left to right) and eventually spoiling the accuracy of the solution along the whole specimen (see Figure 5(b)).
In case of unstructured meshing, the methods here compared furnish the results summarized in Figure  6 in terms of stress field at t = 1.8 · 10 −2 s. As expected, the FE solution wildly oscillates, with ripples of increased amplitude as compared to the structured mesh case; the standard FE-FCT scheme is no more able to fully damp oscillations, whereas the new one works properly and prevents the birth of spurious maxima/minima. It must be noticed that the small delay of the new FE-FCT method in the attainment of the analytical solution (characterized by sharp discontinuities), as shown in the two insets of Figure 6, may be avoided or reduced in duration by a fine tuning of Newmark parameters β, γ and FCT parameters η D , η A . To prove it, Figure 7 shown the effect of η D = η A on the aforementioned delay; as far as this feature of the solution is concerned, it can be seen that smaller parameter values strongly reduce the delay, even if the value of the longitudinal stress may get affected. Only results relevant to η D = η A are here reported; on the one hand, if η D > η A the solution is characterized by an excessive diffusion around the discontinuities, with fronts smeared over too many FEs; on the other hand, if η D < η A the solution has a tendency to instability, characterized by ripple amplitudes continuously increasing in time.
For comparison purposes, outcomes of a two-dimensional simulation are depicted in Figure 8 in terms of contour plots of the longitudinal stress before (Figure 8(a)) and after (Figure 8(b)) reflection at the rear free surface. In this two-dimensional case a much coarser space discretization has been adopted, with 6-node triangular elements featuring a characteristic size ℓ 2D e = L 100 . It can be seen that the stress wave maintain sharp fronts, always perpendicular to the propagation direction independently of local mesh features. Because of the coarse discretization, it is expected that spurious oscillations are here characterized by an increased amplitude when compared to the one-dimensional case. Figure 9, where the time history of the longitudinal stress at x = 0.5 m is reported, testifies the performance of the proposed FE-FCT scheme: even though the ratio between the characteristic sizes of FEs is ℓ 1D e ℓ 2D e = 5, the discrepancy between the two solutions amounts at most to 2%. In a second example, to start assessing the efficiency of the FE-FCT scheme in the presence of inner boundaries, it has been assumed that E B = 2 · 10 6 Pa at fixed values of all the other material parameters; perfect adherence between the two phases was also assigned. FE and FE-FCT results are shown in Figure  10 in terms of stress field at t = 0.5 · 10 −2 s and t = 1.3 · 10 −2 s. At t = 0.5 · 10 −2 s (Figure 10(a)), since the leading front of the stress wave has not impinged yet upon the interface (placed at x = 0.5 m) the results are similar to those of Figure 5(a), showing again that the FCT scheme properly damps oscillations (as before, the analytical solution is not reported here since it can not be distinguished from the FE-FCT one). Moreover, when waves are partially reflected by the inner surface the solution still proves accurate, even without explicitly accounting for the acoustic impedance of each phase as done in [9].
In a third example, the behavior of a homogeneous two-phase specimen (E A = E B = 10 6 Pa) containing a compliant interface has been studied. In the case E Γ = 5 · 10 7 N/m 3 , results in terms of the stress field at t = 0.9 · 10 −2 s are depicted in Figure 11; even in the presence of a compliant interface, it can be seen that the new FE-FCT method correctly reproduces the response of the body without spuri-ous oscillations. To further assess the accuracy of the method, the time history of the normal traction τ across Γ is shown in Figure 12. To increase the duration of the transitory stage and to get insights into the performance of the FE-FCT method, the interface stiffness has been decreased to E Γ = 5 · 10 6 N/m 3 , and the loading interval ∆t τ has been slightly increased accordingly, so as to ensure completion of interface opening till [u] ⋆ =τ E Γ . Once again, the FE solution looks oscillatory; this behavior becomes a serious issue in case of nonlinear, irreversible interface modes, since spurious unloadings and overshoots may cause a drift of the computed solution away from the exact one. The FE-FCT solution instead presents only an oscillation just after the leading front of the stress wave has impinged upon the interface; afterwards, the FCT scheme filters out possible subsequent ripples. This is also shown in Figure 13, where results of one-dimensional and two-dimensional simulations are compared; in the two-dimensional case, the same space discretization depicted in Figure 8 has been adopted. According to previous results, the higher oscillation amplitude around t = 0.005 s in the two-dimensional simulation is caused by the coarser mesh; subsequently, the FCT scheme completely damps ripples.
The small time discrepancy (amounting to ∆t ⋆ ∼ = 10 −5 s in the one-dimensional case) between FE-FCT and exact solutions is due to the dispersive properties of the FE solution and can be explained as follows. While the analytical solution is characterized by sharp wave fronts, in the FE-FCT one discontinuities are smeared over the resolved length-scale, which is typically one element wide; ∆t ⋆ is about the time interval required by the stress wave to travel across ∆x (or the corresponding node spacing in the unstructured mesh case) at speed c A = c B .

Concluding remarks
In this paper we have proposed a finite element flux-corrected transport (FE-FCT) method for the simulation of stress wave propagation in homogeneous as well as heterogeneous bodies.
Through a comparative assessment of existing remedies for the oscillatory behavior of standard FE simulations of moving discontinuities, the FCT algorithm emerged as a very promising tool. It consists of post-processing the velocity field, as furnished by a customary FE scheme, via diffusive and antidiffusive fluxes, the latter being appropriately limited in size. The main limitations of present implementations, already foreseen e.g. in [8,9], are the need of structured space discretizations, which entails an adhoc treatment of irregular boundaries of the body, and the former knowledge of the acoustic impedance of each phase in case of waves traveling inside heterogeneous bodies. This second issue prevents the simulation of the dynamics of polycrystals (i.e. highly heterogeneous solids) in the presence of defective, compliant interphase boundaries, since acoustic impedances turn out to be influenced by the interphase compliances and by the angle between the wave propagation direction and each phase boundary. We have here proposed a new FCT scheme, which is able to deal with irregular and compliant phase boundaries; this new scheme relies upon an algorithmic length-scale, used to define the size of the nodal support within which the velocity field is corrected through the diffusive/antidiffusive fluxes. Close to phase boundaries, alike close to boundaries of the whole body, the aforementioned local support is properly re-shaped: this permits to account for possible interphase compliances.
Through one-dimensional and two-dimensional simulations of stress waves propagating along slender bodies (rods), it has been shown that the newly proposed FE-FCT method is effective and robust. Results, never showing macroscopic fluctuations in the stress and velocity fields, turned out to be very accurate when compared to available analytical solutions, and independent of the (unstructured) space discretizations adopted.