Next Article in Journal
A Note on Generalization of Combinatorial Identities Due to Gould and Touchard
Next Article in Special Issue
Solutions for Some Mathematical Physics Problems Issued from Modeling Real Phenomena: Part 1
Previous Article in Journal
Metaheuristic-Based Hyperparameter Tuning for Recurrent Deep Learning: Application to the Prediction of Solar Energy Generation
Previous Article in Special Issue
Ekeland Variational Principle and Some of Its Equivalents on a Weighted Graph, Completeness and the OSC Property
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Systematic Approach to Standard Dissipative Continua

Fraunhofer Institute for Ceramic Technologies and Systems IKTS, Winterbergstr. 28, 01277 Dresden, Germany
Axioms 2023, 12(3), 267; https://doi.org/10.3390/axioms12030267
Submission received: 31 January 2023 / Revised: 25 February 2023 / Accepted: 27 February 2023 / Published: 4 March 2023
(This article belongs to the Special Issue Principles of Variational Methods in Mathematical Physics)

Abstract

:
Many isothermal dissipative continuum problems can be formulated in a variational setting using the concept of “standard dissipative continua”. A major advantage of this approach is that complex problems can be cast into a compact, thermodynamically consistent formulation based on a single space–time continuous functional together with a corresponding variational principle. Formulating the problem in terms of a functional provides an immediate avenue for performing spatial and temporal discretization, which are the prerequisites for a numerical solution. Within the present contribution, a novel systematic approach to standard dissipative formulations is proposed, with the main goal being the development and implementation of generic procedures and algorithms for the formulation as well as the computational solution of a subset of isothermal dissipative continuum problems. In order to demonstrate the capabilities of the approach, its application to example problems is discussed.
MSC:
35A15; 49S05; 49-04; 65M60; 74F99; 68V99

1. Introduction

Starting from classical continuum theories, recent decades have seen substantial developments towards the modeling of nonlinear dissipative multi-physics phenomena. Examples include electromechanical coupling in ferroelectric materials [1], magnetoelectromechanical coupling in multiferroic materials [2] as well as electrochemomechanical coupling in batteries [3] and hydrogels [4]. Although most continuum concepts are, in principle, readily extended to account for such multiphysics effects, the formulation of thermodynamically consistent theories and their solution by numerical methods becomes an increasingly tedious task as the number of coupled physics and, therefore, unknown variables increases. This is particularly of concern for:
1.
The model formulation in terms of a number of partial differential (or differential algebraic) equations, which may often be combined in many different ways to form a final set of space–time continuous equations;
2.
The mathematical analysis of the space–time continuous equations, e.g., with regard to the well posedness of the problem and stability properties;
3.
The development of consistent, robust and efficient spatial and temporal discretization schemes;
4.
The numerical solution of the space–time discrete problem.
An established approach to reduce the mentioned difficulties is to formulate continuum models in terms of variational principles based on a single functional of the unknowns instead of partial differential equations, so that the solution of the problem furnishes an extremum or a saddle point of the functional. This approach often provides a compact statement of the problem, which clearly exposes its structure, and thereby facilitates the mathematical analysis as well as the development of robust and efficient discretization and numerical solution procedures. The latter has, e.g., been demonstrated by Simo and Honein [5], Ortiz and Stainier [6], and Hackl and Fischer [7] in the context of plasticity and viscoplasticity, by Yang et al. [8] in the context of thermomechanically coupled problems, and by Carstensen et al. [9] and Ortiz and Repetto [10] in the context of micro-structure formation. Of particular interest in this regard is the concept of “standard dissipative continua” (see, e.g., Halphen and Nguyen [11], Germain et al. [12]), which effectively involves minimum or saddle point formulations, with thermodynamic consistency being ensured by choosing a so-called dissipation functional to satisfy certain conditions. Examples illustrating the wide applicability of this concept are, e.g., the early work of Biot [13] on fluid transport in porous media as well as a multitude of contributions covering the homogenization of inelastic micro-structures [14,15], gradient-enhanced plasticity [16,17,18,19,20], rate-dependent ferroelectricity [21], magnetostriction with hysteresis [22], electro-magneto-mechanics [2], fluid transport in porous media [23], Cahn–Hilliard-type and classical diffusion [24,25], and fracture in lithium-ion batteries [26].
Regardless of whether differential (algebraic) equations or a functional together with a variational principle are chosen as the point of departure, the mere implementation effort required to obtain validated numerical software remains a substantial obstacle, especially in the case of multiphysics problems. Due to this aspect, many sophisticated continuum formulations are often only demonstrated for highly simplified situations (e.g., one-dimensional problems), while their application to real-world examples remains elusive. This aspect provides the main motivation for the present contribution. In particular, a systematic approach to standard dissipative formulations is proposed, which covers the steps from problem formulation in the space–time continuous setting through its discretization to implementation.
The plan of the paper is as follows: Firstly, the standard dissipative space–time continuous problem formulation is formalized without reference to particular physics. Secondly, simple approaches for the spatial and temporal discretization of the space–time continuous problem are described. Thirdly, a preliminary implementation based on the open source finite element library deal.II (Alzetta et al. [27], Bangerth et al. [28]) is briefly discussed. Finally, the usefulness of the approach is shown by the application to two multiphysics examples. The first example addresses the buckling behavior of a viscoelastic dielectric elastomer subjected to an electric voltage, while the second example concerns the swelling/shrinking of a hydrogel upon a change of exterior ion concentrations. Both examples include numerical convergence studies with regard to spatial and temporal discretization in order to show that the expected numerical properties indeed translate to practical computations.

2. Formalization of the Space–Time Continuous Problem

In this section, a generic formulation for standard dissipative continua is proposed, with the scope being limited to isothermal problems with elliptic or parabolic structure. It is emphasized in this context that more general standard dissipative formulations than the one presented below are possible. However, certain simplifications are unavoidable in view of implementation aspects. The proposed formulation essentially represents a trade-off between the generality of the formulation and complexity of the implementation into numerical software. Furthermore, it is remarked that most of the theoretical developments are, at least as far as the general ideas are concerned, in some form already available in the pertinent literature. Here, the intended contribution is to compile the theory in a form suitable for subsequent implementation.
The contents of this section are summarized as follows: Firstly, the computational domain including subdomains, boundaries and interfaces is described together with the notation required for bookkeeping purposes. Secondly, the notions of “generalized fields” and “dependent fields” are introduced, and the form of the functionals considered in this work is specified. Despite the fact that these aspects are of subordinate importance from a physical standpoint, they are indispensable with regard to the desired problem-independent implementation of the approach into computer code. Thirdly, the unknowns involved in the statement of the problem are grouped into several categories based on their role in the variational formulation. Fourthly, the space–time continuous incremental potential is introduced, followed by a statement of the space–time continuous incremental variational principle assumed in the present work. Finally, ways to systematically reduce the number of unknowns involved in the problem formulation are described.

2.1. Computational Domain

The Euclidean d-space R d with d = 2 or d = 3 is considered for a range of time T = [ t 0 , t e ) , with t 0 and t e being the initial and the final times, respectively. In what follows, X R d denotes a spatial location and t T an instant of time. As illustrated in Figure 1, R d is assumed to be split into (i) the “environment” Ω 0 ; (ii) the bounded domain Ω , which is further subdivided into N Ω “domain portions” Ω α labeled by α ; and (iii) the bounded, sufficiently smooth d 1 dimensional interface Σ , which is further subdivided into N Σ “interface portions” Σ γ labeled by γ . It is generally assumed that the boundary of each domain portion Ω α is a part of Σ . This means that Ω α Σ = Ω α . Moreover, it is understood that Ω and Σ have no points in common, i.e., Σ Ω = . In the sequel, the union of Σ and Ω will be referred to as the “computational domain”. This computational domain is generally assumed to be kept fixed. This means that it does not change with time. Together with the interface Σ , a piece-wise continuous unit normal field N ( X ) is introduced. Using this unit normal field, the ‘+’ and ‘−’ sides of the interface at a location X Σ are defined such that N ( X ) points from the − to the + side.
Remark 1.
  • An interface portion Σ γ need not necessarily be aligned with an interface between different domain portions Ω α or an interface between a domain portion Ω α and the environment Ω 0 . This means that a “slit” within a domain portion is permissible, as can be seen in the interface portion Σ 1 in Figure 1.
  • A more complete description would also involve lower-dimensional entities (lines, line portions, and points in three dimensions; points in two dimensions). However, these are firstly rarely relevant to practical examples, and secondly, these are more difficult to implement. Therefore, they are ignored below.
  • The introduction of different domain and interface portions concerns in particular multiphysics phenomena, with different physics being operative in different regions, and, possibly, combined with interfacial effects.
  • Keeping the domain of computation fixed does not imply that a Lagrangian or Eulerian viewpoint is adopted. In fact, even arbitrary Lagrangian Eulerian formulations are possible.

2.2. Generalized Fields, Dependent Fields and the Form of Functionals

Before discussing the standard dissipative formulation considered in this work,
(i)
the notions of “generalized fields” and “dependent fields” are introduced, and
(ii)
the form a functional can take within this work is narrowed down.

2.2.1. Generalized Fields

A generalized field G = G ( X , t ) comprises:
1.
A number of N G , Ω time-dependent scalar fields, which are formally defined on the entire computational domain, but which may only be non-zero on Ω . These fields are collected in the vector G Ω = G Ω ( X , t ) and will be referred to as “domain-related fields”. The latter may represent, e.g., mechanical displacements, velocities, pressures, electric fields, etc.
2.
A number of N G , Σ time-dependent scalar fields, which are formally defined on the entire computational domain, but which may only be non-zero on Σ . These fields are collected in the vector G Σ = G Σ ( X , t ) and will be referred to as “interface-related fields”. The latter may represent, e.g., curvatures, electrical surface charges, interfacial reaction rates, etc.
3.
A number of N G , C time-dependent scalars. These scalars are collected in the vector G C = G C ( t ) ; and they may represent quantities such as integral electric currents, displacements and rotations of rigid domain portions, etc.
With these definitions,
G = G Ω G Σ G C
is as a vector of dimension N G = N G , Ω + N G , Σ + N G , C .
The derivative of G with respect to X , i.e., the gradient, is introduced according to
G = G Ω Σ G Σ 0 ,
with Σ indicating the surface gradient, and G Ω = 0 if X Σ and Σ G Σ = 0 if X Ω . It is noted that G may be interpreted as an N G × d matrix.
The derivative of G with respect to t, i.e., the rate, is denoted by
G ˙ = G ˙ Ω G ˙ Σ G ˙ C .
Two generalized fields A and B may be concatenated to form a single generalized field C . In this case, it is understood that the components are rearranged such that the domain-related fields, the interface-related fields and the scalars are contiguous again. This means that:
C = A B = A Ω B Ω A Σ B Σ A C B C = C Ω C Σ C C .
This procedure can be applied recursively to concatenate more than two generalized fields.
Remark 2.
Vector/tensor fields are considered to be represented by an appropriate number of scalar fields, with the individual scalar fields typically representing vector/tensor coordinates with regard to a particular basis system.

2.2.2. Dependent Fields

