Mixed methods for viscoelastodynamics and topology optimization

ABSTRACT. A truly-mixed approach for the analysis of viscoelastic structures and continua is presented. An additive decomposition of the stress state into a viscoelastic part and a purely elastic one is introduced along with an Hellinger-Reissner variational principle wherein the stress represents the main variable of the formulation whereas the kinematic descriptor (that in the case at hand is the velocity field) acts as Lagrange multiplier. The resulting problem is a Differential Algebraic Equation (DAE) because of the need to introduce static Lagrange multipliers to comply with the Cauchy boundary condition on the stress. The associated eigenvalue problem is known in the literature as constrained eigenvalue problem and poses several difficulties for its solution that are addressed in the paper. The second part of the paper proposes a topology optimization approach for the rationale design of viscoelastic structures and continua. Details concerning density interpolation, compliance problems and eigenvalue-based objectives are given. Worked numerical examples are presented concerning both the dynamic analysis of viscoelastic structures and their topology optimization.


INTRODUCTION
iscoelasticity is a constitutive feature of materials that finds applications in several areas such as structural engineering (concrete), biomedics (soft tissues) and also synthetic materials and polymers fabrication [1].From a modeling viewpoint, the most classic approach to viscoelasticity is based on hereditary integral formulations that rely on the interpretation of viscoelastic materials as materials with memory: the current stress is based on a time integral of the strain history.Given such a complex constitutive response, the design of viscoelastic structures is quite a challenging task that calls for a robust numerical model and a rationale design approach.Within this scenario, object of this paper is the proposal of a truly-mixed finite element method for the analysis of viscoelastic structures coupled to a topology optimization approach with the ultimate goal of the development of an integrated analysis-design tool for viscoelastic systems.For this paper's sake, eigenvalue-based objectives shall be pursued within a few numerical applications that are based on the purely elastic cases discussed in [2].As to the mixed finite-element formulation, the truly-mixed setting in a general Hellinger-Reissner framework has been selected mainly because of the need to pursue an accurate approximation of the stress field, feature that is not shared by any standard displacement-based approach.Furthermore, such a formulation is well known to pass the inf-sup condition [3] even in the limiting case of incompressible materials with no further tricks or enrichments and this could represent a V key issue when modeling and designing rubber devices such as base isolators.The adopted formulation resembles the one presented in [4] in that an additive decomposition of the stress field is considered that finds its motivation if the adoption of a Maxwell-type phenomenological model.However, a weak enforcement of the stress tensor symmetry is considered in [4] as opposed to a strong approach adopted herein.The Arnold-Winther finite element is used as to the discretization of the stress field [5].To the author's knowledge, it has never been used for the analysis of viscoelastic structures whereas a few applications in a purely elastic framework are available in the literature [6].Optimal design with respect to eigenvalues is a topic that has received much attention in the last decades and [7,8] and [9] may be cited among the most interesting contributions in this respect.However, optimal design of nonclassically damped viscoelastic structures is far from being a mature topic and further research seems to be in order.To gain insight into such a complex problem, a viscoelastic thin-beam truly-mixed formulation is introduced that has the merit to allow a deep understanding of the spectral properties of the discretized structure so as to end up with convincing optimal topologies.The first three eigenvalues shall be object of optimization and a strict relation between the optimal density distribution and the relevant eigenmode shapes clearly determined.Though preliminary, not many such results are available in the literature and should open the way to more complex results concerning two-dimensional systems.

FORMULATION AND DISCRETIZATION (CONTINUUM CASE)
oal of this section is the definition of a truly-mixed variational formulation for two-dimensional viscoelastic continua.Based on a parallel phenomenological model, compatibility, equilibrium and viscoelastic laws are written in strong form so as to allow the definition of a truly mixed variational formulation of Hellinger-Reissner type.Eventually, the Arnold-Winther finite-element is introduced along with some technical details concerning its implementation.

Strong form
Reference is made to Fig. 1 for the standard viscoelastic solid model that is the basis for the continuum and thin-beam models to be developed hereafter.As usual when adopting Hellinger-Reissner variational principles, compliance tensors relating strains to stresses are introduced that allow one to write where 0 E A and 0 V A are the elastic and viscous compliance tensors of the viscoelastic component, 1

E
A is the elastic compliance tensor that is in parallel to the viscoelastic one and v is the velocity field.One should notice that a stressvelocity formulation is being used that presents several advantages over more classical stress-displacement approaches, including the ease with which dynamic effects may be considered in the analysis.Therefore, compatibility relations are written in terms of strain velocities as whereas the dynamic equilibrium reads Truly-mixed variational formulation By observing that the total stress  may be additively decomposed as 0 1      , the continuous variational formulation of the problem at hand may be obtained by eliminating the strain tensor  in Eq. ( 1) and (2) testing the resulting equation by two virtual stress fields 0 1 ,   and the equilibrium Eq. ( 3) by a virtual velocity field w so as to write: In more compact form Eq. ( 4) may be rewritten as It should be emphasized that static and dynamic models differ only in that the mass matrix M may be neglected or considered in statics or dynamics, respectively.Eq. ( 5) amounts to a system of ordinary algebraic-differential equations (ADE) in the former case, of differential equations in the latter.

