Model-Free Data-Driven Viscoelasticity in the Frequency Domain

We develop a Data-Driven framework for the simulation of wave propagation in viscoelastic solids directly from dynamic testing material data, including data from Dynamic Mechanical Analysis (DMA), nano-indentation, Dynamic Shear Testing (DST) and Magnetic Resonance Elastography (MRE), without the need for regression or material modeling. The problem is formulated in the frequency domain and the method of solution seeks to minimize a distance between physically admissible histories of stress and strain, in the sense of compatibility and equilibrium, and the material data. We metrize the space of histories by means of the flat-norm of their Fourier transform, which allows consideration of infinite wave trains such as harmonic functions. Another significant advantage of the flat norm is that it allows the response of the system at one frequency to be inferred from data at nearby frequencies. We demonstrate and verify the approach by means of two test cases, a polymeric truss structure characterized by DMA data and a 3D soft gel sample characterized by MRE data. The examples demonstrate the ease of implementation of the Data-Driven scheme within conventional commercial codes and its robust convergence properties, both with respect to the solver and the data.


Introduction
A number of techniques are presently available for characterizing viscoelastic materials in the frequency domain, including Dynamic Mechanical Analysis (DMA) [1], nanoindentation [2,3], dynamic shear testing (DST) [4,5] and magnetic resonance elastography (MRE) [6], among others. The systematic application of these techniques results in large material data sets that challenge analytical representation and computational utilization. Two main strategies arise in response to this wealth of material data: i) Modeling the data using conventional rheological materials or, more recently, neuralnetwork representations and machine learning regression [7], then using the calibrated material models in the solution of boundary-value problems; ii) Using the data directly in the solution of suitably-formulated boundaryvalue problems in the spirit of Data-Driven (DD) mechanics [8], thereby bypassing the material-modeling step altogether. The main advantage of the DD paradigm is that is makes direct use of the data, all the data, and nothing but the data. In particular, no loss of information is incurred with respect to the material data set, no modeling biases or errors are introduced and no ad hoc choices need to be made regarding the functional form of the material model. In addition, if material behavior is properly sampled with increasing fidelity, the DD solution may be expected to converge to the solution of the underlying, and unknown, material response [9].
The application of the DD paradigm to history-dependent materials has been investigated in [10]. The history dependence of the material response requires characterizing entire state histories ( (x, t), σ(x, t)) of stress and strain at material points x in the solid. The collection of all such state histories defines a material-history data set D that characterizes the material response of the solid. In addition, the state histories are subject to the constraints of compatibility and equilibrium, both in the interior and at the boundary. The collection of all histories satisfying those physical constraints defines a constraint set E of admissible state histories. Evidently, the actual histories undergone by the solid are those that lie at the intersection of D and E, i. e., histories that are physically admissible, in the sense of compatibility and equilibrium, and are contained in the material-history data set characterizing material behavior.
In this work, we consider the case, often encountered in practice, in which the full material behavior of a viscoelastic material is not known, but only partially characterized by empirical material-history data sets D h . For instance, such data sets could be the result of DMA [1], in which case the sampled histories ( (x, t), σ(x, t)) are harmonic in time, or could be the result of creep and relaxation tests, or computed from micromechanics [11,12], or determined otherwise. In general, the intersection of the empirical materialhistory data set D h and the constraint set E is likely to be empty, i. e., no material history recorded in D h is likely to be compatible and in equilibrium with the applied loading history. In order to circumvent this difficulty, we relax the notion of solution and seek to determine admissible state histories in E that are closest, in the sense of a suitable metric, to the empirical material-history data sets D h . In this manner, we expect to identify approximating state histories ( h (x, t), σ h (x, t)) that converge to the exact state history ( (x, t), σ(x, t)) as the sequence D h of sets converges to D in some appropriate manner.
Here, we focus on questions of practical implementation of DD linear viscoelasticity and demonstrate convergence with respect to material data by means of selected numerical tests. A first question of crucial importance concerns the choice of metric used to select approximate state histories. We specifically consider empirical material-history data sets D h generated by frequency-domain characterization techniques such as DMA and MRE. In particular, the set D h contains uniform harmonic state histories ( (t), σ(t)) defined for all times. Evidently, integral norms fail to be defined for such histories and cannot be taken as a basis for DD analysis. Instead, we note that the Fourier transforms (ˆ (ω),σ(ω)) of the sampled state histories are Dirac-deltas concentrated at the applied frequency ω, with complex amplitudes related by a complex modulus E. A natural empirical representation of such behavior is by means of material data sets D h consisting of pairs (ω i , E i ) of applied frequencies ω i and corresponding complex moduli E i . In addition, a canonical metric for computing the distance between two such data points is supplied by the flat norm [13,14], which measures the discrepancy in both amplitude and frequency in the Fourier representation. A crucial property of the flat norm is that it allows the response of the material at one frequency to be compared to data at nearby frequencies. This property is particularly useful when dealing with sparse material data sets D h containing data at a finite, possibly small, number of applied frequencies.
The article is organized as follows. Following an introductory section aimed at defining the problem and establishing notation, Section 2, an implementation of the frequency-domain DD linear-viscoelastic formalism is presented in Section 3. In Section 4, we illustrate the DD approach by means of two examples of application: i) A truss structure made of a polymer characterized by DMA data, representative of polymer-based additive manufacturing technology [15]; and ii) a three-dimensional viscoelastic gel characterized by MRE data, commonly used as phantom in the field of ultrasonic biomechanics [16]. In both cases, we demonstrate the convergence of the DD approach with respect to the data. We conclude with a summary and outlook in Section 5.