For each generalized field G , a number of N g , Ω “domain-related dependent fields” and a number of N g , Σ “interface-related dependent fields” are introduced. These dependent fields are collected in the vectors g Ω and g Σ , respectively. In this section, the convention will be followed that upper case letters denote generalized fields, while the corresponding lower case letters denote the associated dependent fields.
The purpose of the dependent fields is the introduction of variables, in terms of which the functionals governing the standard dissipative formulation can be conveniently expressed without direct use of spatial derivatives. A typical example for a dependent field is the mechanical strain field, which is used to introduce the strain energy density in continuum mechanics.
The dependent fields g Ω and g Σ are assumed to be linearly related to the generalized field G . In particular,
( g Ω ) λ = ϵ = 1 N G , Ω a λ ϵ G , Ω ( G Ω ) ϵ + b λ ϵ G , Ω · ( G Ω ) ϵ + ι = 1 N G , C c λ ι G , Ω ( G C ) ι + d λ G , Ω
if X Ω and ( g Ω ) λ = 0 otherwise, with a λ ϵ G , Ω , c λ ι G , Ω and d λ G , Ω being constants, b λ ϵ G , Ω constant vectors of dimension d, and ( Y ) λ denoting the component λ of a vector Y . Similarly,
( g Σ ) ν = η = 1 N G , Σ a ν η G , Σ ( G Σ ) η + b ν η G , Σ · ( Σ G Σ ) η + ϵ = 1 N G , Ω a ν ϵ G + ( G Ω | + ) ϵ + b ν ϵ G + · ( G Ω | + ) ϵ + ϵ = 1 N G , Ω a ν ϵ G ( G Ω | ) ϵ + b ν ϵ G · ( G Ω | ) ϵ + ι = 1 N G , C c ν ι G , Σ ( G C ) ι + d ν G , Σ
if X Σ and ( g Σ ) ν = 0 otherwise, with a ν η G , Σ , a ν ϵ G + , a ν ϵ G , c ν ι G , Σ and d ν G , Σ being constants, and b ν η G , Σ , b ν ϵ G + , b ν ϵ G constant vectors of dimension d. In the latter equation as well as in the following, the notation | + / indicates that reference is made to the limit value as points X Σ are approached from the + and the − side of the interface Σ , respectively.
Remark 3.
  • The restriction of g Ω and g Σ to a dependence on first spatial derivatives of the generalized field G may seem rather limiting. However, the formulation to be described does enable the effective incorporation of higher-order derivatives through the introduction of auxiliary quantities. In fact, including higher spatial derivatives would, in a finite element context, require higher-order continuous finite elements; and the latter are rarely used in practice due to the superiority of mixed formulations eliminating the higher spatial derivatives through auxiliary quantities.
  • The choice of the dependent fields is usually subject to a substantial degree of freedom and will typically be determined by what is most convenient for the definition of the governing functionals.
  • The constants and constant vectors appearing in (5) and (6) may be made dependent upon X without affecting the discussion below. Furthermore, nonlinear dependencies of the dependent fields on the generalized fields and derivatives thereof may be allowed for . However, such dependencies are omitted here as they have not yet been implemented (nonetheless, they can be incorporated through constraints, which can be implemented using Lagrangian multipliers).

2.2.3. Form of Functionals

Let H [ A α ; B β , t ] be a functional of a number of generalized fields A α indexed by α and a number of generalized fields B β indexed by β , with the semicolon indicating that the quantities to the right are regarded as parameters when the variations of the functional H are computed. In the sequel it is assumed that every such functional is of the specific form
H [ A α ; B β , t ] = Ω h Ω ( a α , Ω ; b β , Ω , X , t ) d V + Σ h Σ ( a α , Σ ; b β , Σ , X , t ) d S + h C ( A α , C ; B β , C , t ) ,
where d V is the volume element, d S the interface element, and h Ω , h Σ and h C are functions.
Remark 4.
The previous assumption regarding the form of a functional will eventually ensure that the incremental variational principle to be stated below implies a set of partial differential (algebraic) equations and equations without spatial dependency (which may, however, involve integrals).

2.3. Unknown Variables

Conceptually, three different types of unknowns are considered, leading to the definition of the three following generalized fields:
1.
A generalized “state variable” Q , which is used to describe the state within the computational domain. Q may, besides classical “thermodynamic state variables”, also comprise a number of “auxiliary state variables”, which will be used in the usual way as Lagrangian multipliers enforcing constraints on the admissible thermodynamic states (incompressibility constraints, electroneutrality constraints, etc.). In particular, Q is assumed to be of the format
Q = U P ,
where U comprises the thermodynamic state variables, while P comprises the auxiliary state variables/Lagrangian multipliers. Q is subject to an initial condition of the form
Q ( X , t 0 ) = Q 0 ( X ) ,
where Q 0 may need to satisfy certain consistency requirements as will be remarked later.
2.
A generalized “process variable” V , which may be interpreted as an auxiliary quantity used to constrain the ways the thermodynamic state can change in that only certain combinations of the rates Q ˙ and V ˙ are admissible. V is, without loss of generality, subjected to the initial condition
V ( X , t 0 ) = V 0 ( X ) = 0 .
3.
A generalized Lagrangian multiplier M , which will be used to actually incorporate constraints for Q ˙ and V ˙ .
Remark 5.
  • The explicit distinction between state variables and process variables appears to be a feature distinguishing the present standard dissipative formulation from what is available in the literature to date, although it is acknowledged that, e.g., the works of Biot [13] and Miehe et al. [23], Miehe et al. [24], Miehe et al. [25] implicitly contain such a distinction. Physically, the introduction of process variables may, e.g., be motivated by the following two simple examples: (i) diffusion processes, in which molar concentrations are changed by species fluxes; and (ii) chemical reactions, in which the amounts of the substances involved are changed by chemical reaction processes.
  • Ultimately, only the rate of the process variable, V ˙ , will enter the formulation. Hence, it would seem more natural to directly introduce the rate as an unknown. However, the viewpoint adopted here proves beneficial from the implementation point of view.

2.4. Space–Time Continuous Incremental Potential

The space–time continuous problem formulation in terms of an incremental variational principle (herein, the term “incremental” refers to a formulation in terms of rate-type quantities, and does not refer to a time-discrete formulation as it is often the case in the literature) is based on
1.
A Helmholtz free energy functional Ψ [ Q ] . Specifically, the Helmholtz free energy functional is assumed to be of the form
Ψ [ Q ] = Ψ U [ U ] + Ψ P [ P , U ] ,
with Ψ P being linear in its first argument. In general, Ψ will be assumed to be Fréchet differentiable. It is remarked in this context that the constraints for the thermodynamic states U are obtained by requiring the first variation of Ψ P with respect to P to be zero; and Ψ U represents the physical Helmholtz free energy for all thermodynamic states consistent with these constraints.
2.
A power functional P [ U ˙ , V ˙ ; U , t ] accounting for the influence of the environment on the computational domain under consideration. P is assumed to be linear in the pair ( U ˙ , V ˙ ) .
3.
A dissipation functional Δ ˚ [ U ˙ , V ˙ ; U , t ] , which is related to the description of dissipative effects. It is assumed that Δ ˚ is (i) non-negative; (ii) convex in the pair ( U ˙ , V ˙ ) ; and (iii) satisfies the condition Δ ˚ [ 0 , 0 ; U , t ] = 0 . The latter conditions are, within the standard dissipative framework, sufficient for the satisfaction of the second law of thermodynamics as will be remarked later.
4.
A Lagrangian multiplier functional L [ M , U ˙ , V ˙ ; U , t ] incorporating constraints for U ˙ and V ˙ . L is assumed to be of the form
L = L 1 [ M , U ˙ , V ˙ ; U , t ] + L 2 [ M ; U , t ] ,
where L 1 is bilinear with regard to M and the pair ( U ˙ , V ˙ ) and L 2 is linear in M .
Finally, the space–time continuous incremental potential
Π [ Q ˙ , V ˙ , M ; Q , t ] = Ψ ˙ [ Q , Q ˙ ] + Γ [ U ˙ , V ˙ , M ; U , t ]
is introduced, where
Γ [ U ˙ , V ˙ , M ; U , t ] = P [ U ˙ , V ˙ ; U , t ] + Δ ˚ [ U ˙ , V ˙ ; U , t ] + L [ M , U ˙ , V ˙ ; U , t ] .

2.5. Space–Time Continuous Incremental Variational Principle

Using the previously introduced functionals, the space–time continuous incremental variational principle now requires that Q ˙ , V ˙ , M furnishes a stationary point of Π at every instant of time. This means that the problem can be formulated as
Find the unknowns ( Q , V , M ) such that Q ( X , t 0 ) = Q 0 , and V ( X , t 0 ) = 0 , and
0 = δ Π [ Q ˙ , V ˙ , M , δ Q ˙ , δ V ˙ , δ M ; Q , t ] = δ Ψ [ Q , δ Q ˙ ] + δ Γ [ U ˙ , V ˙ , M , δ U ˙ , δ V ˙ , δ M ; U , t ]
for all [ δ Q ˙ , δ V ˙ , δ M ] and all t ( T t 0 ) ,
where the notation is used here and henceforth that the first variation of a functional I [ x 1 , , x n ; y 1 , , y n ] is δ I [ x 1 , , x n , δ x 1 , , δ x n ; y 1 , , y n ] .
Remark 6.
  • Equation (15) may be seen as an extension of Biot’s equation, cf. Equation (2.14) in Biot [29] and Equation (2.14) in Biot [30].
  • For the general case, partial differential algebraic equations are involved and, therefore, (15) cannot generally be used as a “point-wise” equation to determine the rates Q ˙ , V ˙ based on the current thermodynamic state Q . Rather, neighboring instants of time need to be taken into consideration by differentiating the equations implied by (15) with respect to time. This aspect becomes particularly evident when considering the example of elasticity (which does not involve dissipative processes). In this case, (15) implies the mechanical equilibrium condition, which only involves the mechanical displacement and not the rate of the mechanical displacement. Therefore, differentiation with respect to time is necessary in order to obtain an equation in the rates.
  • In order for the problem to be well posed, further conditions need to be satisfied. Although a general discussion of this aspect without reference to a particular problem appears to be difficult, a few further remarks are made in the following:
    -
    The constraints represented by the functionals Ψ P and L must be consistent in that they are neither redundant nor contradictory.
    -
    Non-convex Ψ U and/or domains of admissible thermodynamic states may be a severe issue with regard to unique solvability and stability, as can be seen in the discussion in Abed-Meraim and Nguyen [31].
    -
    The initial conditions need to be consistent in that Q 0 must be such that a stationary point of Π exists at t = t 0 .
  • The problem formulation in terms of the first variation of Π implies the assumption of sufficient regularity of the functionals involved. This becomes a crucial aspect in the case of problems including rate-independent processes. In this case, the dissipation functional Δ ˚ contains contributions which are positively homogeneous of degree one in the rates. The associated non-differentiability of Δ ˚ invalidates (15), and hence it may appear that the formulation cannot be used for these cases. However, the situation can be remedied by employing the techniques of convex analysis and, in particular, replacing regular differentiation by the use of subdifferentials (as can be seen, e.g., in Rockafellar [32], pp. 213 ff.), thus extending the formulation to rate-independent processes.
  • The above formulation in terms of an incremental variational principle is motivated by the description of reversible isothermal processes, in which case the concepts of equilibrium thermodynamics apply and the free energy of the system under consideration is minimized at all times. In order to account for the actual dissipative nature of the problem under consideration, “generalized dissipative forces” are included into the formulation, which are mathematically characterized by the Fréchet derivative of Δ ˚ with respect to U ˙ and V ˙ . These generalized dissipative forces may be interpreted as those external loads, which would have to be applied by an external agency at a certain instant of time t to bring the system instantaneously into equilibrium. This approach is abstract in that these external loads are fictitious and cannot be applied in reality. However, accepting the approach based on the experience that it has been succesfully applied to a wide range of problems in the past, the situation is formally transformed from a problem of non-equilibrium thermodynamics into one of equilibrium thermodynamics, so that the well-known concepts of equilibrium thermodynamics can be used, and a mathematically sound formulation results, as can also be seen in Germain et al. [12]. Moreover, by ensuring that δ Δ ˚ ( U ˙ , V ˙ , U ˙ , V ˙ ; U , t ) 0 for any choice of ( U ˙ , V ˙ ) , it is guaranteed that the generalized dissipative forces are always associated with a negative power (i.e., a positive dissipation), such that the second law of thermodynamics is satisfied a priori and need not be discussed further. The assumptions made above ( Δ ˚ is non-negative, convex with respect to the pair ( U ˙ , V ˙ ) and satisfies the condition Δ ˚ [ 0 , 0 ; U , t ] = 0 ) indeed imply that δ Δ ˚ ( U ˙ , V ˙ , U ˙ , V ˙ ; U , t ) 0 .