Imposing Dirichlet and Neumann boundary conditions
Within the adopted truly mixed formulations, Dirichlet boundary conditions on the velocity field are imposed weakly thanks to the line integral that appears on the right hand side of Eq. ( 4)1 and (4)2, and in fact, taking no action on a boundary amounts to imposing a homogeneous condition 0, on u v   .Conversely, Neumann traction boundary conditions n t    are to be imposed strongly.However, when an additive stress decomposition is adopted, as is the case herein, i.e.
 nor 1  are known but their sum.Therefore a Lagrange multiplier approach is adopted to enforce weakly a condition of type 0 1 ( ) n t      .Therefore, the resulting linear algebraic-differential system takes the form that may be compactly rewritten as with an obvious definition of vectors and matrices.

The Arnold-Winther finite element
The triangular Arnold-Winther finite element used in this paper is the lowest-order of the family of finite elements introduced in the pioneering paper [5].Fig. 2 shows the relevant degrees of freedom that may be listed as follows.
Arnold-Winter Stress DOFS Dispacement DOFS As to the stresses, one should notice that the symmetry of the stress tensor is imposed strongly so that the components to be approximated are 11 22 12 , ,    and one ends up with 24 degrees-of-freedom: , , As to the displacements, the two components x u and y u are linear on each element and globally discontinuous.

Implementation details
No doubt that implementing the Arnold-Winther finite element represents a severe challenge, especially if the goal of minimizing the memory storage is pursued as it should if one recalls that this element is far more expensive than more conventional ( , ) u p elements and other ( ) H div elements such as the Johnson-Mercier element [10].A possible implementation of the Arnold-Winther element is proposed in [6] that exploits indirect evaluations of the relevant stiffness matrices.Within the present paper, a different semi-analytical approach has been followed inspired by classical isoparametric elements according to which stress shape functions are first computed analytically on a parent triangular domain.Though not being isoparametric, the Arnold-Winther element enjoys the well-known Piola transformation property where B is the Jacobian of the affine transformation between the reference and actual configuration.This allows to compute the stiffness matrix in any actual deformed configuration.

FORMULATION AND DISCRETIZATION (THIN BEAM CASE)
he present section mimics the content of previous one but for the fact that it focuses on a Bernoulli viscoelastic beam.

Truly mixed variational form and FEM discretization
The phenomenological model of Fig. 1 is now used to express the adopted viscoelastic behavior in terms of bending moments M and dual curvatures  .Eq. ( 1) is then rewritten in the form and coupled to the equilibrium equation to arrive at the following Hellinger-Reissner truly-mixed formulation: find , which is apparently the thin beam counterpart to Eq. ( 4).As to the finite element discretization, the bending moment is approximated by means of cubic Hermite polynomials so that each node has two degrees of freedom, i.e. the bending moment itself and its first spatial derivative, i.e. the shear value.By analogy with the continuum case, the velocities are interpolated by elementwise linear, globally discontinuous polynomials.For later use, Eq. ( 14) is rewritten in compact form as where and The class of eigenvalue based optimal design problems considered herein may be written as where ( ) F s is a scalar-valued objective function depending on the eigenvalues s of the problem, (20)2 is the eigenproblem that enters the optimization procedure as a constraint and (20) 3 and (20) 4 are global and local material density bounds, respectively.A classical SIMP model relating elastic and viscoelastic moduli to the third power of the material density is considered whereas a linear dependence between the mass density and the material density is adopted, i.e. the two densities coincide as a matter of fact.The solution of the topology optimization problem is sought by means of a sequential quadratic programming approach that has been implemented in Matlab in a self-consistent code that include the finite-element analysis part as well.

Numerical solution of the Algebraic-Differential problem
The solution of the (forced) algebraic-differential equation that governs the dynamics of the viscoelastic problem is far from being trivial, especially in the 2D case wherein the algebraic part of the system has considerable dimension because of the need to impose several Cauchy boundary conditions on all the unconstrained sides of the domain.For this paper's sake, a second-order accurate numerical approximation consisting of a trapezoidal rule followed by a 2-step backward difference scheme, as suggested in [4], has been adopted:

NUMERICAL STUDIES
Analysis and topology optimization of a clamped/simply-supported beam clamped/simply-supported beam of length 8 L  is investigated first.The initial design consists of a uniform cross-section and density.If x denotes the local density of the beam (that may be interpolated element wise or the node level), the following SIMP-like interpolations are adopted: where 3 p  and the minimum and maximum values adopted in the numerical simulations are respectively given in Tab.Before introducing a few choices for ( ) F s to be tested numerically, a few comments on the spectral properties of the problem at hand are in order.Notice in fact that neither is the system classically damped, nor enjoys the classical structure of state-space structural dynamics whose eigenproperties are well known.Reference is made to the doubly-clamped beam but similar considerations apply to all constraint cases as well.If NUMEL denotes the number of elements, the total dimension N of the problem may be written and decomposed as