Dynamic viscoelasticity
For definiteness, we restrict attention to finite-dimensional discrete models. The models under consideration comprise m material points, or structural members in the setting of structural mechanics, whose state at time t is characterized by strains (t) ≡ { e (t) ∈ R d , e = 1, . . . , m} and stresses σ(t) ≡ {σ e (t) ∈ R d , e = 1, . . . , m}, where d is the stress and strain dimension.
2.1. Governing equations. The state of the system is subject to compatibility and dynamic equilibrium constraints of the form e (t) = B e u(t) + g e (t), e = 1, . . . m, (1a) where u(t) ∈ R n is the array of unconstrained nodal displacements at time t, with n the spatial dimension, f (t) ∈ R n is a nodal force array at time t, g e (t) ∈ R d is an initial strain at time t resulting from prescribed displacement boundary conditions, w e are positive weights, B e ∈ L(R n , R d ) is the straindisplacement matrix for material point e and M ∈ R n×n is a mass matrix. The constraints (1) can also be expressed in compact form as with B = (B 1 , . . . , B m ) ∈ R N ×n , N = md, (t) = ( 1 (t), . . . , m (t)) ∈ R N , g(t) = (g 1 (t), . . . , g m (t)) ∈ R N , σ(t) = (σ 1 (t), . . . , σ m (t)) ∈ R N and W ∈ R N ×N is a diagonal matrix of weights such that W σ = {w 1 σ 1 , . . . , w m σ m }.
Classically, problem (1) is closed by assuming a hereditary dependence of the local stress history σ e (t) on the local strain history e (t) at every material point e. The mapping F e maps local histories of strain to local histories of stress and must be causal, i. e., σ e (t) = σ e (t) if σ e = F e ( e ), σ e = F e ( e ), and e (s) = e (s) for all s ≤ t.

Fourier representation.
Dynamic viscoelasticity techniques, such as Dynamic Mechanical Analysis (DMA) [1], are common means of characterizing viscoelastic behavior as a function of frequency and temperature. DMA consists of applying a harmonic force to a test specimen and measuring the resulting strain response. Evidently, DMA contemplates histories of strain and stress that extend over infinite time. In order to account for such histories, we adopt a Fourier representation. Formally, an application of the Fourier transform to (2) giveŝ whereˆdenotes the Fourier transform of a function and ω is the frequency. We note that the Fourier transforms of infinite wave trains can be general measures. For instance, the Fourier transform of a harmonic time history is a Dirac measure centered at the harmonic frequency. Therefore, in general the relations (3) must be interpreted in the sense of measures.
, J e (t) to 0 for t < 0, is the relaxation function and * is the convolution operator. We may regard (4) as a mapping that gives the stress σ(t) and time t as a function of the past history of strain. Evidently, (4) satisfies causality by construction.
Taking the Fourier transform of (4) and applying the convolution theorem, we obtain are the storage and loss moduli, respectively. We note that the relation (6) separates by frequency in the linear case. These moduli satisfy the Kramers-Kronig relations which are necessary and sufficient for causality.