2.6. Reduced Space–Time Continuous Incremental Variational Principles

In the following, two straightforward possibilities to reduce the set of unknown variables are described. It is emphasized that further possibilities exist, which are not discussed here. Furthermore, it is noted that the space–time discretization methods as well as the implementation discussed in the next sections are equally applicable to both, the original and the reduced formulations.

2.6.1. Reduction of the Set of Process Variables and Lagrangian Multipliers

The introduction of the generalized process variable V and the generalized Lagrangian multiplier M often leads to an excessive number of unknown fields. However, in many cases, it is possible to systematically eliminate V and M , or rather subsets thereof, before actually numerically solving the problem. In particular, the space–time continuous incremental functional Π only depends on V ˙ , but not on V itself. Furthermore, Π depends on V ˙ and M only through the contribution Γ . Thus, a reduced functional
Γ red [ U ˙ , V ˙ 0 , M 0 ; U , t ] = inf V ˙ 1 sup M 1 Γ [ U ˙ , V ˙ , M ; U , t ]
may be introduced by joint partial minimization and maximization, where
V = V 0 V 1 M = M 0 M 1 .
Then, δ Γ can be replaced by δ Γ red in (15) to obtain an equivalent formulation in the reduced set of variables ( Q , V 0 , M 0 ) .
Remark 7.
The selection of the variables V 0 and M 0 typically depends on the properties of the dissipation functional Δ ˚ . However, numerical considerations (e.g., with respect to the choice of finite elements) may also be important in this regard.

2.6.2. Reduction of the Set of State Variables

Another possibility to reduce the set of variables exists if certain state variables factor into Ψ , while Γ is independent of these variables. In order to discuss this case, it is assumed that
Q = U 1 U 0 P 0 = U 1 Q 0 ,
and that Γ neither depends on U 1 nor on U ˙ 1 , meaning that there exists a functional Γ red such that
Γ [ U ˙ , V ˙ , M ; U , t ] = Γ red [ U ˙ 0 , V ˙ , M ; U 0 , t ] .
Then, if Ψ is additionally convex in U 1 , a reduced thermodynamic potential
Ψ red [ Q 0 ] = inf U 1 Ψ [ Q ]
may be introduced; and δ Ψ and δ Γ can be replaced by δ Ψ red and δ Γ red in (15) to obtain an equivalent formulation in the reduced set of variables ( Q 0 , V , M ) .
Remark 8.
Using this approach, the linearity of Ψ with regard to P may be lost, which does, however, not have consequences for the applicability of (15).

3. Space–Time Discrete Formulation

The spatial discretization of the fields involved in Q , V and M is assumed to be performed using the finite element method in this work. In this regard, standard methods can be used, which will not be further discussed here. However, it is pointed out that the satisfaction of the inf sup condition is typically the most difficult part of spatial discretization if the variational problem to be solved is of saddle point type, see, e.g., Bathe [33] and Auricchio et al. [34].
With regard to temporal discretization, a number of simple and straightforward-to-implement one step schemes have been discussed and analyzed in a somewhat narrower context in a previous work [35]. These methods can be readily extended to the current formulation; and for the self-containedness of the present work, these will be briefly described below, although the reader is referred to [35] for details.
In the following, it is assumed throughout that the time interval [ t 0 , t e ] is divided by the discrete time points t n = t 0 + n Δ t , where n = 0 , 1 , N and Δ t = ( t e t 0 ) / N (however, the formulations are not restricted to uniform time increments). The numerically approximated values of Q , V and M at time t n are denoted by Q n , V n and M n , respectively; and the initial values at time t 0 are given by Q 0 ( X ) and V 0 ( X ) = 0 .

3.1. Variationally Consistent Method

The “variationally consistent method” is based on a straightforward discretization of the space–time continuous incremental functional Π . In particular, a time-discrete counterpart Δ Π n n + 1 of Π is introduced according to
Δ Π n n + 1 [ Q n + 1 , V n + 1 , M n + 1 ] = Ψ [ Q n + 1 ] Ψ [ Q n ] + Δ t Γ U n + 1 U n Δ t , V n + 1 V n Δ t , M n + 1 ; U n , t n + 1 .
Then, the stationary point of Δ Π n n + 1 is searched during each time step to obtain Q n + 1 , V n + 1 and M n + 1 based on the solution of the previous time step. The corresponding stationarity condition leads to the time-discrete variational formulation
Find the unknowns ( Q n + 1 , V n + 1 , M n + 1 ) such that
0 = δ Ψ [ Q n + 1 , δ Q n + 1 ] + δ Γ U n + 1 U n Δ t , V n + 1 V n Δ t , M n + 1 , U n + 1 U n Δ t δ U n + 1 , δ V n + 1 , Δ t δ M n + 1 ; U n , t n + 1
for all ( δ Q n + 1 , δ V n + 1 , δ M n + 1 ) and all n { 0 , 1 , N 1 } ,
which is a consistent approximation of (15). According to the analysis in [35], the method represents a first-order accurate time integrator, provided that the the problem is sufficiently regular.
Remark 9.
As a result of discretizing Π, the variational structure of the space–time continuous problem is preserved, which makes the method particularly attractive from the mathematical point of view. This concerns, e.g., the fact that the resulting finite element systems are symmetric.

3.2. The α -Family

The “ α -family” is based on the introduction of a real parameter α [ 0 , 1 ] , where α = 0 corresponds to a forward Euler type scheme, α = 1 / 2 to a Crank–Nicolson type scheme, and α = 1 to a backward Euler type scheme. Moreover, the following quantities are defined:
Q n α = ( 1 α ) Q n + α Q n + 1
M n α = ( 1 α ) M n + α M n + 1
t n α = ( 1 α ) t n + α t n + 1 .
Then, the time-discrete version of (15) reads
Find the unknowns ( Q n + 1 , V n + 1 , M n α ) such that
0 = + ( 1 α ) δ Ψ [ Q n , δ Q n + 1 ] + α δ Ψ [ Q n + 1 , δ Q n + 1 ] + δ Γ U n + 1 U n Δ t , V n + 1 V n Δ t , M n α , δ U n + 1 , δ V n + 1 , Δ t δ M n α ; U n α , t n α
for all ( δ Q n + 1 , δ V n + 1 , δ M n α ) and all n { 0 , 1 , N 1 } .
For a sufficiently regular problem, the method is second-order accurate for α = 1 / 2 or otherwise first-order accurate [35]. However, the method is non-convergent for α < 1 / 2 if the system of equations implied by (15) is of differential algebraic nature. Furthermore, no associated time-discrete incremental potential exists for the method if Γ exhibits a parametric dependence on U , so that the finite element systems are unsymmetric in this case.

3.3. The Modified α -Family

The “modified α -family” is a “predictor-corrector” method. In particular, the quantities Q ^ n + 1 , V ^ n + 1 , M ^ n α are introduced, which are “predictions” of Q n + 1 , V n + 1 and M n α . During the predictor step, the problem
Find the unknowns ( Q ^ n + 1 , V ^ n + 1 , M ^ n α ) such that
0 = + ( 1 α ) δ Ψ [ Q n , δ Q ^ n + 1 ] + α δ Ψ [ Q ^ n + 1 , δ Q ^ n + 1 ] + δ Γ U ^ n + 1 U n Δ t , V ^ n + 1 V n Δ t , M ^ n α , δ U ^ n + 1 , δ V ^ n + 1 , Δ t δ M ^ n α ; U n , t n α
for all ( δ Q ^ n + 1 , δ V ^ n + 1 , δ M ^ n α )
is solved. Based on the resulting predicted values, the time-discrete problem is then
Find the unknowns ( Q n + 1 , V n + 1 , M n α ) such that
0 = + ( 1 α ) δ Ψ [ Q n , δ Q n + 1 ] + α δ Ψ [ Q n + 1 , δ Q n + 1 ] + δ Γ U n + 1 U n Δ t , V n + 1 V n Δ t , M n α , δ U n + 1 , δ V n + 1 , Δ t δ M n α ; U ^ n α , t n α
for all ( δ Q n + 1 , δ V n + 1 , δ M n α ) and all n { 0 , 1 , N 1 } , where U ^ n α = ( 1 α ) U n + α U ^ n + 1 .
For many (but not all) cases, the method provides the same accuracy as the α -family, while retaining the symmetric and variational structure, so that corresponding time-discrete incremental potentials exist for both steps.

4. Implementation