Velocity dofs Moment dofs
Viscoleastic law dofs The 2 NUMEL  constitutive laws (two devices in parallel per element) bring into the formulation an equal number of purely real eigenvalues, whereas moments and velocities dofs, globally amounting to 4 2 NUMEL   , are associated to 2 NUMEL  pairs of complex and conjugate eigenvalues and 2 null eigenvalues (that become 3 in the case of a clampedsupported beam).As to the objective function ( ) F s in Eq. (20), the following three cases are considered: The first objective is very classical and aims to the maximization (of the imaginary part of) the lowest eigenvalue.The second objective aims to maximize the spectral band between the first two adjacent eigenvalues to avoid dangerous mode superpositions.The third function follows the same streamline as the second but involves the third eigenvalue as well.For the three optimal design problems with objectives 1 2 3 , , F F F , Fig. 3, 4 and 5 respectively present the convergence path (on the left) and the optimal density distribution along the beam axis (on the right).Fig. 6 and 7 show the first three velocity and moment eigenmodes, respectively, and demonstrate the strict correlation between the optimal density at convergence and the spectral properties of the beam.Tab. 3 shows a few values for the three design cases.The three objective functions 1 2 3 , , F F F are respectively dominated by the first, second and third eigenmode and therefore each design maximizes its own objective (boldface values on the diagonal in Tab. 3) but, with no surprise, happen to lead to small values when tested on a different objective than the one used in the optimization (first three columns in Tab. 3).To assess accuracy and convergence of the proposed approach reference is made to Fig. 8, 9 and 10 that, for the objectives 1 2 3 , , F F F respectively, show the optimal values taken on by the objective functions versus the number of elements used for the discretization (on the left) and relevant optimal designs (on the right).The method is clearly robust in that optimal solutions are clearly mesh insensitive and even using a coarse mesh of 8 elements a good approximation to the exact optimal solution (that is of course analytically unknown and is herein assumed to coincide with the numerical one computed with a mesh of 64 elements) is found.Of course, the higher the mode number that dominates the objective function the finer the mesh needed to find an accurate solution: one may in fact see that the solution using 8 elements (Fig. 8 on the right) nearly coincides with the exact one when the objective function is 1 F that depends on the first eigenmode only, whereas the accuracy decreases when solving the optimal design problems governed by 2 F and 3 F that involve higher eigenmodes.However, this should not be seen as a limitation of the proposed optimization approach in that when a coarse mesh is not capable to approximate accurately the eigenmodes for a given material density, the solution of the optimization procedure cannot be any better.

Analysis of a plane cantilever
The plane cantilever of Fig. 11, uniformly loaded on the right side is considered.The spatially uniform load is modulated in time by means of the function As to the physical properties of the structure the following values are adopted as to the isotropic elastic and viscous compliance tensors that are defined in terms of Lamé constants:      Finally, Fig. 17 shows the time-variation of the stress xx  at the upper-left corner of the cantilever.

CONCLUSIONS
truly-mixed variational formulation for the analysis of viscoelastic continua and thin beams has been presented that uses stresses (moments) as primary variables and velocities (instead of displacements) as Lagrange multipliers.A topology optimization method for viscoelastic devices has then been developed and applied to thin beams within a few representative numerical simulations.Ongoing extensions include applications to two dimensional viscoelastic devices with respect to eigenvalue-based design objectives.

Figure 2 :
Figure 2: Degrees-of-freedom of the Arnold-Winther stress element.
the three components of the stress tensor 11 22 12 , ,    at each vertex of the triangle (9 dofs),  the moments of order zero and one of the traction vector n   along each edge of the triangle (12 dofs),  the averages of the components of the stress tensor over the triangle, i.e. 11 22 12

T
EIGENVALUE-BASED TOPOLOGY OPTIMIZATION AND SOLUTION OF THE DIFFERENTIAL-ALGEBRAICPROBLEMEigenvalue based topology optimizationhe abstract unforced equation of motion for the two systems intriduced above reads where

Figure 3 :
Figure 3: Objective function 1: path to convergence and optimal density.

Figure 4 :
Figure 4: Objective function 2: path to convergence and optimal density.

Figure 5 :
Figure 5: Objective function 3: path to convergence and optimal density.

Figure 8 : 1 sF
Figure 8: 1 s F : Convergence vs number of elements (left) -Optimal densities vs number of elements (right).

Figure 9 : 2 sF
Figure 9: 2 s F : Convergence vs number of elements (left) -Optimal densities vs number of elements (right).

Figure 10 : 3 sF
Figure 10: 3 s F : Convergence vs number of elements (left) -Optimal densities vs number of elements (right).
further considered.Fig.12, 13 and 14 show the maps of the stress components , ,

Figure 15 :
Figure 15: Evolution of xx  (normal stress) with time.

AFigure 17 :
Figure 17: Evolution of xx  (normal stress) with time at the upper left corner of the cantilever.

Table 1 :
SIMP interpolation: minimum and maximum values.