2.4.
Steady state under monochromatic transduction. As a particular example of application, suppose that the forcing is harmonic, where F ∈ C n and G ∈ C N are constant complex amplitudes and Ω is the transduction frequency. Harmonic forcing of this type is induced, e. g., as a result of monochromatic transduction. Then, we may solve (3) by setting with U ∈ C n , E ∈ C N and S ∈ C N also constant complex amplitudes such that Suppose that the complex modulus E(Ω) of the material is known exactly. Then, and (11) reduces to which is a monochromatic form of (8) and can be solved directly for the steady-state displacement amplitude.

Data-Driven reformulation
An essential difficulty with the classical approach just described is that, whereas the field equations (3) are known exactly and are uncertainty-free, the material law (5) is known only empirically through experimental data. We wish to reformulate the problem in such a way that predictions can be made directly from the field equations and the available empirical data. We additionally wish such Data-Driven (DD) predictions to be approximating in the sense of converging to the classical solution when the empirical data converges, in some appropriate sense to be made precise, to an underlyingand unknown-material law of the form (5).
Likewise, we may regard (2), or its Fourier representation (3) as defining a constraint set of admissible histories The aim then is to find all histories ( (t), σ(t)) in the intersection D ∩ E, i. e., all histories that are admissible and satisfy the constitutive relations of the material.

Empirical data.
Suppose that the complex modulus E(ω) of a linear viscoelastic material, and therefore its data set D, eq. (14), is not known exactly but only approximately through empirical data. In the present work, we specifically consider local material data of the Dynamical Material Analysis (DMA) type. Thus, the behavior of every material point e is characterized by a set P = {(ω i , E i ), i = 1, . . . , N }, consisting of measured pairs of frequencies and corresponding complex moduli. Alternatively, we may regard the data as a collection of experimentally observed harmonic histories of stress and strain (16) e (t) = Ee iω i t , σ e (t) = Se iω i t , S = E i E, i = 1, . . . , N. Collecting these histories, we obtain a local material data set of the form where Z loc denotes the phase space of local histories ( e (t), σ e (t)) of material point e. The corresponding global material data set is, then, , σ e (t)) ∈ D loc , e = 1, . . . , m}.

3.3.
The data-driven problem with empirical material data. When the material behavior is characterized by means of an empirical data set D such as just described, the intersection D ∩ E is likely to be empty, i. e., there is no history ( (t), σ(t)) that is both admissible and in the material data set, and the notion of solution needs to be extended. Following [10], we seek instead a pair of histories, an admissible history ( (t), σ (t)) ∈ E and a material history ( (t), σ (t)) ∈ D, that are as close as possible to each other.
A suitable notion of distance between dynamical histories, including infinite wave trains whose Fourier transforms are general measures, is supplied by the flat norm [13]. For present purposes, it suffices to characterize the flat-norm distance between harmonic histories. To this end, we begin by introducing a norm in the complex stress-strain space where we write Z e = (E e , S e ), with E e ∈ C d and S e ∈ C d for short and C > 0 is a parameter introduced to even out dimensions. Consider now two local harmonic histories z e (t) and z e (t) with Fourier representation where Z e , Z e ∈ C d ×C d are constant complex amplitudes and δ ω and δ ω are Dirac measures centered at the harmonic frequencies ω and ω , respectively. Then, we identity the distance between z (t) and z (t) with the flat-norm distance ẑ e −ẑ e betweenẑ e andẑ e , which is given explicitly by Prop. A.1 of Appendix A with λ = 1/ω cutoff and ω cutoff a cutoff frequency. As may be seen from Prop. A.1, the flat-norm distance between two Dirac measures quantifies the difference between their amplitudes as well as the distance between their points of application, which indeed defines a natural distance between Dirac measures. In particular the convergenceẑ e →ẑ e requires their amplitudes and their points of application to converge, i. e., Z e → Z e in C d × C d and ω → ω in R, as expected.
We can further extend the flat norm distance to two global harmonic histories, (21)ẑ = 2πZ δ ω ,ẑ = 2πZ δ ω , with Z , Z ∈ C N × C N , in a natural way as The corresponding DD problem is, then, i. e., we seek a material history y(t) and an admissible history z(t) that are closest to each other in the sense of the flat norm. We may regard (23) as an adversarial game with players y(t) and z(t) with goals F (y(t), z(t)) = y(t) − z(t) + I D (y(t)), (24a) where I D and I E are the indicator functions of D and E, respectively, i. e., they are zero in D and E, respectively, and +∞ elsewhere. In the game, y(t) plays to minimize F (·, z(t)), for given z(t), whereas y(t) plays to minimize G(y(t), ·), for given y(t). Thus, the goal of y(t) and z(t) is to remain in D and E, respectively, while simultaneously minimizing their distance, as required.
Apart from being natural for infinite wave trains such as harmonic histories, the flat norm distance affords the great practical advantage of enabling the use of data at frequencies other that the applied frequencies of excitation, with an appropriate penalty accorded to the discrepancy between the two. This property is specially advantageous when the data is sampled at discrete frequencies, which are likely to differ from the applied frequency.
e ) randomly from P. end for ii) Displacement problem. Solve for U (k) such that iv) Data assignment. for all e = 1, . . . , m do iv.a) Find (ω i ,

Numerical examples
In this section, we present two examples of application that illustrate the range and scope of the Data-Driven (DD) framework developed in the foregoing. We place particular focus on the convergence of the scheme with respect to data. Specifically, we consider a sequence P h of data sets providing an increasingly better representation of the exact complex modulus E(ω) and verify that, indeed, the DD solutions (y h (t)), z h (t)) converge to the exact solution y(t) = z(t).

4.1.
Fixed-point iteration. We solve the DD problem (23), or the equivalent game (24), by means of a fixed-point iteration adapted from [8]. We assume that the material is characterized by means of discrete data of the form P = {(ω i , E i ), i = 1, . . . , N }. The solver is summarised in Algorithm 1. The material is modelled as a linear viscoelastic material and is assumed to be exactly characterized by the Prony series measures experimentally by Knauss et al. [17], Fig. 1b. The reference solution for those data is shown in Fig. 2. Suppose that the exact characterization of the material, assumed to be as shown in Fig. 1b, is not known, but, instead, the material is characterized by a sequence of point data sets P h that sample the exact material behavior with increasing fidelity. In practice, the data could be the result of increasingly extensive successive experimental campaigns, or could be the result of diverse experimental campaigns of different degrees of reliability. The noise  in the data could be due to experimental measurement error, variability of the material samples, or combinations thereof. The specific data sets used in calculations are shown in Fig. 3. The data sets are synthetically generated by adding random noise to the true behavior of Fig. 1b. A steady increase in fidelity as achieved by: i) increasing the size of the sample, and ii) simultaneously decreasing the standard deviation of the noise. This type of data convergence has been termed uniform in [8,9]. The data-set sizes used in calculations are 10 3 , 10 4 , 10 5 and 10 6 . The truss problem is solved for each of the point data sets P h using the fixed-point iteration solver of Algorithm 1. The fast convergence of the iterative solver, which is attained in ten or fewer iterations, is remarkable given the combinatorial complexity of the data assignment operation and renders the calculations extremely fast. The normalized flat-norm error in the DD solutions, relative to the exact reference solution, is shown in Fig. 4 as a function of the point-data size. As is evident from the figure, the DD solutions exhibit a clear trend towards convergence at a rate of ∼ 1/3, which illustrates the convergence property of the general scheme.

Insonated gel block.
We demonstrate the ease of implementation of the DD approach in conventional finite-element codes by performing DD finite-element calculations using the commercial code Abaqus (Abaqus/Standard, Dassault Systemes Simulia, France). The problem concerns the response of a sample of agarose gel in the shape of a block under the action of sonic excitation on its boundary. Agarose gel is often used as a phantom material for soft biological tissue [16]. The analysis is structured as that of the preceding section. Fig. 5a shows the finite element model of a block 8 cm in size. The model is comprised of 4,913 nodes and 4,096 8-node hexahedral elements. A 4 cm 2 square on the top surface is subject to sonic pressure of 1 kPa amplitude and frequency Ω = 100Hz, while the bottom surface of the block is fully constrained. The material is assumed to be linear viscoelastic and to  be characterized by experimental data obtained using dynamic shear testing (DST) and magnetic resonance elastography (MRE) adapted from [18], Fig. 5b. The solution for the assumed exact material characterization is shown in Fig. 6. As in the preceding section, we assume next that the exact characterization of the material is not known exactly, but only through a sequence of point-data sets P h of increasing fidelity, Fig. 7. The data is generated synthetically from the exact characterization by the addition of random noise, representing experimental scatter, measuring error, sample-to-sample variability, multiple provenances, or other sources of error. The data sets contain 10 3 , 10 4 , 10 5 and, 10 6 points, respectively, and, simultaneously, exhibit decreasing amounts of scatter. In order to assess convergence, the normalized flat-norm error of the DD solutions relative to the limiting reference solution is plotted against the point-data set size, Fig. 8. The robust trend to convergence of the DD solutions, with convergence rate ∼ 1/2, is evident from the figure. Again remarkably, the fixed-point iteration solver converges in ten iterations or less, which renders the calculations exceedingly fast.

Concluding remarks
The present work showcases one particular way in which history data can be measured, codified and then used in model-free Data-Driven (DD) calculations, namely, through the use of the Fourier transform. This representational paradigm for histories is particularly well-suited in connection with long-term solid and structural dynamics where the histories of interest often take the form of 'infinite wave trains', a harmonic function being the simplest case in point [13]. Other representational paradigms for histories better suited for short-term behavior of materials with memory, such as viscoplastic materials, have been treated elsewhere [10]. It is therefore convenient that a number of experimental techniques supply data pertaining to the Fourier representation, including Dynamic Mechanical Analysis (DMA), nano-indentation, Dynamic Shear Testing (DST) and Magnetic Resonance Elastography (MRE), among others.
In the context of DD methods, in particular, it becomes necessary to measure distances between trial histories, satisfying compatibility and equilibrium, and material histories codified in material-data repositories. The Fourier transforms of such 'infinite histories' are localized measures in Fourier space, a Dirac-delta being the simplest case in point. There are a number of natural ways to measure distance between such measures, most of them variants of the Wasserstein distance from optimal transport [14], which quantify both the difference in amplitude and in the 'center of mass' of the measures. We have chosen to measure distance in Fourier space by means of the flat norm and shown that it lends itself to a straightforward implementation while providing the required control over the approximating solutions. Another advantage of the flat norm in Fourier space is that it allows data at different frequencies to be 'mixed and matched'. This feature is particularly useful when the data is sparse and specific frequencies are sampled only sparingly.
Finally, we emphasize that all DD calculations presented in this work have been directly predicated on material data without any recourse to material modeling of any type at any stage of the calculations. The ability to make quantitative-and convergent-predictions directly from material data is remarkable and sets forth a new model-free paradigm in computational mechanics.
where B ∈ S and B ⊥ in the orthogonal complement S ⊥ .
We state the result as a proposition for ease of reference.
Proposition A.1. Let B ∈ C n , α, β ∈ R. Let S be a proper subspace of C n and S = {Aδ α : A ∈ S}. Suppose that C n is metrized by a norm · derived from a Hermitian inner product. Then, Solving fo µ gives (33), whence (32) follows. We verify that which bears out the assumption (35). The condition Otherwise, the minimizer is A = 0 and the corresponding distance which completes (32).