In order to show the feasibility of the approach, the space–time discrete formulation described above has been preliminarily implemented in two and three dimensions for both sequential and distributed memory parallel computations. In particular, the open source C++ libraries GalerkinTools and IncrementalFE have been developed. The library GalerkinTools allows for the definition of weak forms corresponding to functionals of the generic form (7) (in fact, the library covers slightly more general weak forms including those unsymmetric weak forms required for the α -family of methods for temporal discretization). Furthermore, GalerkinTools provides methods for the automatic assembly of the corresponding finite element systems based on the functionalities of the open source finite element library deal.II (Alzetta et al. [27], Bangerth et al. [28]). In addition, interfaces for the solution of the resulting linear systems of equations are included. For this purpose, the unsymmetric sparse direct UMFPACK solver (Davis [36]), the symmetric and unsymmetric sparse direct PARDISO solvers (Schenk and Gärtner [37], Schenk and Gärtner [38], Karypis and Kumar [39]), the symmetric and unsymmetric sparse direct MUMPS solvers (Amestoy et al. [40,41]) and the symmetric sparse direct MA57 solver (HSL, a collection of Fortran codes for large-scale scientific computation, see http://www.hsl.rl.ac.uk/ (accessed on 5 January 2023)) can currently be used. Based on the library GalerkinTools, the library IncrementalFE then encapsulates the special functionality required for the solution of problems formulated within the standard dissipative framework. In particular, this involves functionalities for the definition of the functionals Ψ and Γ , together with the temporal discretization schemes described in the previous section. Furthermore, IncrementalFE implements a simple Newton–Raphson scheme, optionally combined with a line search, to solve the nonlinear finite element problems arising in each time step. For details regarding the functionalities and implementation of the libraries GalerkinTools and IncrementalFE, the reader is referred to the documentation bundled with these libraries. However, the general sequence of steps needed to define and solve a particular problem is briefly described in the following:
1.
Setup of the finite element mesh, which corresponds to the definition of the computational domain. During this step, the domain Ω is first meshed using either the built-in functionalities of deal.II or by importing a mesh from an external mesh generator. Subsequently, the domain is partitioned into domain portions by assigning “material IDs” to the individual cells of the mesh. This is followed by the definition of the interface Σ based on faces of the previously defined volume-related cells, which also goes along with the definition of the − and + sides of the interface as well as a partitioning of the interface by assigning “material IDs”. As a final step of mesh generation, local mesh refinement may be performed.
2.
Definition of the unknowns Q , V and M of the problem and the corresponding finite element spaces. During this step, the domain-related unknown fields, the interface-related unknown fields and the scalars are defined. This includes the definition of the finite elements to be used for the discretization of the domain-related and interface-related unknowns as well as Dirichlet-type constraints. Furthermore, for each unknown field, the domain portions and interface portions, respectively, are defined where the respective field may be non-zero. In order to allow for the use of vector-valued finite elements, several unknown fields of a particular type may be defined at once.
3.
Definition of the domain-related dependent fields q Ω , v Ω , m Ω and the interface-related dependent fields q Σ , v Σ , m Σ based on the previously defined unknowns Q , V and M according to Equations (5) and (6), respectively.
4.
Definition of Ψ and Γ as functionals of the form (7). It is noted in this context that the definition of the integrals over Ω and Σ , respectively, may be split into several integrals, which are summed up subsequently. This split (i) facilitates the re-usability of the implementation of certain terms (e.g., a Neo-Hookean strain energy density) and thus makes the definition of Ψ and Γ “modular”; (ii) allows for using different numerical quadrature schemes for different terms; and (iii) allows for a better structuring of the problem, which can be internally exploited by the library GalerkinTools to determine which entries in the finite element system matrix are non-zero. The actual definition of a single integral typically involves (i) the implementation of the integrand (which is expressed in terms of the dependent fields relevant to the integral) and its first and second derivatives; (ii) the specification of domain and interface portions, respectively, where the integrand is non-zero; (iii) the definition of the quadrature scheme to be used for evaluating the integral; (iv) the definition of the temporal integration scheme to be used (in fact, one may use different schemes for different terms, although this has not been discussed above); and (v) the definition of parameters (particularly constitutive parameters).
5.
Solving the problem. Currently, time-stepping needs to be implemented manually, although single time steps are automatically handled by the library IncrementalFE.
Remark 10.
  • The entire problem definition is facilitated by the use of C++ as an object-oriented language. As a result, e.g., each unknown field is a separate object, which leads to a rather natural structuring of the code in that it closely resembles the underlying mathematical statement of the problem under study in terms of formulae.
  • Typically, the most time-consuming step of the problem definition is, besides in some cases the mesh generation, the implementation of the integrands for the definition of Ψ and Γ. However, due to the mentioned modular design, the implementation of certain integrands can be re-used. This does, in principle, allow for a library of integrands covering a multitude of situations. Furthermore, automatic differentiation schemes may be used in order to simplify the implementation of the first and second derivatives of the integrands, which are required for the solution of the problem.
  • The library GalerkinTools also allows for the definition of “hidden variables” for each integrand, so that, in principle, classical rate-independent plasticity approaches are covered. For example, for the variationally consistent method, the latter may easily be achieved by locally “pre-minimizing” the time-discrete incremental potential Δ Π n n + 1 with respect to the plastic deformation variables. However, proper testing of classical plasticity approaches is yet to be performed. Therefore, a discussion of these is not included in the present contribution.
  • Implementation-wise, the main difficulties are the definition of functionals of the generic form (7) and the assembly procedures for the associated finite-element systems. These functionalities are provided by GalerkinTools; and, on this basis, the implementation of particular temporal discretization schemes and iterative algorithms for the solution of the nonlinear problems arising in each time step is straightforward. Hence, the library IncrementalFE is rather slim; and if temporal discretization schemes and/or iterative algorithms other than those discussed in this work need to be considered, these can be implemented with little effort (provided that the temporally discretized versions are compatible with the assumed form for functionals). The same applies to the use of other linear solvers than those mentioned above, and, in particular, iterative linear solvers.

5. Examples

In this section, (i) the buckling behavior of a viscoelastic dielectric elastomer subjected to an electric voltage, and (ii) the swelling/shrinking of a hydrogel subjected to a change in exterior ion concentrations are discussed as examples. It is emphasized that the main objective is to demonstrate the feasibility and practicability of the approach, and not to give a detailed exposition of particular physical effects. Therefore, the statement of the problems will be kept concise; and a discussion of the differential (algebraic) equations implied by the formulations is omitted.
It is emphasized that the discussion of the example problems mostly follows the sequence of steps described in Section 4, which highlights the systematic approach to standard dissipative problems proposed in this work. The equations are written using the symbolic vector/tensor notation; and it is tacitly assumed that each vector/tensor field is represented by an appropriate number of domain-related/interface-related scalar fields, which embody vector/tensor coordinates with regard to the standard basis of R 3 .
Moreover, it is noted that both example problems have been chosen to be axisymmetric in order to reduce the computational cost and, therefore, facilitate numerical convergence studies. Despite this fact, the problems will be described in the fully three-dimensional setting for the sake of simplicity. The steps necessary for reduction to the axisymmetric setting should be sufficiently obvious and are, therefore, not described in detail.

5.1. Example 1: Buckling of a Viscoelastic Dielectric Elastomer

The first example problem is illustrated in Figure 2. A circular disc of incompressible viscoelastic dielectric elastomer with radius L and thickness H is considered. The constitutive behavior is assumed to be described by the relations proposed by Hong [42]. The radial displacement is constrained to zero on the entire lateral surface of the disc. Furthermore, the average axial displacement is constrained to zero on the lateral surface. The latter constraint was chosen instead of a homogeneous constraint on the entire surface in order to avoid stress/strain singularities, which would complicate the convergence studies due to the need for excessive local mesh refinement. The bottom surface and the top surface are assumed to be electrodes, with a transient voltage U ( t ) applied across these electrodes. As a result of Maxwell stress, the disc shrinks in the axial direction; and due to the incompressibility of the material, this leads to buckling and large deformations of the disc in the axial direction. In order to induce the “imperfection” needed for a well-defined buckling behavior, the gravitational acceleration g is taken into account and assumed to point downwards in the figure. Fringing fields at the edges of the electrodes are neglected for simplicity.

5.1.1. Computational Domain

The computational domain Ω shown in Figure 2 represents the reference configuration of the disc, which is assumed to be stress-free. It is noted that this reference configuration is not actually assumed at any instant of time due to the action of the gravitational acceleration. The in-plane Cartesian coordinates in Figure 2 are X 1 and X 2 , while the out-of-plane coordinate is X 3 . If a radial coordinate R and an axial coordinate Z of a cylindrical coordinate system are introduced as shown in the figure, the body occupies the region described by 0 R L and H / 2 Z H / 2 . The lateral surface of this disc is Σ 1 , the bottom surface Σ 2 and the top surface Σ 3 .

5.1.2. Unknown Variables

State Variables

The thermodynamic state variables
u : Ω × T R 3 ( X , t ) u ( X , t )
D : Ω × T R 3 ( X , t ) D ( X , t )
U i : Ω × T R 6 ( X , t ) U i ( X , t )
are introduced together with the corresponding initial conditions
u ( X , t 0 ) = u 0 ( X )
D ( X , t 0 ) = D 0 ( X )
U i ( t 0 ) = I ,
with I being the identity tensor of rank 2. The thermodynamic state variables have the following significance:
  • u is the mechanical displacement field. It describes how the placements of material points in the reference configuration are mapped onto their placements in the current configuration. In particular, x ( X , t ) = X + u ( X , t ) , with x ( X , t ) being the placement of the material point X at time t.
  • D is the referential electric displacement field.
  • U i is the (symmetric) inelastic stretch tensor.
The thermodynamic state variables are supplemented by the auxiliary state variables
φ : Ω × T R ( X , t ) φ ( X , t )
p : Ω × T R ( X , t ) p ( X , t )
p i : Ω × T R ( X , t ) p i ( X , t )
l u Z : T R t l u Z ( t )
together with the corresponding initial conditions
φ ( X , t 0 ) = φ 0 ( X )
p ( X , t 0 ) = p 0 ( X )
p i ( X , t 0 ) = p 0 i ( X )
l u Z ( t 0 ) = l 0 u Z .
The auxiliary state variables have the following significance:
  • φ is the electric scalar potential, which acts as a Lagrangian multiplier enforcing the condition that there is no free space charge density in the interior of the domain.
  • p is a pressure type variable, which acts as a Lagrangian multiplier enforcing the condition that the elastomer does not exhibit any overall local volume change.
  • p i is a further pressure type variable, which acts as a Lagrangian multiplier enforcing the condition that the elastomer does not exhibit any inelastic local volume change.
  • l u Z is used to constrain the average axial displacement of the disc on the lateral surface Σ 1 to zero.
The practical determination of the initial values u 0 ( X ) , D 0 ( X ) , φ 0 ( X ) , p 0 ( X ) , p 0 i ( X ) and l 0 u Z will be discussed later.

Process Variables and Lagrangian Multipliers

For this first example, neither process variables nor Lagrangian multipliers are needed.

Dirichlet Constraints

The state variables u and D are constrained by the Dirichlet conditions
u · N = 0 on Σ 1
D · N = 0 on Σ 1 .

5.1.3. Dependent Fields

Here and in the following, only those dependent fields are listed explicitly, which are “non-trivial” in that they are not equal to an unknown field itself or the gradient/divergence/curl of an unknown field (the divergence actually involves three unknown scalar fields in three dimensions since a vector field is formally represented by an appropriate number of scalar fields; and the curl involves three unknown scalar fields and three dependent fields). In particular, the only non-trivial domain-related dependent fields used below represent the mechanical deformation gradient
F = u + I ,
where the latter equation must be expanded into coordinate-wise relations in order to be consistent with the generic form of domain-related dependent fields according to Equation (5).

5.1.4. The Functionals

As noted above, the functionals describing the variational problem are chosen such that the constitutive behavior of the viscoelastic dielectric elastomer is described by the constitutive model of Hong [42].

Helmholtz Free-Energy Functional Ψ

The Helmholtz free energy functional Ψ [ Q ] = Ψ U [ U ] + Ψ P [ P , U ] is composed of the contributions
Ψ U [ U ] = Ω 1 2 ϵ D · F · F · D + μ e 2 tr F · F 2 ln ( J ) 3          + μ i 2 tr U i 1 · F · F · U i 1 2 J J i 3 + ρ 0 g ( u ) Z d V
and
Ψ P [ P , U ] = Ω φ · D + p 1 J + p i 1 J i d V Σ 1 l u Z ( u ) Z d S .
In these equations, J = det F , J i = det U i , ( u ) Z is the axial displacement, g is the gravitational acceleration, ρ 0 is the mass density of the elastomer in the reference state, ϵ is the permittivity of the material and μ e and μ i are elastic constants.

Power Functional P

The power functional takes the form
P [ U ˙ , V ˙ ; U , t ] = Σ 3 U ( t ) D ˙ · N d S .
In the following, the ramp loading
U ( t ) = U e t t 0 t e t 0
with the maximum voltage U e will be considered.

Dissipation Functional Δ ˚

For the dissipation functional, the form
Δ ˚ [ U ˙ , V ˙ ; U , t ] = Ω η tr ( d · d ) d V
is assumed, where η is the viscosity and d = sym F · U i 1 · U i ˙ · F 1 is the inelastic stretching.

5.1.5. Reduced Formulations

Variant 1: Vector Potential Formulation

A first possibility to reduce the set of variables is to express D in terms of a vector potential A according to
D = × A .
By doing so, D becomes a dependent field, while
A : Ω × T R 3 ( X , t ) A ( X , t )
becomes a state variable. As a consequence of the identity · ( × A ) = 0 , the Lagrangian multiplier term φ · D drops out in (47), such that φ is no longer required as an unknown variable. Furthermore, the Dirichlet constraint D · N = 0 on Σ 1 is transformed to
A = A S N × e 2 on Σ 1 ,
where e 2 is the unit normal vector in the axial direction, and
A S : T R t A S ( t )
is another state variable. For details regarding this boundary condition, the reader is referred to [43]. In general, gauging is required for the uniqueness of the vector potential in three dimensions, as can also be seen in [43]. However, for the axisymmetric case considered herein, gauging can be achieved by requiring that only the circumferential component of A is non-zero, which is assumed henceforth.

Variant 2: Scalar Potential Formulation

Another possibility is to apply the procedure described in Section 2.6.2 to eliminate D . A technical difficulty occurring in this context is that the normal electric displacement D · N on the surfaces Σ 2 and Σ 3 cannot be eliminated. To resolve this issue, the additional state variable
ω : ( Σ 2 Σ 3 ) × T R ( X , t ) ω ( X , t )
is introduced, which represents D · N . Then,
Ψ U , red [ U ] = Ω ϵ 2 φ · F 1 · F · φ + μ e 2 tr F · F 2 ln ( J ) 3           + μ i 2 tr U i 1 · F · F · U i 1 2 J J i 3 + ρ 0 g ( u ) Z d V ,
Ψ P , red [ P , U ] = Ω p 1 J + p i 1 J i d V Σ 1 l u Z ( u ) Z d S Σ 2 Σ 3 φ ω d S ,
and
P [ U ˙ , V ˙ ; U , t ] = Σ 3 U ( t ) ω ˙ d S ,
while Δ ˚ remains unchanged. In this formulation, ω ˙ acts as a Lagrange multiplier enforcing the boundary conditions
φ = 0 on Σ 2
φ = U ( t ) on Σ 3 .
By imposing the latter conditions directly as Dirichlet constraints for φ , the state variable ω is eliminated from the set of unknown variables.

5.1.6. Normalization

In order to improve the conditioning of the finite element system matrices and to remove units from the simulations, it is advantageous to normalize the equations. Here, the normalization proposed by Hong [42] is used:
X ˜ = X H , t ˜ = t μ i η , ˜ ( ) = ( ) H , ( ) ˙ ˜ = ( ) ˙ η μ i , u ˜ = u H , D ˜ = D μ ϵ , U ˜ i = U i , A ˜ = A H μ ϵ , φ ˜ = φ H ϵ μ , p ˜ = p μ , p ˜ i = p i μ , l ˜ Z u = l Z u μ , t ˜ 0 = t 0 μ i η , t ˜ e = t e μ i η , U ˜ e = U e H ϵ μ , H ˜ = 1 , L ˜ = H L , μ ˜ e = μ e μ , μ ˜ i = μ i μ , ϵ ˜ = 1 , η ˜ = μ i μ , ρ ˜ 0 g ˜ = ρ 0 g H μ ,
where μ = μ e + μ i .

5.1.7. Parameters

The parameters used in the numerical example below are listed in Table 1.

5.1.8. Spatial Discretization, Temporal Discretization

The α -family with α = 1 / 2 is used for temporal discretization. The spatial discretization is done in the two-dimensional axisymmetric setting, thus reducing the number of scalar fields required for the representation of u and D to two each, and for U i to four. In the case of the vector potential A , only a single scalar field is required since only the circumferential component is non-zero. The particular quadrilateral finite elements employed for the analysis are compiled in Table 2. In this table, the nomenclature of the deal.II library [44] is used. It is noted that the use of the pair RT 1 DGQ 1 for D φ is empirical since the usual proof of stability cannot be readily extended to the axisymmetric case. In contrast, positive results concerning the stability of the pair Q 2 P 1 for u p for the axisymmetric Stokes problem are available in the literature [45]. Moreover, it is noted that the support points of the elements used for U i and p i have been chosen to coincide with the quadrature points; and, as a consequence, these variables become “local” variables similar to the plastic strain in classical elastoplasticity (in fact, these variables can be eliminated at the quadrature point level, although this was not performed herein).

5.1.9. Initial Values

As mentioned earlier, the initial values for the state variables and auxiliary state variables need to be determined before the calculation. Due to the action of the gravitational acceleration, the initial state is inhomogeneous and must, therefore, be determined numerically. This has for simplicity been achieved by taking a very small time step of 1 · 10 8 t ˜ e with α = 1 before the actual computation, such that any increments in U i remain negligibly small. A more rigorous approach would be to explicitly constrain U ˙ i to zero for determining the initial values. This has, however, not been implemented yet.

5.1.10. Example Calculations and Results

The implementation of the constitutive equations and the numerical procedures were first validated against the results obtained by Hong [42], as can be seen in Appendix A. Subsequently, spatial and temporal convergence studies have been performed for the formulations in terms of D φ (“mixed formulation”), A (“vector potential formulation”) and φ (“scalar potential formulation”). The starting point for these simulations is to take N t = 16 equally spaced time increments for the entire simulation together with a coarse finite element mesh consisting of five rectangular cells in the radial direction and a single rectangular cell in the axial direction, see Figure 3 for an illustration of the coarse finite element mesh. Then, the number of time increments is increased to N t = 16 · 2 m t , where m t is the refinement cycle. Similarly, the finite element mesh is uniformly refined m h times (in this work, the term “uniform refinement” refers to splitting a rectangular cell into four equally shaped rectangular cells by introducing new vertices at the midpoints of the edges and at the center of the cell). In particular, calculations with m t = 0 , 1 , , m t , max ( m t , max = 6 ) and m h = 0 were performed to study the convergence behavior in time; and calculations with m h = 0 , 1 , , m h , max ( m h , max = 6 ) and m t = 0 were performed to study the spatial convergence behavior. For the finest mesh with m h = 6 , this results in 1 , 231 , 108 degrees of freedom for the scalar potential formulation, 1 , 395 , 075 degrees of freedom for the mixed formulation, and 1 , 231 , 109 degrees of freedom for the vector potential formulation. However, it is emphasized that the computational cost cannot be directly related to the number of unknowns, since the system structure as well as the linear solver employed also play a major role in this regard; and the study of these aspects is beyond the scope of the present contribution.
With regard to the convergence behavior in time, the error e m t is, for each m t = 0 , 1 , , m t , max 1 , evaluated according to
e m t = ˜ u ˜ m t ˜ u ˜ m t , max L 2 .
For simplicity, no other unknowns are included into the error. However, it has also been checked that the other unknowns show the expected convergence behavior.
An analogous approach is utilized to assess the spatial convergence behavior. In particular, the error e m h is, for each m h = 0 , 1 , , m h , max 1 , evaluated according to
e m h = ˜ u ˜ m h ˜ u ˜ m h , max L 2 .
The deformed shape of the disc is shown in Figure 3 together with the radial inelastic stretch. It can be noticed that substantial inelastic deformation of approximately 10 % takes place at the lateral surface of the body. Figure 4 shows the results of the convergence studies. In particular, Figure 4a shows the temporal convergence behavior, while Figure 4b shows the spatial convergence behavior. It is seen that in both cases the expected rates of convergence of k t = 2 and k h = 2 are observed for all three formulations, which confirms the validity of the approach. It is interesting that all three formulations essentially exhibit the same convergence behavior, a fact which is likely to be attributed to the relatively simple electric field distribution in the disc.

5.2. Example 2: Swelling/Shrinking of a Hydrogel

The second example problem is illustrated in Figure 5. A circular cylinder of hydrogel material with radius H and height H is considered. The hydrogel is assumed to consist of a polymeric backbone material with charged groups attached to it, and an electrolyte composed of water and the ions of a fully dissociated monovalent salt. The polymeric backbone material and water are assumed to be incompressible, and the volume change associated with a change in the concentrations of cations and anions is neglected. It is further assumed that there are no voids in the hydrogel and, hence, the volume fractions of polymeric backbone material and water must locally add up to 1. The constitutive behavior of the hydrogel is assumed to be described by the relations proposed by Acartürk [4], with the simplification of local electroneutrality (i.e., the local net space charge density is assumed to be zero). The displacement is fixed to zero at the center of the bottom surface of the cylinder, while all other surfaces are mechanically free. No matter is allowed to flow across the bottom and lateral surfaces. In contrast, water and ions may be exchanged across the top surface, which is in contact with an external electrolytic solution. The latter is composed of water and a fully dissociated monovalent salt of the same type as in the hydrogel, with a prescribed transient homogeneous molar concentration of both cations and anions, of c ¯ ext ( t ) . Starting from equilibrium, the external ion concentration is changed and the response of the hydrogel in terms of swelling/shrinking deformation and change in local water and ion concentrations is tracked.

5.2.1. Computational Domain

The computational domain Ω shown in Figure 5 represents the reference configuration of the cylinder, which coincides with the configuration at the initial time t = t 0 . The in-plane Cartesian coordinates in Figure 2 are again X 1 and X 2 , while the out-of-plane coordinate is X 3 . In terms of the radial coordinate R and the axial coordinate Z of a cylindrical coordinate system, the body occupies the region described by 0 R H and 0 Z H . The lateral surface of the cylinder is Σ 1 , the bottom surface Σ 2 , and the top surface Σ 3 .

5.2.2. Unknown Variables

State Variables

The thermodynamic state variables
u : Ω × T R 3 ( X , t ) u ( X , t ) .
c 2 H O : Ω × T R ( X , t ) c 2 H O ( X , t )
c i : Ω × T R ( X , t ) c i ( X , t )
are introduced, where i { A + , B } . The corresponding initial conditions read
u ( X , t 0 ) = 0
c 2 H O ( X , t 0 ) = c 2 H O 0
c i ( X , t 0 ) = c 0 i ,
with c 0 H 2 O and c 0 i being constants since the initial state is homogeneous. The thermodynamic state variables have the following significance:
  • u is the mechanical displacement field.
  • c 2 H O is the referential molar concentration of water (with regard to the total volume in the reference state).
  • c i are the referential molar concentrations of the A + and B ions (with regard to the total volume in the reference state).
These thermodynamic state variables are supplemented by the auxiliary state variables
p : Ω × T R ( X , t ) p ( X , t )
φ : Ω × T R ( X , t ) φ ( X , t )
together with the corresponding initial conditions
p ( X , t 0 ) = p 0
φ ( X , t 0 ) = φ 0 ,
where p 0 and φ 0 are constants. The auxiliary state variables have the following significance:
  • p is a pressure type variable, which acts as a Lagrangian multiplier enforcing the condition that the volume fractions of water and polymeric backbone material must add up to 1 everywhere.
  • φ is the electric scalar potential, which enforces the local electroneutrality condition.
It is noted that the constants c 2 H O 0 , c 0 i , p 0 , and φ 0 cannot be chosen independently. Rather, they are interrelated through the assumption that the initial state at time t = t 0 is equilibrated.

Process Variables

The process variables involved in the problem are
I 2 H O : Ω × T R 3 ( X , t ) I 2 H O ( X , t )
I i : Ω × T R 3 ( X , t ) I i ( X , t )
These have the following significance:
  • I 2 H O is the time-integrated referential flux of water in the hydrogel.
  • I i are the time-integrated referential fluxes of the ionic species i { A + , B } .

Lagrangian Multipliers

The Lagrangian multipliers
η 2 H O : Ω × ( T t 0 ) R ( X , t ) η 2 H O ( X , t )
η i : Ω × ( T t 0 ) R ( X , t ) η i ( X , t )
φ s : Σ 3 × ( T t 0 ) R ( X , t ) φ s ( X , t )
are introduced. These have the following significance:
  • η 2 H O and η i are used to incorporate the respective balance equations for water and the ionic species A + and B .
  • φ s is used to ensure that no net charge flows locally across Σ 3 . It is noted in this context that the assumption of no local net charge flow implies a certain behavior of the external solution bath. Other assumptions are possible (e.g., no global net charge flow), however, these are in any case only an approximation to the real behavior; and for the proper resolution of the behavior at the boundary, the exterior solution bath needs to be explicitly modeled.

Dirichlet Constraints

The unknowns u , I 2 H O , I i and η B are subject to the following Dirichlet type constraints
u = 0 at X = ( 0 , 0 , 0 )
I 2 H O · N = 0 on Σ 1 Σ 2
I i · N = 0 on Σ 1 Σ 2
η B = η 0 B at X = ( 0 , 0 , 0 ) .
Here, η B is constrained to the arbitrary constant η 0 B at the origin in order to ensure the uniqueness of the involved scalar potential type quantities, which are only determined up to a constant by the model.

5.2.3. Dependent Fields

As before, only the “non-trivial”-dependent fields are listed; and the mechanical deformation gradient
F = u + I
is the only field of this type.

5.2.4. The Functionals

Helmholtz Free Energy Functional Ψ

The part Ψ U of the Helmholtz free energy functional is assumed to take the form
Ψ U [ U ] = Ω ψ el ( F ) d V + Ω c 2 H O ψ fl c i c 2 H O d V ,
with the first integral representing the elastic-free energy contribution of the polymeric backbone and the second integral the free energy contributions of the fluid. In particular,
ψ el = μ 2 tr F · F 2 ln ( J ) 3 + λ ( 1 n 0 ) 2 J 1 1 n 0 ln J n 0 1 n 0 ,
ψ fl = i c i c 2 H O R T ln c i c 2 H O 1
are assumed. Here, J = det F again, λ and μ are Lamé’s constants, n 0 = 1 c 2 H O 0 V 2 H O m is the volume fraction of polymeric backbone material in the reference state, where V 2 H O m denotes the molar volume of water, R is the gas constant and T is the temperature.
The part Ψ P of the Helmholtz free energy functional reads
Ψ P [ P , U ] = Ω p n 0 + c 2 H O V 2 H O m J d V + Ω F φ i z i c i + z c d V ,
where the first integral enforces the condition that the volume fractions of fluid and polymeric backbone material add up to 1, and the second integral incorporates the electroneutrality condition. In the latter term, F is Faraday’s constant, z i and z are the respective charges of the ionic species A + , B and the charged groups attached to the backbone polymer in multiples of the elementary charge, and c is the (constant) referential molar concentration of the charged groups.

Power Functional P

For the power functional, the form
P [ U ˙ , V ˙ ; U , t ] = Σ 3 μ ¯ 2 H O ( t ) I ˙ 2 H O · N + i μ ¯ i ( t ) I ˙ i · N d S
is assumed, where
μ ¯ 2 H O ( t ) = 2 R T c ¯ ext ( t ) V 2 H O m
μ ¯ i ( t ) = R T ln c ¯ ext ( t ) V 2 H O m .
For the purposes of this example, the ramp function
c ¯ ext ( t ) = c ¯ 0 ext + c ¯ e ext c ¯ 0 ext t t 0 t e t 0
is used for c ¯ ext ( t ) , where c ¯ 0 ext is the initial concentration of ions in the external solution at time t = t 0 , and c ¯ e ext the final concentration at time t = t e .

Lagrangian Multiplier Functional L

The Lagrangian multiplier functional is assumed to take the form
L [ M , U ˙ , V ˙ ; U , t ] = Ω η 2 H O · I ˙ 2 H O c ˙ 2 H O d V + Ω i η i · I ˙ i c ˙ i d V + Σ 3 F φ s i z i I ˙ i · N d S .
where the first two integrals incorporate the respective balance equations for water and the two ion species, and the third integral enforces the condition that the local net charge flow across Σ 3 is zero.

Dissipation Functional Δ ˚

For the dissipation functional, the form
Δ ˚ [ U ˙ , V ˙ ; U , t ] = Ω δ ˚ I ˙ 2 H O , I ˙ i ; F , c i , c 2 H O d V
is assumed. In this relation,
δ ˚ = c 2 H O V 2 H O m J i R T 2 D i c i I ˙ i c i c 2 H O I ˙ 2 H O · F · F · I ˙ i c i c 2 H O I ˙ 2 H O = c 2 H O V 2 H O m det F + R T 2 D 2 H O c 2 H O I ˙ 2 H O · F · F · I ˙ 2 H O
is the local dissipation function, where D i are the respective diffusivities of the ionic species in water, and D 2 H O is a material parameter characterizing the Darcy type “resistance” associated with the motion of water through the polymeric backbone material.
For later use, the dual local dissipation function is introduced according to
ϕ ( η 2 H O , η i ; F , c i , c 2 H O ) = inf I ˙ 2 H O , I ˙ i η 2 H O · I ˙ 2 H O + i η i · I ˙ i + δ ˚ I ˙ 2 H O , I ˙ i ; F , c i , c 2 H O = J c 2 H O V 2 H O m i D i c i 2 R T η i · F 1 · F · η i + D 2 H O c 2 H O 2 R T = × η 2 H O + i c i c 2 H O η i · F 1 · F · η 2 H O + i c i c 2 H O η i .

5.2.5. Reduced Formulation

By using the procedure described in Section 2.6.1, the variables I 2 H O and I i can be eliminated on Ω . This leads to the reduced sets of variables
V 0 = I 2 H O N I N i M 0 = M ,
where
I 2 H O N : Σ 3 × T R ( X , t ) I 2 H O N ( X , t )
I N i : Σ 3 × T R ( X , t ) I N i ( X , t )
represent I 2 H O · N and I i · N , respectively, on Σ 3 . The corresponding reduced functional is obtained from Equation (16) using the dual local dissipation function (95). It reads
Γ red [ U ˙ , V ˙ 0 , M 0 ; U , t ] = Ω η 2 H O c ˙ 2 H O + i η i c ˙ i d V Ω ϕ ( η 2 H O , η i ; F , c i , c 2 H O ) d V Σ 3 η 2 H O μ ¯ 2 H O ( t ) I ˙ 2 H O N d S Σ 3 i η i F φ s z i μ ¯ i ( t ) I ˙ N i d S .
Similarly to example 1, I 2 H O N and I N i can be eliminated by imposing appropriate Dirichlet conditions for η 2 H O and η i . However, due to the involvement of φ s in these conditions, the resulting equations are more complicated to implement and hence, I 2 H O N and I N i are kept as variables.

5.2.6. Normalization

All quantities are normalized using the reference values c * = 1 mol / L , D * = 1 · 10 4   cm 2 / s , R * = 8.3144   J / ( mol K ) , T * = 298   K , F * = 96 , 485.33   C / mol and L * = 1 mm . In particular, the normalized quantities are
X ˜ = X L * , t ˜ = t D * L * 2 , ˜ ( ) = ( ) L * , ( ) ˙ ˜ = ( ) ˙ L * 2 D * , u ˜ = u L * , c ˜ 2 H O = c 2 H O c * , c ˜ i = c i c * , p ˜ = p R * T * c * , φ ˜ = φ F * R * T * , I ˜ 2 H O = I 2 H O c * L * , I ˜ i = I i c * L * , I ˜ 2 H O N = I 2 H O N c * L * , I ˜ N i = I N i c * L * , η ˜ 2 H O = η 2 H O R * T * , η ˜ i = η i R * T * , φ ˜ s = φ s F * R * T * , t ˜ 0 = t 0 D * L * 2 , t ˜ e = t e D * L * 2 , c ¯ ˜ 0 ext = c ¯ 0 ext c * , c ¯ ˜ e ext = c ¯ e ext c * , H ˜ = H L * , c ˜ 0 H 2 O = c 0 H 2 O c * , c ˜ 0 i = c 0 i c * , p ˜ 0 = p 0 R * T * c * , φ ˜ 0 = φ 0 F * R * T * , μ ˜ = μ R * T * c * , λ ˜ = λ R * T * c * , V ˜ 2 H O m = V 2 H O m c * , n ˜ 0 = n 0 , z ˜ i = z i , z ˜ = z , c ˜ = c c * , D ˜ 2 H O = D 2 H O D * , D ˜ i = D i D * , R ˜ = R R * , F ˜ = F F * , T ˜ = T T * .

5.2.7. Parameters

The set of parameters used below is listed in Table 3. The constitutive parameters have been taken from Table 5.1 in Acartürk [4] and normalized. It is noted in this context that some of the parameters are defined differently in the latter work than in the present work, so that some conversion between the parameters is involved. This applies in particular to D 2 H O , which is equal to R T k F / ( γ F R V 2 H O m ) in the work of Acartürk [4] (the quantities k F and γ F R are defined in [4]).
The initial values appearing in Table 3 were computed from the condition that the initial state is an equilibrium state, which leads to the following equations
c 2 H O 0 = 1 n 0 V 2 H O m
c 0 A + = 1 2 z c z A + + 1 2 | z A + | z c 2 z A + z B 2 c ¯ 0 ext c 2 H O 0 V 2 H O m 2
c 0 B = 1 2 z c z B + 1 2 | z B | z c 2 z A + z B 2 c ¯ 0 ext c 2 H O 0 V 2 H O m 2
p 0 = R T V 2 H O m c 0 A + + c 0 B c 2 H O 0 2 V 2 H O m c ¯ 0 ext
φ 0 = R T F z A + ln c ¯ 0 ext c 2 H O 0 V 2 H O m c 0 A + .
The corresponding normalized values given in Table 3 are subject to rounding. However, in the simulations, values accurate to within numerical accuracy have been used. This also applies to λ ˜ , μ ˜ and D ˜ 2 H O .

5.2.8. Spatial Discretization and Temporal Discretization

Again, the α -family with α = 1 / 2 is used for temporal discretization, and the spatial discretization is performed in the two-dimensional axisymmetric setting. The particular quadrilateral finite elements employed for the analysis are compiled in Table 4. No elements are given for I 2 H O and I i in this table because only the reduced formulation is considered in the numerical example below. It is noted that there are no interface terms coupling the unknowns associated with the discontinuous finite elements of type DGP 1 between neighboring cells. This means that the unknowns associated with these elements are local to the cell and can be eliminated with little cost by the direct linear solver. Another possibility would be to use discontinuous Lagrange elements of type DGQ 2 for c 2 H O , c i and φ with the support points coinciding with the quadrature points. However, to fully benefit from this approach where the respective unknowns are local to the quadrature points, changes to the implementation would be necessary, which have not been implemented yet.

5.2.9. Example Calculations and Results

The implementation of the constitutive equations and the numerical procedures was first validated against a result obtained by Acartürk [4], as can be seen in Appendix B. Subsequently, spatial and temporal convergence studies have been performed for the reduced formulation. In the latter case, the starting point is to take N t = 4 equally spaced time increments for the entire simulation together with a coarse finite element mesh consisting of a single cell. With regard to the definition and evaluation of errors, exactly the same approach as in Section 5.1.10 is used. However, for this example, the convergence behavior in time is studied with m t , max = 8 and m h = 3 , while the spatial convergence study is conducted with m h , max = 7 and m t = 2 . In the latter case, this results in 577 , 033 unknowns for the finest mesh with m h = 7 .
The deformed shape of the hydrogel cylinder is shown in Figure 6 together with the cation concentration c ˜ A + and the outline of the undeformed cylinder; whilst Figure 7a,b show the results of the temporal and spatial convergence study, respectively. The expected rates of convergence of k t = 2 and k h = 2 are clearly observed as for the first example.

6. Concluding Remarks

In the present contribution, a systematic approach to the thermodynamically consistent modeling of isothermal dissipative continuum problems using the concept of standard dissipative continua is proposed and applied to two example problems. The formulation does explicitly allow for different spatial regions, which are associated with different physics, as well as for the description of interfacial processes, which may be coupled to the physics in the volume. An important conceptual feature of the approach is the distinction between state variables, process variables and Lagrangian multipliers. While the state variables have the usual significance in that they describe the thermodynamic state of the system under consideration, the process variables are used to constrain the ways in which the thermodynamic state can change; and the latter constraints are incorporated with the help of the Lagrangian multipliers. Due to the variational structure of the formulation, a systematic elimination of variables is possible under certain circumstances, thus reducing the number of unknown variables, which is an important aspect for numerical solution procedures.
The variational approach is also beneficial in view of the temporal and spatial discretization of the space–time continuous problem. With regard to temporal discretization, it is, at least for a certain class of problems, possible to use generic discretization approaches with known convergence properties. Similarly, when using the finite element method for spatial discretization, much of the classical theory applies.
Albeit a large amount of formalization is clearly involved in the presented formulation, it is expected that the associated systematic approach substantially reduces the time penalty associated with implementing complex multiphysics models. At the same time, it likely reduces the potential for error both, during the formulation of multiphysics models as well as during the implementation into program code. In particular, it appears possible to solve a large class of problems using a single, well-tested implementation.
Some of the above aspects were demonstrated by considering two example problems. In this context, it was empirically shown that optimal rates of convergence in space and time can be achieved. For the purposes of this work, the example problems were kept simple. However, applications to more complex problems will be described in forthcoming contributions.

Funding

Financial support from the Deutsche Forschungsgemeinschaft (DFG) under Grants STA 1593/1-1 and STA 1593/2-1 is gratefully acknowledged.

Data Availability Statement

The source code for the library GalerkinTools and the corresponding documentation are available from the following repository: https://github.com/sebastian-stark/GalerkinTools (accessed on 26 February 2023); The source code for the library IncrementalFE and the corresponding documentation are available from the following repository: https://github.com/sebastian-stark/IncrementalFE (accessed on 26 February 2023); The source code used for the second example is available from the following repository: https://github.com/sebastian-stark/hydrogel_model_concentration_change (accessed on 26 February 2023); The source code used for the first example is in part available from the author upon request (the full source code cannot be made publicly available due to the use of a software library, which cannot be made available for legal reasons).

Acknowledgments

The author wishes to thank B.D. Reddy for an invitation for a research visit at the Centre for Research in Computational and Applied Mechanics at the University of Cape Town, during which part of the work underlying the present contribution was conducted.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A. Validation of the Implementation for Example 1

In order to verify the correct implementation of the constitutive relations for example 1, some results were compared to those obtained by Hong [42], with Figure A1 showing a representative example. In particular, a laterally free plate of viscoelastic dielectric elastomer is considered, with the fringing fields neglected. An electric voltage is applied across the electrodes on the bottom and top faces of this plate. The resulting nominal electric field in the plate is E ˜ = φ ˜ / H ˜ , with H ˜ = 1 being the normalized thickness of the plate. The nominal electric field is increased at a rate of d E ˜ / d t ˜ = 0.2 until instability takes place. The material parameters used by Hong [42] coincide with the parameters given in Table 1 (note that the definition of the viscosity used by Hong [42] differs by a factor of 2 from the one used in the present work). For the comparisons, the total time interval was split into 1024 equally spaced time increments, and the α -family with α = 1 / 2 was used for temporal discretization. Since the fields are homogeneous, the use of a single finite element is sufficient. Figure A1 shows the results for the total lateral stretch λ and the inelastic lateral stretch λ i . It is seen that the results obtained with all three formulations are indistinguishable; and these results are consistent with the results of Hong [42], which confirms the correct implementation for at least the case of homogeneous uniaxial fields.
Figure A1. Comparison between the present formulation and the results obtained by Hong [42].
Figure A1. Comparison between the present formulation and the results obtained by Hong [42].
Axioms 12 00267 g0a1

Appendix B. Validation of the Implementation for Example 2

In order to verify the correct implementation of the constitutive relations for example 2, an example problem described by Acartürk [4] is considered. The problem involves a block of hydrogel material with a width of 0.2 mm and a height of 1 mm in a plane strain setting. The normal mechanical displacement is constrained to zero on the lateral surfaces of the block as well as on its bottom surface, while the top surface is mechanically free. No matter is allowed to flow across the lateral surfaces and the bottom surfaces, while the top surface is in contact with a solution bath in the same way as shown in Figure 5. The ion concentration is linearly ramped from c ¯ ˜ ext = 0.9 to c ¯ ˜ ext = 0.1 between time t ˜ 0 = 0 and t ˜ 1 = 0.5 . Subsequently, the external ion concentration is kept constant until the end of the simulation at time t ˜ e = 200 . In the present work, the calculation was performed with 80 finite elements in the axial direction. For the time stepping, an initial grid was created with 8 equally spaced time increments for the linear ramp and another 8 time increments for the subsequent step with constant external concentration. In the case of the latter, no equal time increment spacing was chosen. Rather, the time increment size was increased by a factor of 3 between subsequent time increments. For the actual computation, this initial temporal grid was uniformly refined six times, thus resulting in 1024 time increments in total. The resulting numerical solution is compared to the variant (SIb) discussed in Acartürk [4]. For this variant, the net electric current is constrained to zero; and for the particular one-dimensional problem with a monovalent salt considered herein, this is equivalent to the requirement of local electroneutrality used in this work. It is noted in this context that a number of different formulations were discussed by Acartürk [4], which should theoretically be equivalent. However, different results were obtained by Acartürk [4]; and this suggests that there are either errors in the implementation or other numerical issues for some of the formulations. The choice to consider the variant (SIb) for comparison is motivated by the fact that the results of this variant agree very well with the results of the present formulation. In particular, Figure A2 shows a comparison of the normal displacement u ˜ Z of the top surface of the hydrogel block over time. Excellent agreement is evident for this case. Figure A3 shows comparisons for the cation concentrations (with regard to the actual volume of the electrolyte) at different instants of time. Again, good agreement is observed.
Figure A2. Comparison of the results for the normal displacement of the top surface of the hydrogel block [4].
Figure A2. Comparison of the results for the normal displacement of the top surface of the hydrogel block [4].
Axioms 12 00267 g0a2
Figure A3. Comparison of the results for the cation concentration along the axial coordinate Z ˜ at different times [4].
Figure A3. Comparison of the results for the cation concentration along the axial coordinate Z ˜ at different times [4].
Axioms 12 00267 g0a3

References

  1. Landis, C. Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics. J. Mech. Phys. Solids 2002, 50, 127–152. [Google Scholar] [CrossRef]
  2. Miehe, C.; Rosato, D.; Kiefer, B. Variational principles in dissipative electro-magneto-mechanics: A framework for the macro-modeling of functional materials. Int. J. Numer. Methods Eng. 2011, 86, 1225–1276. [Google Scholar] [CrossRef]
  3. Bucci, G.; Chiang, Y.M.; Carter, W.C. Formulation of the coupled electrochemical-mechanical boundary-value problem, with applications to transport of multiple charged species. Acta Mater. 2016, 104, 33–51. [Google Scholar] [CrossRef] [Green Version]
  4. Acartürk, A.Y. Simulation of Charged Hydrated Porous Materials. Ph.D. Thesis, University of Stuttgart, Stuttgart, Germany, 2009. [Google Scholar]
  5. Simo, J.C.; Honein, T. Variational Formulation, Discrete Conservation Laws, and Path-Domain Independent Integrals for Elasto-Viscoplasticity. J. Appl. Mech. 1990, 57, 488–497. [Google Scholar] [CrossRef]
  6. Ortiz, M.; Stainier, L. The variational formulation of viscoplastic constitutive updates. Comput. Methods Appl. Mech. Eng. 1999, 171, 419–444. [Google Scholar] [CrossRef]
  7. Hackl, K.; Fischer, F.D. On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2008, 464, 117–132. [Google Scholar] [CrossRef]
  8. Yang, Q.; Stainier, L.; Ortiz, M. A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids. J. Mech. Phys. Solids 2006, 54, 401–424. [Google Scholar] [CrossRef]
  9. Carstensen, C.; Hackl, K.; Mielke, A. Non-convex potentials and microstructures in finite-strain plasticity. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2002, 458, 299–317. [Google Scholar] [CrossRef]
  10. Ortiz, M.; Repetto, E.A. Nonconvex energy minimization and dislocation structures in ductile single crystals. J. Mech. Phys. Solids 1999, 47, 397–462. [Google Scholar] [CrossRef]
  11. Halphen, B.; Nguyen, Q.S. Sur les matériaux standard généralisés. J. De Mécanique 1975, 14, 39–63. [Google Scholar]
  12. Germain, P.; Nguyen, Q.S.; Suquet, P. Continuum Thermodynamics. J. Appl. Mech. 1983, 50, 1010–1020. [Google Scholar] [CrossRef]
  13. Biot, M.A. Theory of Finite Deformations of Porous Solids. Indiana Univ. Math. J. 1972, 21, 597–620. [Google Scholar] [CrossRef]
  14. Miehe, C. Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation. Int. J. Numer. Methods Eng. 2002, 55, 1285–1322. [Google Scholar] [CrossRef]
  15. Miehe, C.; Schotte, J.; Lambrecht, M. Homogenization of inelastic solid materials at finite strains based on incremental minimization principles. Application to the texture analysis of polycrystals. J. Mech. Phys. Solids 2002, 50, 2123–2167. [Google Scholar] [CrossRef]
  16. Miehe, C. A multi-field incremental variational framework for gradient-extended standard dissipative solids. J. Mech. Phys. Solids 2011, 59, 898–923. [Google Scholar] [CrossRef]
  17. Miehe, C. Variational gradient plasticity at finite strains. Part I: Mixed potentials for the evolution and update problems of gradient-extended dissipative solids. Comput. Methods Appl. Mech. Eng. 2014, 268, 677–703. [Google Scholar] [CrossRef]
  18. Miehe, C.; Welschinger, F.; Aldakheel, F. Variational gradient plasticity at finite strains. Part II: Local-global updates and mixed finite elements for additive plasticity in the logarithmic strain space. Comput. Methods Appl. Mech. Eng. 2014, 268, 704–734. [Google Scholar] [CrossRef]
  19. Miehe, C.; Mauthe, S.; Hildebrand, F. Variational gradient plasticity at finite strains. Part III: Local-global updates and regularization techniques in multiplicative plasticity for single crystals. Comput. Methods Appl. Mech. Eng. 2014, 268, 735–762. [Google Scholar] [CrossRef]
  20. McBride, A.T.; Reddy, B.D.; Steinmann, P. Dissipation-consistent modelling and classification of extended plasticity formulations. J. Mech. Phys. Solids 2018, 119, 118–139. [Google Scholar] [CrossRef]
  21. Miehe, C.; Rosato, D. A rate-dependent incremental variational formulation of ferroelectricity. Int. J. Eng. Sci. 2011, 49, 466–496. [Google Scholar] [CrossRef]
  22. Miehe, C.; Kiefer, B.; Rosato, D. An incremental variational formulation of dissipative magnetostriction at the macroscopic continuum level. Int. J. Solids Struct. 2011, 48, 1846–1866. [Google Scholar] [CrossRef]
  23. Miehe, C.; Mauthe, S.; Teichtmeister, S. Minimization principles for the coupled problem of Darcy-Biot-type fluid transport in porous media linked to phase field modeling of fracture. J. Mech. Phys. Solids 2015, 82, 186–217. [Google Scholar] [CrossRef]
  24. Miehe, C.; Hildebrand, F.E.; Böger, L. Mixed variational potentials and inherent symmetries of the Cahn-Hilliard theory of diffusive phase separation. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2014, 470, 20130641. [Google Scholar] [CrossRef] [Green Version]
  25. Miehe, C.; Mauthe, S.; Ulmer, H. Formulation and numerical exploitation of mixed variational principles for coupled problems of Cahn-Hilliard-type and standard diffusion in elastic solids. Int. J. Numer. Methods Eng. 2014, 99, 737–762. [Google Scholar] [CrossRef]
  26. Miehe, C.; Dal, H.; Schänzel, L.M.; Raina, A. A phase-field model for chemo-mechanical induced fracture in lithium-ion battery electrode particles. Int. J. Numer. Methods Eng. 2016, 106, 683–711. [Google Scholar] [CrossRef]
  27. Alzetta, G.; Arndt, D.; Bangerth, W.; Boddu, V.; Brands, B.; Davydov, D.; Gassmoeller, R.; Heister, T.; Heltai, L.; Kormann, K.; et al. The deal.II Library, Version 9.0. J. Numer. Math. 2018, 26, 173–183. [Google Scholar] [CrossRef]
  28. Bangerth, W.; Hartmann, R.; Kanschat, G. deal.II—A General Purpose Object Oriented Finite Element Library. ACM Trans. Math. Softw. 2007, 33, 24/1–24/27. [Google Scholar] [CrossRef] [Green Version]
  29. Biot, M.A. Theory of Stress-Strain Relations in Anisotropic Viscoelasticity and Relaxation Phenomena. J. Appl. Phys. 1954, 25, 1385–1391. [Google Scholar] [CrossRef]
  30. Biot, M.A. Variational Principles in Irreversible Thermodynamics with Application to Viscoelasticity. Phys. Rev. 1955, 97, 1463–1469. [Google Scholar] [CrossRef] [Green Version]
  31. Abed-Meraim, F.; Nguyen, Q.S. A quasi-static stability analysis for Biot’s equation and standard dissipative systems. Eur. J. Mech.—A/Solids 2007, 26, 383–393. [Google Scholar] [CrossRef] [Green Version]
  32. Rockafellar, R.T. Convex Analysis; Princeton University Press: Princeton, NJ, USA, 1970. [Google Scholar]
  33. Bathe, K.J. The inf−sup condition and its evaluation for mixed finite element methods. Comput. Struct. 2001, 79, 243–252. [Google Scholar] [CrossRef] [Green Version]
  34. Auricchio, F.; da Veiga, L.B.; Brezzi, F.; Lovadina, C. Mixed Finite Element Methods. In Encyclopedia of Computational Mechanics Second Edition; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 2017; pp. 1–53. [Google Scholar] [CrossRef] [Green Version]
  35. Stark, S. On a certain class of one step temporal integration methods for standard dissipative continua. Comput. Mech. 2021, 67, 265–287. [Google Scholar] [CrossRef]
  36. Davis, T.A. Algorithm 832: UMFPACK V4.3—An Unsymmetric-pattern Multifrontal Method. ACM Trans. Math. Softw. 2004, 30, 196–199. [Google Scholar] [CrossRef]
  37. Schenk, O.; Gärtner, K. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Gener. Comput. Syst. 2004, 20, 475–487. [Google Scholar] [CrossRef]
  38. Schenk, O.; Gärtner, K. On fast factorization pivoting methods for symmetric indefinite systems. Electron. Trans. Numer. Anal. 2006, 23, 158–179. [Google Scholar]
  39. Karypis, G.; Kumar, V. A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM J. Sci. Comput. 1998, 20, 359–392. [Google Scholar] [CrossRef]
  40. Amestoy, P.R.; Duff, I.S.; Koster, J.; L’Excellent, J.Y. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling. SIAM J. Matrix Anal. Appl. 2001, 23, 15–41. [Google Scholar] [CrossRef] [Green Version]
  41. Amestoy, P.R.; Guermouche, A.; L’Excellent, J.Y.; Pralet, S. Hybrid scheduling for the parallel solution of linear systems. Parallel Comput. 2006, 32, 136–156. [Google Scholar] [CrossRef] [Green Version]
  42. Hong, W. Modeling viscoelastic dielectrics. J. Mech. Phys. Solids 2011, 59, 637–650. [Google Scholar] [CrossRef]
  43. Stark, S.; Semenov, A.; Balke, H. On the boundary conditions for the vector potential formulation in electrostatics. Int. J. Numer. Methods Eng. 2015, 102, 1704–1732. [Google Scholar] [CrossRef]
  44. deal.II Documentation. Available online: https://www.dealii.org/current/doxygen/deal.II/ (accessed on 28 January 2023).
  45. Ruas, V. Mixed finite element methods with discontinuous pressures for the axisymmetric Stokes problem. ZAMM—J. Appl. Math. Mech./Z. Für Angew. Math. Und Mech. 2003, 83, 249–264. [Google Scholar] [CrossRef]
Figure 1. Domain under consideration.
Figure 1. Domain under consideration.
Axioms 12 00267 g001
Figure 2. Illustration of example problem 1.
Figure 2. Illustration of example problem 1.
Axioms 12 00267 g002
Figure 3. Coarse finite element mesh in the reference configuration together with radial inelastic stretch in the deformed configuration at t = t e (result obtained with m t = 0 and m h = 6 ).
Figure 3. Coarse finite element mesh in the reference configuration together with radial inelastic stretch in the deformed configuration at t = t e (result obtained with m t = 0 and m h = 6 ).
Axioms 12 00267 g003
Figure 4. Convergence results for example 1: (a) temporal convergence; and (b) spatial convergence.
Figure 4. Convergence results for example 1: (a) temporal convergence; and (b) spatial convergence.
Axioms 12 00267 g004
Figure 5. Illustration of example problem 2.
Figure 5. Illustration of example problem 2.
Axioms 12 00267 g005
Figure 6. Deformed shape of the hydrogel cylinder at time t ˜ e together with cation concentration c ˜ A + and outline of the undeformed shape (result obtained with m t = 2 and m h = 7 ).
Figure 6. Deformed shape of the hydrogel cylinder at time t ˜ e together with cation concentration c ˜ A + and outline of the undeformed shape (result obtained with m t = 2 and m h = 7 ).
Axioms 12 00267 g006
Figure 7. Convergence results for example 2: (a) temporal convergence, (b) spatial convergence.
Figure 7. Convergence results for example 2: (a) temporal convergence, (b) spatial convergence.
Axioms 12 00267 g007
Table 1. Parameters used for example 1.
Table 1. Parameters used for example 1.
QuantityValue
t ˜ 0 0
t ˜ e 2.25
U ˜ e 0.45
L ˜ 30
μ ˜ e 0.5
μ ˜ i 0.5
η ˜ 0.25
ρ ˜ 0 g ˜ 1 · 10 5
Table 2. Finite elements used for example 1.
Table 2. Finite elements used for example 1.
VariableFinite Element
u Q 2
D RT 1
U i DGQ 2
A Q 2
φ DGQ 1 ( Q 2 for scalar potential formulation)
p DGP 1
p i DGQ 2
Table 3. The parameters used for the hydrogel example.
Table 3. The parameters used for the hydrogel example.
QuantityValue
t ˜ 0 0
t ˜ e 5
c ¯ ˜ 0 ext 0.9
c ¯ ˜ e ext 0.1
H ˜ 1
c ˜ 2 H O 0 41.509
c ˜ 0 A + 0.75415
c ˜ 0 B 0.60415
p ˜ 0 0.011077
φ ˜ 0 0.11088
μ ˜ 0.040360
λ ˜ 0.016144
V ˜ 2 H O m 0.018069
n ˜ 0 0.25
z ˜ A + 1
z ˜ B 1
z ˜ 1
c ˜ 0.15
D ˜ 2 H O 14.673
D ˜ A + 0.05
D ˜ B 0.08
R ˜ 1
F ˜ 1
T ˜ 1
Table 4. Finite elements used for example 2.
Table 4. Finite elements used for example 2.
VariableFinite Element
u Q 2
c 2 H O , c i , p, φ DGP 1
I 2 H O N , I N i Q 2
η 2 H O , η i Q 2
φ s Q 2
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Stark, S. A Systematic Approach to Standard Dissipative Continua. Axioms 2023, 12, 267. https://doi.org/10.3390/axioms12030267

AMA Style

Stark S. A Systematic Approach to Standard Dissipative Continua. Axioms. 2023; 12(3):267. https://doi.org/10.3390/axioms12030267

Chicago/Turabian Style

Stark, Sebastian. 2023. "A Systematic Approach to Standard Dissipative Continua" Axioms 12, no. 3: 267. https://doi.org/10.3390/axioms12030267

APA Style

Stark, S. (2023). A Systematic Approach to Standard Dissipative Continua. Axioms, 12(3), 267. https://doi.org/10.3390/axioms12030267

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop