
Stephan

Teichtmeister

Variational Methods for
Dissipative Multifield Problems

in Solid Mechanics

8
Publication series of the
Institute of Applied Mechanics (IAM)


Variational Methods for Dissipative Multifield
Problems in Solid Mechanics

Von der Fakultät Bau- und Umweltingenieurwissenschaften der Universität
Stuttgart zur Erlangung der Würde eines Doktor-Ingenieurs (Dr.-Ing.)

genehmigte Abhandlung

vorgelegt von

Stephan Benjamin Teichtmeister

aus Graz

Hauptberichter:
Mitberichter:

Prof. Dr.-Ing. Marc-André Keip
Prof. Dr. rer. nat. Klaus Hackl

Assoc. Prof. Christian Linder, Ph.D.

Tag der mündlichen Prüfung: 11. September 2020

Institut für Mechanik (Bauwesen) der Universität Stuttgart

2021


Publication series of the Institute of Applied Mechanics (IAM), Volume 8
Institute of Applied Mechanics
University of Stuttgart, Germany, 2021

Editors:

Prof. Dr.-Ing. Dr. h.c. W. Ehlers
Dr.-Ing. Dipl.-Math. techn. F. Fritzen
Prof. Dr.-Ing. M.-A. Keip
Prof. Dr.-Ing. H. Steeb

Organisation und Verwaltung:

Institut für Mechanik (Bauwesen)
Lehrstuhl für Materialtheorie
Universität Stuttgart
Pfaffenwaldring 7
70569 Stuttgart
Tel.: +49 (0)711 685-66378
Fax: +49 (0)711 685-66347

c© Stephan Benjamin Teichtmeister
Institut für Mechanik (Bauwesen)
Lehrstuhl für Materialtheorie
Universität Stuttgart
Pfaffenwaldring 7
70569 Stuttgart
Tel.: +49 (0)711 685-66378
Fax: +49 (0)711 685-66347

Alle Rechte, insbesondere das der Übersetzung in fremde Sprachen, vorbehalten. Ohne
Genehmigung des Autors ist es nicht gestattet, dieses Heft ganz oder teilweise auf foto-
mechanischem Wege (Fotokopie, Mikrokopie) zu vervielfältigen.

ISBN 978-3-937399-56-0 (D 93 Stuttgart)


Abstract

The present work deals with the construction, analysis and numerical implementation

of constitutive variational principles for the description of dissipative material behavior

at small and large deformations. Such variational formulations are particularly suitable

for coupled multifield problems and allow to model the interactions of different physical

fields in a mathematically elegant way. In this context, special attention is paid to the

versatility and flexibility of variational principles. A broad class of dissipative processes in

solids are treated which are often due to structural irreversible changes on the microscale.

In this sense, there is a close link to the laws of thermodynamics and the principles of

material theory. After a short introduction to these, a general recipe for constructing

variational principles of simple rheological models is presented. This will already show

essential concepts of the approach. After that, the focus is put on crack propagation in

ductile materials. The phenomenon is treated within the framework of innovative phase

field approaches which are understood as specific gradient-damage theory, but with roots

in fracture mechanics. The coupling is done with a model of gradient-plasticity in order

to scale not only the damage but also the plastic zone. Furthermore, a micromorphic

extension is incorporated which, on the numerical side, makes the capture of spatial

elastic-plastic transitions easier. Next, the phase field method is also used to model crack

growth in anisotropic brittle materials. The basic idea is an anisotropic modification of

the purely geometrically motivated crack surface density based on the theory of tensor

invariants. The resulting structural tensors represent the microstructure of the material.

A modified length scale is obtained, which subsequently leads to an anisotropic Griffith

critical energy release rate. The latter causes crack propagation directions that differ

drastically from those in isotropic materials. This is shown in several numerical examples

that are consistent with experimental findings documented in the literature. Another

part of the thesis discusses a new minimization principle for the coupled problem of fluid

transport in fully saturated poroelastic media. It is based on the identification of the

fluid flow vector as canonical variable which, together with the deformation, determines

the energy and dissipation functionals. Here, particular attention is paid to the advan-

tages of the associated finite element formulation which is not restricted by the discrete

inf-sup stability condition. Formally, various saddle point formulations can be deduced

from this minimization principle, which underlines its canonical nature. The last part

of the thesis will present a general variational framework for gradient-extended continua,

which also takes into account the strong coupling with temperature evolution. A generic

internal variable is introduced which, together with its first gradient and the entropy,

enters the constitutive functions. The flexibility of the formulation is shown by means of

several modeling examples. In particular, nonisothermal processes are considered within

the Cahn-Hilliard diffusion, gradient-damage and gradient-plasticity. The regularizing

effect of such gradient theories is shown by means of plastic shear band evolution. In ad-

dition, a reduced minimization principle of thermo-gradient-plasticity is used to predict

the localization of elasto-plastic fields in a tensile bar.


Zusammenfassung

Die vorliegende Arbeit befasst sich mit der Konstruktion, Analyse und numerischen Im-

plementation von konstitutiven Variationsprinzipien zur Beschreibung dissipativen Mate-

rialverhaltens bei kleinen und großen Deformationen. Solche variationellen Formulierun-

gen eignen sich besonders für gekoppelte Mehrfeldprobleme und bilden die Interaktionen

verschiedener physikalischer Größen in einer mathematisch eleganten Art und Weise ab.

In diesem Rahmen wird ein besonderes Augenmerk auf die Vielseitigkeit und Einsatz-

flexibilität von Variationsprinzipien gelegt. Es wird eine breite Klasse von dissipativen

Vorgängen in Festkörpern behandelt, welche oft auf strukturelle irreversible Änderungen

auf der Mikroskala zurückzuführen sind. In diesem Sinne gibt es eine enge Verknüpfung

mit den Gesetzen der Thermodynamik und den Prinzipien der Materialtheorie. Nach einer

kurzen Einführung in diese, wird ein allgemeines Rezept zur Konstruktion von Variati-

onsprinzipien anhand einfacher rheologischer Modelle dargelegt. Hierdurch werden bereits

wesentliche Konzepte aufgezeigt. Danach wird das Augenmerk auf die Rissausbreitung

in duktilen Werkstoffen gelegt. Die Behandlung geschieht im Rahmen innovativer Pha-

senfeldansätze, welche als spezifische Gradienten-Schädigungstheorie verstanden werden

können, jedoch eng mit Konzepten der Bruchmechanik verbunden sind. Die Kopplung

erfolgt mit einem Modell der Gradienten-Plastizität um neben der Schädigungszone auch

die plastische Zone skalieren zu können. Des Weiteren wird eine mikromorphe Erweite-

rung eingearbeitet, welche auf der numerischen Seite vor allem das Erfassen der elastisch-

plastischen räumlichen Übergänge erleichtert. Ebenso im Rahmen der Phasenfeldmethode

wird als Nächstes ein Modell zur Beschreibung des Risswachstums in anisotropen spröden

Materialien vorgestellt. Die grundlegende Idee ist dabei eine anisotrope Modifikation der

rein geometrisch motivierten Rissflächendichte basierend auf der Theorie der Tensorinva-

rianten. Die sich dabei ergebenden Strukturtensoren bilden die Mikrostruktur des Materi-

als ab. Man erhält eine modifizierte Längenskala, welche in weiterer Folge eine anisotrope

kritische Griffith-Energiefreisetzungsrate bedeutet. Letztere bewirkt Rissausbreitungsrich-

tungen, die sich fundamental von denen in isotropen Materialien unterscheiden. Dies wird

in einigen numerischen Beispielen gezeigt, die im Einklang mit experimentellen Befunden

aus der Literatur sind. Ein weiterer Abschnitt diskutiert ein neues Minimierungsprinzip

für das gekoppelte Problem des Fluidtransports in vollgesättigten poroelastischen Medi-

en. Es basiert auf der Identifikation des Fluidflussvektors als kanonische Variable, welche

zusammen mit der Deformation das Energie- und Dissipationsfunktional bestimmt. Hier

wird insbesondere auf Vorteile der zugehörigen finite Elemente Formulierung eingegangen,

welche nicht durch die diskrete inf-sup Stabilitätsbedingung eingeschränkt ist. Rein formal

lassen sich des Weiteren verschiedene Sattelpunktformulierungen aus dem Minimierungs-

prinzip herleiten, was dessen kanonische Natur unterstreicht. Als Letztes wird eine allge-

meine Variationsstruktur für gradientenerweiterte Kontinua vorgestellt, welche die starke

Kopplung mit der Evolution der Temperatur berücksichtigt. Es wird eine generische in-

terne Variable eingeführt, welche zusammen mit ihrem ersten Gradienten und der Entro-

pie Eingang in konstitutive Funktionen findet. Die Flexibilität der Formulierung wird

anhand von Modellbeispielen gezeigt. Im Speziellen wird auf nicht-isotherme Prozesse in

der Cahn-Hilliard Diffusion, Gradienten-Schädigung und Gradienten-Plastizität eingegan-

gen. Der regularisierende Effekt solcher Gradiententheorien wird anhand der Formation

plastischer Scherbänder gezeigt. Zusätzlich wird ein reduziertes Minimierungsprinzip der

Thermo-Gradienten-Plastizität herangezogen, um die Lokalisation des elasto-plastischen

Feldes in einem Stab unter Zugbelastung vorherzusagen.


Danksagung

Die vorliegende Arbeit entstand während meiner Tätigkeit als wissenschaftlicher Mitar-
beiter am Institut für Mechanik (Bauwesen), Lehrstuhl für Materialtheorie der Universität
Stuttgart.

An erster Stelle möchte ich mich herzlich bei Prof. Dr.-Ing. Christian Miehe bedanken, der
leider im Jahr 2016 verstorben ist. Die gemeinsame Zeit mit ihm war sehr produktiv und
lehrreich und ich denke behaupten zu können, dass ich durch ihn die Kontinuumsmechanik
erst richtig begonnen habe zu verstehen. Er war nicht nur ein Förderer, sondern auch
jemand, mit dem man sich über alles Mögliche unterhalten konnte, wenn er nicht gerade
im Sinne des wissenschaftlichen Fortschritts einer sich selbst auferlegten Stresssituation
ausgesetzt war. Ich möchte die (leider nur kurze) prägende Zeit mit ihm nicht missen.

Ein besonderer Dank gilt Prof. Dr.-Ing. Marc-André Keip für die Unterstützung in schwie-
rigen Zeiten und Betreuung dieser Arbeit. Für das Interesse an der Arbeit, die Übernahme
des Koreferats und die anregenden Diskussionen möchte ich mich bei Prof. Dr. rer. nat.
Klaus Hackl und Assoc. Prof. Christian Linder, Ph.D. bedanken.

Ich habe mich über die Jahre in Stuttgart sehr wohl gefühlt, und das vor allem wegen des
sehr guten Arbeitsklimas, der Toleranz, des Zusammenhalts und der großen Hilfsbereit-
schaft innerhalb der Gruppe. Für die angenehme Büroatmosphäre, anregenden Diskussio-
nen und die produktive Zeit möchte ich Steffen Mauthe danken, ebenso Elten Polukhov.
Auch danken möchte ich Fadi Aldakheel, Daniel Vallicotti und Daniel Kienle für die sehr
gute Zusammenarbeit. Nicht weniger zum Gelingen dieser Arbeit haben Lukas Böger, Aref
Nateghi, Felix Göküzüm, Matthias Rambausek, Ashish Sridhar, Omkar Nadgir, Khiem
Nguyen und Nikolai Arnaudov beigetragen. Weiters möchte ich mich auch bei unseren bei-
den Sekretärinnen Nadine Steinecke und Leonie Fischer, sowie meinen beiden Kollegen
Lukas Eurich und Patrick Schröder vom Lehrstuhl für Kontinuumsmechanik bedanken.
Besonders dankbar bin ich auch für die sehr schöne und angenehme Zeit in der Wohnge-
meinschaft mit Aref Nateghi und meinen ehemaligen MitbewohnerInnen Daniel Malcher,
Emma Le Floch und Chiara Boccelatto. Eine solche harmonische und respektvolle Wohn-
gemeinschaft seh ich nicht als selbstverständlich an.

Ohne dieses persönliche Umfeld wäre diese Arbeit nicht möglich geworden. Es sind dabei
sehr gute Freundschaften entstanden und ich hoffe, dass diese trotz Entfernung aufrecht-
erhalten werden können.

Zuletzt möchte ich mich ganz herzlichen bei meiner Familie – meiner Mama Gabriele,
meinem Papa Günter und meiner Schwester Vanessa – bedanken, vor allem für ihre Ge-
duld, große Hilfsbereitschaft und Unterstützung. Ohne sie wäre ich nie so weit gekommen.
Auch den anderen Familienmitgliedern und Bekannten, besonders meiner Oma Maria und
meiner Patentante Brigitte, sei hier gedankt.

Stephan Teichtmeister


Contents i

Contents

I Introduction 1

1. Motivation and Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2. Continuum Mechanics and Material Theory . . . . . . . . . . . . . . . . . . . . . 7

2.1. Geometry of Finite Deformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.1. Fundamental Mappings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.2. Stretch and Convected Metric Tensors . . . . . . . . . . . . . . . . . . . . . 9

2.1.3. Velocity Gradient, Rate-of-Deformation and Spin Tensors . . . . . . . 10

2.2. Stress Tensors and Heat Flux Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.2.1. Stress Power Expressions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2.2. Heat Flux Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3. Physical Balance Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3.1. Balance of Mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3.2. Balance of Linear Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.3. Balance of Angular Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.4. Balance of Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.3.5. Balance of Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.4. Concepts of Material Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.1. Principle of Equipresence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.2. Principle of Determinism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.3. Concept of Internal Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.4. Principle of Local Action . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.5. Principle of Objectivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4.6. Second Law of Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.4.7. Principle of Material Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3. Variational Principles for Rheological Models . . . . . . . . . . . . . . . . . . . . 23

3.1. Elementary Formulation of Basic Rheological Devices . . . . . . . . . . . . . . . 23

3.1.1. Variational Principles for Elastic Device . . . . . . . . . . . . . . . . . . . . 23

3.1.2. Variational Principles for Viscous Device . . . . . . . . . . . . . . . . . . . 26

3.1.3. Variational Principles for Friction Device . . . . . . . . . . . . . . . . . . . 29

3.2. General Setting for Rheological Models . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.3. Modeling Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35


ii Contents

3.3.1. Visco-Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.3.2. Perfect Elasto-Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.3. Elasto-Plasticity with Kinematic Hardening . . . . . . . . . . . . . . . . . 39

3.3.4. Visco-Elasto-Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4. Variational Phase Field Approach to Brittle Fracture . . . . . . . . . . . . . 45

4.1. The Fundamental Variational Theory of Brittle Fracture . . . . . . . . . . . . . 45

4.1.1. Fundamental Variational Statement of Brittle Fracture . . . . . . . . . 45

4.1.2. Regularized Variational Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2. Variational Gradient Damage Approaches to Brittle Fracture . . . . . . . . . . 47

4.2.1. Formulation Without Threshold . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.2.2. Formulation With Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

5. List and Short Summaries of Appended Papers . . . . . . . . . . . . . . . . . . 49

5.1. Paper A: Phase-field modelling of ductile fracture: a variational gradient-
extended plasticity-damage theory and its micromorphic regularization . . . 49

5.2. Paper B: Phase field modeling of fracture in anisotropic brittle solids . . . 49

5.3. Paper C: Aspects of finite element formulations for the coupled problem
of poroelasticity based on a canonical minimization principle . . . . . . . . . 49

5.4. Paper D: A Variational Framework for the Thermomechanics of Gradient-
Extended Dissipative Solids. Formulation, Model Problems and Stability
Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

II Appended Papers 55


Part I

Introduction


3

1. Motivation and Overview

The goal of this work is to present rate-type and incremental variational formulations
for smooth as well as nonsmooth dissipative processes in solids that undergo small and
large deformations. Dissipative problems usually have a multifield character which means
that beside the mechanical deformation, also other physical fields arise as primary vari-
ables. Examples for these other physical fields are the temperature, plastic or damage
variables, species concentrations or fluxes and electric polarization or magnetization fields.
The fundamental question is, how different physical quantities interact with each other
or, in other words, which physical coupling effects occur and how they can be modelled
accordingly within a thermodynamically consistent framework. This is where rate-type
or incremental variational principles take a special role. Roughly speaking, these formu-
lations are based on two constitutive functions only, namely an energy function and a
dissipation potential function. These functions are formulated in terms of a set of canon-
ical variables and their rates, respectively, representing the current multifield state of a
material point. Then, the optimization of a rate-type or incremental global variational
potential governs the full dissipative response of the material at current time. At this
point we want to mention the works Simo & Honein [46], Ortiz & Stainier [42],
Carstensen et al. [9] and Miehe [27] among many others. Variational principles
have a strong versatile character and offer a great flexibility. In this sense, special focus
should be put on canonical formulations which represent the most fundamental problem
statements. If such principles are at hand, other physical variables of interest such as
dissipative driving forces can be included via Legendre transformations. Additionally, if
the incremental variational principle has a pure minimization structure and a nonconvex
potential is at hand, an analysis of material and structural instabilities is possible by
elaborating weak convexity conditions stemming from direct methods in the calculus of
variations, see Dacorogna [10]. With these methods, the formation of microstructures
observed in experiments or buckling phenomena can be predicted. Another important
feature of variational principles is their inherent symmetry which directly influences the
numerical solution by the finite element method. Especially, within a typical Newton-
Rapshon iteration step, the finite element stiffness matrix is governed by the Hessian
of the underlying potential density and is therefore symmetric. On the other hand, the
residuum vector is determined by the Jacobian of this potential density.

The thesis consists of four scientific papers. Each article investigates one of the
following coupled dissipative processes by multifield variational principles and includes
the numerical implementation by the finite element method: (A) evolution of fracture in
solids undergoing small elastic and large plastic deformations, (B) evolution of fracture in
anisotropic brittle materials, (C) diffusion in largely deforming poroelastic media and (D)
evolution of macro- and microstructural fields strongly coupled to the evolving thermal
state in gradient-extended solids undergoing small and large deformations.

The formulations elaborated for (A)–(B) fall into the category of recently developed
phase field models for fracture. These models regularize sharp crack discontinuities within
a pure continuum approach. The fundamental basis for these models is the seminal work
of Francfort & Marigo [14] who propose an energetic minimization principle for the
Griffith-type theory of brittle fracture in isotropic solids. Their formulation overcomes
drawbacks of the classical Griffith theory of linear fracture mechanics, such as the ne-
cessities of an initial crack and pre-knowledge of the crack path. The regularized setting


4 1 Motivation and Overview

of Francfort & Marigo’s variational formulation for fracture is considered in Bourdin

et al. [7] and is inspired by the work of Ambrosio & Tortorelli [3]. The regulariza-
tion comes along with the introduction of an auxiliary variable which can be considered
as a phase field that interpolates between an unbroken and a fully broken state. An im-
portant feature of this model is the treatment of the irreversibility in a time incremental
setting by evolving internal Dirichlet-type conditions on the phase field, see also Kuhn

& Müller [21]. Later, Miehe et al. [33, 32] formulated a phase field model for brit-
tle fracture based on rigorous thermodynamic arguments and the principle of maximum
dissipation. In big contrast, this model uses a local inequality constraint on the rate of
the fracture phase field and can hence be seen as a gradient-damage model, however with
roots in fracture mechanics, see also Amor et al. [4] and Pham et al. [43]. In this
thesis, extensions of the approach of Miehe towards (A) fracturing in ductile materials
and (B) fracture propagation in anisotropic brittle materials are considered. Former are
based on a specific gradient damage formulation which is strongly coupled with models
of elasto-plasticity, see also Ambati et al. [2], Alessi et al. [1] and Miehe et al.

[34]. Key aspects are the coupling with a gradient-plasticity model offering a scaling of
the width of damage and plastic zones, the formulation of a work density function that
splits into energetic and dissipative parts and a micromorphic extension in the sense of
Forest [13]. The modeling of fracturing in anisotropic solids is based on a modification
of the crack surface density function by making use of the theory of tensor invariants.
This leads to definitions of characteristic structural tensors of second and fourth order.
As an overall consequence, an anisotropic Griffith’s critical energy release rate is obtained.

The incremental minimization formulation associated with (C) is the most funda-
mental statement for the coupled problem of Darcy-Biot-type fluid transport in porous
media. As suggested in Miehe et al. [35], the key for attaining a minimization struc-
ture is the definition of rate-of-energy and dissipation potential functionals in terms of
the fluid flux. In this thesis, a special focus is put on conforming and nonconforming
minimization-based finite element designs. They show specific advantages over finite el-
ement formulations based on the classical saddle point principle of poroelasticity which
determines the current deformation and fluid pressure. To be more specific, finite element
formulations of saddle point principles are constrained by the discrete inf-sup condition
which, especially for the classical saddle point principle of poroelasticity, is difficult to
handle with while keeping computational efficiency in mind, see Brezzi & Fortin [8].
A violation of the discrete inf-sup condition often results in unphysical oscillations in
primary fields that propagate inside the domain. Stabilized methods are able to remove
spurious modes but depend on a stabilization parameter that is in general not easy to
identify, see Preisig & Prévost [44]. In contrast, minimization-based finite element
formulations are not affected by the discrete inf-sup condition. However, as a disadvan-
tage one has to point out possible locking behavior in the incompressible fluid limit. As
shown, this can be overcome by a selective reduced integration method which is connected
to an extended saddle point principle that relaxes overconstraints.

Finally, for (D) we develop a versatile incremental variational formulation that is ap-
plicable to a wide range of model problems. Especially, this part of the thesis presents an
extension of the general framework of gradient-extended dissipative solids at large strains
by Miehe [30] to nonisothermal processes. Starting point is the identification of the en-
tropy and the entropy rate, respectively, as fundamental thermal variables which enter the
internal energy and canonical dissipation potential functions. Rigorously distinguishing


5

between thermodynamic quantities dual to the entropy and entropy rate, a (generalized)
Legendre transformation of the canonical dissipation potential function yields as opti-
mality condition an entropy evolution equation that must render the structure of the
energy equation. This results in a condition for the “thermally dual” dissipation poten-
tial function which is in line with the fundamental findings of Yang et al. [49] who
deal with local thermo-visco-plasticity. In addition, extended variational principles that
include dissipative driving forces as well as threshold mechanisms are constructed in this
part of the thesis. The versatility of the formulation is underlined by modeling examples
such as Cahn-Hilliard diffusion, gradient-damage and gradient-plasticity strongly coupled
with temperature evolution. For latter model problem, the formation of a shear band is
numerically analysed. Finally, for a specific gradient-plasticity model as well as its micro-
morphic extension, an analytical stability analysis of homogeneous elasto-plastic solutions
in a tensile bar undergoing thermal softening is performed.

The introductory part of the thesis is organized as follows: In Chapter 2, the basics of
continuum mechanics and thermodynamics as well as material theory for solid materials
undergoing large deformations are summarized. A strong accent is put on geometrically
exact formulations which are indispensable when, for example, formulating large strain
plasticity theories. In Chapter 3, the basic recipe for setting up rate-type variational
formulations is shown by means of several rheological models that are combinations of
three basic devices. These three devices are the elastic spring, the viscous dashpot and
the friction device. Even though the models considered are quite simple, the essence
of the variational approach, which is the key subject of the thesis, can be shown quite
nicely. The basics of the variational phase field approach to fracture in brittle solids is
shortly presented in Chapter 4. Here, the fundamental differences of existing models
are discussed and illustrated by means of a tension test of a 1D bar. In Chapter 5, the
appended papers are listed and shortly summarized.


7

2. Continuum Mechanics and Material Theory

In this chapter, we shortly review and summarize some important aspects of contin-
uum mechanics and material theory. To describe large deformations of solid bodies and
their states of stress, we use concepts from tensor calculus on manifolds. In this sense,
a main aspect is the elaboration of the tensors’ mapping properties. The basis for this
chapter are lecture notes of Miehe [31] who taught a course on geometrical methods of
continuum mechanics. There exists a wide literature on the mechanics and thermodynam-
ics of solids and we refer to the books Truesdell & Noll [48], Malvern [24], Ogden

[41], Marsden & Hughes [25], Holzapfel [19], Haupt [18] and Gurtin et al. [16],
just to name a few. In addition, we want to point out the book Krawietz [20] which
is written in German. Especially in Marsden & Hughes [25] concepts of differential
geometry on manifolds are used – for this topic we for example refer to Lee [22].

2.1. Geometry of Finite Deformations

In this section, we give a summary on kinematics and introduce the four fundamental
mappings, namely the point, tangent, area and volume maps. In addition, we define
convected metric tensors and deformation rates.

2.1.1. Fundamental Mappings. Let B ⊂ R
ω with dimension ω ∈ {2, 3} be the

reference placement of a material body B into the Euclidean space. The motion of the
body within a time interval T = (0, te) ⊂ R+ is described by the vector-valued deformation
map

ϕ :

{
B × T → R

ω

(X, t) 7→ x = ϕ(X, t) ,
(2.1)

which maps at time t ∈ T points X ∈ B from the reference configuration B to points
x ∈ S of the current configuration S = ϕt(B). The deformation map is restricted by
the pointwise condition (2.4) which makes it one-to-one. We consider the reference as
well as the current configuration as Riemannian C∞-manifolds with coordinate systems
{XA}A=1,ω and {xa}a=1,ω. At every point X ∈ B, the associated differential operators
EA = ∂

∂XA with A = 1, ω span a vector space called tangent space TXB. The cor-
responding dual space at X ∈ B is called cotangent space T ∗

XB and is formed by the
differential operators EA = dXA with A = 1, ω. Latter are defined by the duality prod-
uct 〈EA,EB〉 = δAB in terms of the Kronecker symbol. Likewise, at every point x ∈ S the
associated differential operators ea = ∂

∂xa and ea = dxa with a = 1, ω span the tangent
and cotangent spaces TxS and T ∗

xS, respectively. We denote the Euclidean metric ten-
sors associated with the reference and current configurations by G = GABE

A ⊗EB and
g = gabe

a ⊗ eb. Throughout this text, we assume for simplicity Cartesian coordinates on
both manifolds, i.e. the components of the metric tensors are represented by Kronecker
symbols GAB = δAB and gab = δab.

Next, we define the deformation gradient as the tangent to the deformation map

F (X , t) = Dϕ(X, t) or F a
A =

∂ϕa

∂XA
. (2.2)

It is a mixed-variant two-point tensor and maps material tangent vectors T ∈ TXB to


8 2 Continuum Mechanics and Material Theory

spatial tangent vectors t ∈ TxS, i.e.

F :

{
TXB → TxS
T 7→ t = F T ,

(2.3)

see Figure 1. In addition, we have the Jacobian

J(X, t) = detF (X, t) > 0 (2.4)

which has to be positive everywhere in order to rule out penetration of matter. Note
carefully, that the term gradient for F is geometrically unprecise. To see this, consider
a transformation of coordinates {XA} → {X̄A} and {xa} → {x̄a} on the reference and
current configurations. We obtain by applying the chain rule

F a
A =

∂ϕa

∂XA
=

∂

∂XA
xa(ϕ̄b(X̄B)) =

∂xa

∂x̄b
∂ϕ̄b

∂X̄B

∂X̄B

∂XA
(2.5)

and identify the components of the deformation gradient in the (possibly curvelinear)
coordinate systems {X̄A} and {x̄a} as F̄ a

A = ∂ϕ̄a/∂X̄A. Hence, in contrast to components
of real gradients in curvelinear coordinate systems, the components of F do not contain
Christoffel symbols of second kind.

As another fundamental mapping, we consider the one which maps material area
elements A ∈ T ∗

XB to spatial area elements a ∈ T ∗
xS. Considering distinct material

tangent vectors T 1, T 2 ∈ TXB, we define A = T 1×T 2 = NA in terms of a material unit
normal vector N ∈ T ∗

XB and likewise a = FT 1 × FT 2 = n a in terms of a spatial unit
normal vector n ∈ T ∗

xS. Note that normal vectors are usually considered as one-forms
and as such live in the cotangent spaces. The relationship between A and a is given by
the cofactor of the deformation gradient

cof F :

{
T ∗
XB → T ∗

xS
A 7→ a = JF−T A

(2.6)

which is known as Nanson’s formula. Accordingly, the linear map F−T : T ∗
XB → T ∗

xS is
denoted as normal map.

Finally, we construct a material volume element V = 〈T 1,T 2×T 3〉 via three disctinct
tangent vectors T 1, T 2, T 3 ∈ TXB and a spatial volume element v = 〈FT 1,FT 2×FT 3〉.
The relationship between both is given by the Jacobian

J :

{
R+ → R+

V 7→ v = J V
(2.7)

and is referred to as volume map. At this point, we want to point out that the determinant
of a second order tensor equals the determinant of its component matrix in any coordinate
system only if the second order tensor is an endomorphism. The deformation gradient F
is a two-point tensor, and hence we have detF 6= det[F a

A] in general. However, since the
underlying geometry is Euclidean, we have TXB ∼= R

ω and TxS ∼= R
ω and we can see F

as endomorphism if we choose coordinate systems such that the metric coefficients are the
same for every pairing (X,ϕ(X, t)), e.g. Cartesian coordinates. Hence, detF = det[F a

A]
if the coordinate systems {XA} and {xa} are Cartesian, but

detF = det[F̄ a
A]

√
det[ḡab]√
det[ḠAB]

(2.8)

for general (possibly curvelinear) coordinate systems {X̄A} and {x̄a}, see Marsden &

Hughes [25].


2.1 Geometry of Finite Deformations 9

replacemen

X x

B S

F = Dϕ(X, t)

ϕ(X, t)
R

ω

T t

Figure 1: Geometry of finite deformations. The deformation map ϕ(X , t) maps at time
t points X ∈ B to points x ∈ S. The tangent F = Dϕ to this map is called deformation
gradient and maps tangent vectors T ∈ TXB to tangent vectors t ∈ TxS.

2.1.2. Stretch and Convected Metric Tensors. Let a material reference direction
M ∈ TXB with ‖M‖G = 1 be given. The stretch vector λM measures the change of the
deformation in direction M via the directional derivative. We get

(λM )a =
d

dǫ

∣∣∣∣
ǫ=0

ϕa(XA + ǫMA) = F a
AM

A or λM = FM . (2.9)

Due to the mapping properties of the deformation gradient, we see that λM ∈ TxS is
an object that lives in the tangent space associated with a point x ∈ S of the current
configuration. The stretch λM associated with the direction M is defined as the norm of
the stretch vector

λM = ‖λM‖g =
√

〈λM , gλM〉 =
√

(λM )a(λM )bδab . (2.10)

Inserting (2.9) yields when using the definition of the transposition of a second order
tensor

λM =
√
〈FM , gFM〉 =

√
〈M ,CM〉 =

√
MAMBCAB = ‖M‖C (2.11)

in terms of the so called right Cauchy-Green tensor

C = F TgF or CAB = F a
AF

b
B δab . (2.12)

This tensor is symmetric and positive definite and relation (2.11) shows that C : TXB →
T ∗
XB plays the role of a metric on the reference configuration. Especially, we say that

the right Cauchy-Green tensor is the convected current metric obtained by a pull-back
operation

C = ϕ∗
t (g) = F TgF , (2.13)

see the commutative diagram in Figure 2a).

Alternative to the spatial definition (2.9) of the stretch, we consider a material one and
start with considering the spatial reference direction m ∈ TxS with ‖m‖g = 1. Based on

λM = FM = λMm (2.14)

we obtain 1/λM M = F−1m and introduce the material stretch vector

Λm = F−1m or (Λm)A = (F−1)Aam
a (2.15)


10 2 Continuum Mechanics and Material Theory

TXBTXB TxSTxS

T ∗
XBT ∗

XB T ∗
xST ∗

xS

C g

FF

F−TF−T

G c

a) b)

Figure 2: Convected metric tensors. a) The right Cauchy-Green tensor C = F T gF is the
pull-back of the spatial metric g. b) The Finger tensor c = F−TGF−1 is the push-forward
of the material metric G.

associated with the directionm. Due to the mapping properties of the inverse deformation
gradient, we see that Λm ∈ TXB is an object that lives in the tangent space associated
with a point X ∈ B of the reference configuration. The inverse stretch can now be
expressed by the norm of the material stretch vector

1

λM
= ‖Λm‖G =

√
〈Λm,GΛm〉 =

√
(Λm)A(Λm)BδAB . (2.16)

Inserting (2.15) gives, when using the definition of the transposition of a second order
tensor, the alternative representation

1

λM
=

√
〈F−1m,GF−1m〉 =

√
〈m, cm〉 =

√
mambcab = ‖m‖c (2.17)

in terms of the so called Finger tensor

c = F−TGF−1 or cab = (F−1)Aa(F
−1)Bb δAB . (2.18)

This tensor is again symmetric and positive definite and relation (2.17) shows that
c : TxS → T ∗

xS plays the role of a metric on the current configuration. We say that
the Finger tensor is the convected reference metric obtained by a push-forward operation

c = ϕt∗(G) = F−TGF−1 , (2.19)

see the commutative diagram in Figure 2b).

2.1.3. Velocity Gradient, Rate-of-Deformation and Spin Tensors. We define
the spatial velocity field

v(x, t) = (
∂

∂t
ϕ(X, t)) ◦ ϕ−1(x, t) (2.20)

which gives for every point x ∈ S a velocity vector v ∈ TxS. The spatial velocity gradient
field is then defined as1

l(x, t) = ∇v(x, t) = (
∂

∂t
F (X, t) ◦ ϕ−1(x, t))F−1(x, t) (2.21)

1We denote by ˙(•) = ∂
∂t
(•)(X , t) the material time derivative.


2.2 Stress Tensors and Heat Flux Vectors 11

and gives for every point x ∈ S a velocity gradient

l = Ḟ F−1 or lab = Ḟ a
A(F

−1)Ab . (2.22)

It is a mixed-variant second order tensor with the mapping property l : TxS → TxS.
Especially, the spatial velocity gradient determines the rate of the stretch vector via

λ̇M = l λM or (λ̇M)a = lab(λM)b . (2.23)

To define symmetric as well skewsymmetric parts of the spatial velocity gradient, we lower
the first index by the spatial metric tensor and obtain

gl = l♭ = d+w with d = sym(gl) and w = skew(gl) (2.24)

or in index representation

δacl
c
b = dab + wab with dab =

1

2
[ lab + (l♭T )ab ] and wab =

1

2
[ lab − (l♭T )ab ] . (2.25)

The covariant second order tensors d : TxS → T ∗
xS and w : TxS → T ∗

xS are called rate-of-
deformation and spin tensors, respectively. Former plays an important role in the spatial
formulation of incremental material laws since it is related to the Lie derivative of the
spatial metric tensor

d = 1
2
£vg (2.26)

with £vg = ġ + l♭ + l♭T and ġ = 0 . Hence, in contrast to the spatial velocity gradient l
defined in (2.22), the rate-of-deformation tensor represents an objective deformation rate.
The corresponding commutative diagram is depicted in Figure 3a). To gain some physical
insight, consider a spatial unit vector m ∈ TxS and let the rate-of-deformation tensor d
be given. Then, the rate of the logarithmic stretch associated with the spatial direction
m is given by

˙lnλM = 〈m,dm〉 . (2.27)

To interpret the spin tensor, consider a spatial unit direction m ∈ TxS which now should
be an eigenvector of the mixed-variant rate-of-deformation tensor

g−1dm = ǫmm or δabdbcm
c = ǫmm

a , (2.28)

where the eigenvalue ǫm represents the rate of the logarithmic stretch in the associated
direction m according to (2.27). Then, the mixed-variant spin tensor governs the rate of
this eigenvector

ṁ = g−1wm or ṁa = δabwbcm
c , (2.29)

where clearly 〈ṁ, gm〉 = 0.

2.2. Stress Tensors and Heat Flux Vectors

In this section, we introduce the notion of stress. Especially, different stress tensors
are defined via the concept of conjugate variables and are geometrically related to each
other. In a way similiar to the concept of stress, heat flux vectors are defined on the
current and reference configurations.


12 2 Continuum Mechanics and Material Theory

2.2.1. Stress Power Expressions. One starts with assuming the existence of a trac-
tion vector field t(x, t;n) which depends on a spatial unit vector n. The field t(x, t;n)
represents the reaction force per unit area of a cut surface through the deformed contin-
uum. The orientation of the area element is indicated by the outer normal n. Then, the
Cauchy’s second order stress tensor relates at any point x ∈ S a given unit normal vector
n ∈ T ∗

xS to the traction vector t ∈ TxS, i.e. the traction vector field is

t(x, t;n) = σ(x, t)n or ta = σabnb . (2.30)

Hence, the Cauchy stress tensor is considered as a map σ : T ∗
xS → TxS. It should be

noted, that relation (2.30) relies on the balance of linear momenum. We define at time t
the stress power in the overall body as

∫

S

σ : d dv or

∫

S

σabdab dv (2.31)

in terms of the rate-of-deformation tensor. Using the volume map (2.7), we can alterna-
tively express the stress power in the overall body by

∫

S

σ : d dv =

∫

B

τ : d dV with τ = Jσ . (2.32)

The second order tensor τ is called Kirchhoff stress tensor and differs from the Cauchy
stress tensor just by the Jacobian. The expression

P = τ : d = τ : 1
2
£vg or P = τabdab (2.33)

represents the stress power per unit undeformed volume. We say that the metric tensor
g is the quantity that is energetically conjugate to the Kirchhoff stress tensor τ . Like the
Cauchy stress tensor, τ is a contravariant object defined in the current configuration. Its
pull-back defines the second Piola-Kirchhoff stress tensor

S = ϕ∗
t (τ ) = F−1τF−T or SAB = (F−1)Aa(F

−1)Bb τ
ab (2.34)

which is a map S : T ∗
XB → TXB, see Figure 3b). The stress power (2.33) per unit

undeformed volume can be rewritten as

P = S : 1
2
Ċ or P = SAB(1

2
Ċ)AB (2.35)

and we identify the right Cauchy Green tensor C as the quantity that is energetically
conjugate to the second Piola-Kirchhoff stress tensor S. Note that (2.35) relies on the
symmetry of the Cauchy stress tensor which is an outcome of the balance of angular
momentum, see Section 2.3.3. A partial pull-back of the Kirchhoff stress tensor τ , i.e.
a transformation with respect to the second slot of τ (·, ·), defines the contravariant first
Piola-Kirchhoff stress tensor

P = τF−T or P aA = τab(F−1)Ab (2.36)

which is a map P : T ∗
XB → TxS, see Figure 3b). It takes at any point X ∈ B as input a

unit normal vector N ∈ T ∗
XB, indicating an undeformed area element’s orientation, and


2.2 Stress Tensors and Heat Flux Vectors 13

TXBTXB TxSTxS

T ∗
XBT ∗

XB T ∗
xST ∗

xS

1
2
Ċ d

FF

F−TF−T

S τ
P

a) b)

Figure 3: Rates and stress tensors. a) The Lie derivative £vg of the spatial metric g is the
push-forward of the material time derivative Ċ of the convected metric C. b) The Kirchhoff
stress tensor τ is the push-forward of the second Piola-Kirchhoff stress tensor S. The first
Piola-Kirchhoff stress tensor P is obtained by a partial pull-back of the Kirchoff stress.

gives as result the nominal traction vector T ∈ TxS which represents a reaction force per
unit undeformed area, i.e. the nominal traction vector field is

T (X, t;N) = P (X, t)N or T a = P aANA . (2.37)

Like the deformation gradient, the first Piola-Kirchhoff stress is a two-point tensor. Hence,
it is not surprising that the deformation gradient is the quantity that is energetically
conjugate to the mixed-variant first Piola-Kirchhoff stress

P = gP : Ḟ or P = δabP
bAḞ a

A . (2.38)

In total, we have three equivalent representations of the stress power per unit undeformed
volume

P = gP : Ḟ = S : 1
2
Ċ = τ : 1

2
£vg (2.39)

which are referred to as two-point, material and spatial formulations.

2.2.2. Heat Flux Vectors. In formal analogy to before, we assume the existence
of a scalar heat flux field h(x, t;n) that depends on a spatial unit vector n. The field
h(x, t;n) represents the amount of heat which per unit time flows through unit area of a
cut surface through the deformed continuum body. The orientation of the area element is
indicated by the outer normal n. Then, the Cauchy or spatial heat flux vector relates at
any point x ∈ S a given unit normal vector n ∈ T ∗

xS to the heat flux h, i.e. the spatial
heat flux field is

h(x, t;n) = 〈−q(x, t),n〉 or h = −qana . (2.40)

The spatial heat flux vector field gives for every point x ∈ S a Cauchy heat flux vector
q ∈ TxS. The total amount of heat which per unit time enters/leaves a subregion U ⊂ S
of the deformed body is

∫

U

〈−q,n〉 da =

∫

ϕ−1
t

(U)

〈−JF−1q,N〉 dA (2.41)


14 2 Continuum Mechanics and Material Theory

where we have used Nanson’s formula (2.6). The second integrand in (2.41) represents
the nominal heat flux, that is the amount of heat transferred per unit time and per unit
undeformed area

H(X, t;N) = 〈−J(F−1 ◦ϕ)(q ◦ϕ),N〉 . (2.42)

This motivates the definition of the material heat flux vector

Q = JF−1q or QA = J(F−1)Aa q
a (2.43)

and we see that Q ∈ TXB. To underline the analogy to the different geometric settings
of the stress, Q is often called the Piola-Kirchhoff heat flux vector. It can formally be
defined by doing a pull-back operation on the Kirchhoff heat flux vector Jq that differs
from the Cauchy heat flux vector just by the Jacobian.

2.3. Physical Balance Laws

We consider an arbitrary subregion U ⊂ S of the current configuration and replace
the action of the surrounding part B \ U on ∂U by flux terms, i.e. tractions, heat and
entropy fluxes. The change of extensive physical quantities such as mass, linear and
angular momenta, energy and entropy are not only caused by fluxes, but also by source
and supply terms. These are body forces, heat sources and the entropy production.

2.3.1. Balance of Mass. The total mass associated with the subbody U of the
deformed continuum is defined as

m(U) =

∫

U

ρ dv (2.44)

in terms of the spatial mass density ρ which represents mass per unit deformed volume.
In its material form, the balance of mass states thatm(U) must be identical with the mass
associated with the corresponding subbody of the reference (undeformed) configuration

m(U) = m0(ϕ
−1
t (U)) with m0(ϕ

−1
t (U)) =

∫

ϕ−1
t

(U)

ρ0 dV (2.45)

in terms of the time independent material mass density ρ0 which represents mass per unit
undeformed volume. We immediately conclude from (2.45) the local statement of the
balance of mass

Jρ = ρ0 . (2.46)

If we do not want to refer to the reference configuration, we can write the balance of mass
in the alternative form

d

dt
m(U) = 0 . (2.47)

Then, Reynold’s transport theorem, which relies on the relationship J̇ = J div v, yields
after localization the balance of mass in its spatial form

ρ̇+ ρ div v = 0 . (2.48)

Here, ρ̇ denotes the material time derivative of the spatial mass density. Note that these
statements are the simplest forms of the balance of mass since no source and flux terms
are assumed.


2.3 Physical Balance Laws 15

2.3.2. Balance of Linear Momentum. We define the linear momentum associated
with the subbody U of the deformed continuum and the resultant force loading this
subbody

I(U) =

∫

U

ρv dv and F r(U) =

∫

U

ρb dv +

∫

∂U

t da . (2.49)

Here, b is a body force per unit mass and t = σn the traction. The balance of linear
momentum valid for an inertial frame states

d

dt
I(U) = F r(U) . (2.50)

Then, Gauss theorem applied to the boundary integral and Reynold’s transport theorem
together with the spatial form (2.48) of the mass balance yield after localization the spatial
form of the balance of linear momentum2

ρv̇ = ρb+ divσ or ρv̇a = ρba + σab
,b . (2.51)

Here, v̇ denotes the material time derivative of the spatial velocity. To obtain a material
local form of the balance of linear momentum, we multiply (2.51) by the Jacobian J ,
reparametrize the velocity and body force fields V = v ◦ ϕ and B = b ◦ ϕ, make use of
Piola’s identity3

J divσ = DivP or Jσab
,b = P aA

,A , (2.52)

see Marsden & Hughes [25], and insert the material form (2.46) of the balance of mass.
We obtain the result

ρ0V̇ = ρ0B +DivP or ρ0V̇
a = ρ0B

a + P aA
,A . (2.53)

We see that Cauchy’s stress σ and the first Piola-Kirchhoff stress P have fundamental
meanings with respect to the balance of linear momentum.

2.3.3. Balance of Angular Momentum. With respect to a point O fixed in an
inertial frame, we define the angular momentum associated with the subbody U and the
resulting couple force loading this subbody

DO(U) =

∫

U

x× ρv dv and MO(U) =

∫

U

x× ρb dv +

∫

∂U

x× t da . (2.54)

Here, x denotes the position vector pointing from the reference point O to the continuum
point x ∈ S. The concept of a position vector is ultimately related to Euclidean geometry.
As before, b denotes a body force per unit mass and t = σn the traction. Then, the
balance of angular momentum states

d

dt
DO(U) = MO(U) . (2.55)

2In an arbitrary (possibly curvelinear) coordinate system, the components of the divergence of the
Cauchy stress tensor field are expressed by the covariant derivative (divσ)a = σab

|b. For the special
case of a Cartesian coordinate system, the Christoffel symbols of second kind vanish and we simply have
(divσ)a = ∂σab/∂xb = σab

,b.
3Likewise, the components of the divergence of the first Piola-Kirchhoff stress tensor field in a general

(possibly curvelinear) coordinate system are governed by the covariant derivate (DivP )a = P aA
|A.


16 2 Continuum Mechanics and Material Theory

Its local statement gives the symmetry of the Cauchy stress tensor

σ = σT or σab = σba . (2.56)

Using the geometric relationships for the different stress tensors given in Section 2.2.1, we
obtain for the first and second Piola-Kirchhoff stress tensors the symmetry relations

FP T = PF T and S = ST or F a
AP

bA = P aAF b
A and SAB = SBA . (2.57)

The unsymmetry of the first Piola-Kirchhoff stress tensor is not a surprise. The relation
“P = P T” does not make sense since the tensors P and P T have different mapping
properties.

2.3.4. Balance of Energy. First, we consider the total energy associated with the
subbody U of the deformed continuum. It is the sum of the kinetic and internal energies
defined as

K(U) =

∫

U

1
2
〈ρv, gv〉 dv and E(U) =

∫

U

ρe dv (2.58)

where e is the specific internal energy which represents the internal energy per unit mass.
Second, the total power loading the subbody is the sum of mechanical and thermal con-
tributions

Pm(U) =

∫

U

〈ρb, gv〉 dv +

∫

∂U

〈t, gv〉 da and Pq(U) =

∫

U

ρr dv +

∫

∂U

h da . (2.59)

Here, r models a heat source with dimension power per unit mass and h = 〈−q,n〉 is the
heat flux crossing the boundary of the subbody. Then, the balance of energy or first law
of thermodynamics states

d

dt
[K(U) + E(U) ] = Pm(U) + Pq(U) . (2.60)

Applying Gauss theorem to the boundary integrals in (2.59) and using Reynold’s transport
theorem together with ġ = 0 and the spatial forms (2.48) and (2.51) of the mass and
linear momentum balances yield after localization the spatial form of the energy balance

ρė = σ : gl − div q + ρr or ρė = σabδac l
c
b − qa,a + ρr . (2.61)

Here, ė is the material time derivative of the specific internal energy. Note that σ : gl =
σ : d due to the symmetry of the Cauchy stress tensor. To obtain a material local form
of the balance of energy, we multiply (2.61) by the Jacobian J , reparametrize the specific
internal energy and heat source power fields e0 = e ◦ ϕ and r0 = r ◦ ϕ, recall the stress
power equivalence Jσ : d = gP : Ḟ , make use of Piola’s identity

J div q = DivQ or J qa,a = QA
,A (2.62)

and insert the material form (2.46) of the balance of mass. The result is

ρ0ė0 = gP : Ḟ −DivQ+ ρ0r0 or ρ0ė0 = δabP
bAḞ a

A −QA
,A + ρ0r0 . (2.63)

Note that due to ρ̇0 = 0 we can rewrite this equation such that the material density does
not show up explicitely

˙̃e0 = gP : Ḟ − DivQ+ r̃0 , (2.64)

where ẽ0 and r̃0 denote the internal energy and heat source power per unit undeformed
volume.


2.3 Physical Balance Laws 17

2.3.5. Balance of Entropy. We consider the total entropy associated with the
subbody U of the deformed continuum and the entropy flow as the sum of an entropy flux
across the boundary ∂U , an entropy supply and an entropy source term

S(U) =

∫

U

ρη dv and H(U) =

∫

∂U

h

θ
da+

∫

U

ρr

θ
dv +

∫

U

ργ dv . (2.65)

Here, η is the specific entropy which represents an entropy per unit mass, h = 〈−q,n〉
the heat flux, r the specific heat source power, θ > 0 the absolute temperature and γ
the specific entropy production rate. For simplicity, we assume that the entropy flux
and entropy supply are only related to heat flux and heat source, and not e.g. to mass
transport. Then, the balance of entropy states

d

dt
S(U) = H(U) . (2.66)

Again, applying Reynold’s transport theorem together with the spatial form (2.48) of the
mass balance yields after localization the spatial form of the entropy balance

ργ = ρη̇ −
1

θ
(ρr − div q)−

1

θ2
〈q,∇θ〉 . (2.67)

Here, η̇ is the material time derivative of the specific entropy. In conceptually the same
manner as before, the material local form of the entropy balance is obtained by multiplying
(2.67) by the Jacobian J , reparametrizing the fields η0 = η ◦ ϕ, Θ = θ ◦ ϕ, r0 = r ◦ ϕ,
γ0 = γ ◦ ϕ, using the relationship ∇θ = F −T∇0Θ between the material and spatial
temperature gradients, recalling the definition (2.43) of the material heat flux vector and
using Piola’s identity (2.62). The result is

ρ0γ0 = ρ0η̇0 −
1

Θ
(ρ0r0 −DivQ)−

1

Θ2
〈Q,∇0Θ〉 . (2.68)

Since ρ̇0 = 0, we can rewrite this equation as

γ̃0 = ˙̃η0 −
1

Θ
(r̃0 −DivQ)−

1

Θ2
〈Q,∇0Θ〉 (2.69)

in terms of the entropy rate, heat source power and entropy production rate ˙̃η0, r̃0 and
γ̃0 per unit undeformed volume. Recalling the energy equation (2.64), we can reformulate
the material form (2.69) of the balance of entropy

D = gP : Ḟ +Θ˙̃η0 −
˙̃e0 −

1

Θ
〈Q,∇0Θ〉 (2.70)

in terms of the dissipation D = Θγ̃0 per unit undeformed volume. This equation suggests
a parametrization of the internal energy by the deformation and entropy. However, since
latter variable is difficult to handle practically, one introduces the Helmholtz free-energy
ψ = ẽ0 −Θη̃0 per unit undeformed volume. Then, (2.70) takes the form

D = gP : Ḟ − Θ̇η̃0 − ψ̇ −
1

Θ
〈Q,∇0Θ〉 (2.71)

which suggests a parametrization of the Helmholtz free-energy by the deformation and
temperature.


18 2 Continuum Mechanics and Material Theory

2.4. Concepts of Material Theory

The kinematic relations in Section 2.1 are material independent and so are, up to a
certain point4, the balance equations presented in Section 2.3. In the following section,
we repeat basic principles for the setup of material dependent equations which are needed
for the closure of the continuum-thermodynamical framework.

2.4.1. Principle of Equipresence. A quantity that arises as independent variable
in one constitutive law should also be present in all other constitutive equations unless
physical laws or invariance principles forbid it. As independent variables we identify the
deformation and the temperature.

2.4.2. Principle of Determinism. It is assumed that the thermomechanical mate-
rial response at (X, t) ∈ B × T depends in the most general sense on the full history of
deformation and temperature at all points Y ∈ B. For example, the stress

P (X, t) =
t

P
τ=0

(ϕ(Y , τ),Θ(Y , τ), g(y, τ) ◦ϕ−1(Y , τ);X, t) (2.72)

is determined by the constitutive functional P. Note that we include the spatial metric
tensor g as an additional argument5 which we do for geometric consistency, see the Doyle-
Ericksen formula (2.87). Taking a look at the entropy balance (2.71), we need additional
constitutive functionals that determine at (X, t) ∈ B × T the entropy η̃0, the Helmholtz
free-energy ψ and the heat flux vector Q.

2.4.3. Concept of Internal Variables. It is not practical to calculate e.g. a stress
response at point (X, t) ∈ B×T via (2.72) based on the full history of the deformation and
temperature. Alternatively, a possible path dependence of material behavior is accounted
for by an internal variable field q : B × T → R

δ such that

P (X, t) =P∗
(ϕ(Y , t),Θ(Y , t), g(y, t) ◦ϕ−1(Y , t),q(Y , t);X, t) . (2.73)

The array q has in total δ scalar-valued entries and consists of objects of any tensorial rank
which are defined in the reference configuration. Likewise, we replace the dependency on
the full history of deformation and temperature in all the other constitutive functionals by
including this internal variable field. Note that for all members of q we need to formulate
evolution equations which must fulfill the second law of thermodynamics, see below.

2.4.4. Principle of Local Action. The thermomechanical material response at
(X, t) ∈ B×T is assumed to depend on variables only in an open neighborhoodNX ⊂ B of
X ∈ B. We do Taylor expansions of the deformation, temperature and internal variables,
for example

ϕ(Y , t) = ϕ(X, t) +Dϕ(X, t) (Y −X) + o(‖Y −X‖2G) . (2.74)

For a material of grade n with n ≥ 1, we specify the neighborhood NX by truncating
these Taylor series after the n+ 1 terms. Then, in case of a material of grade one n = 1,
the stress at (X, t) ∈ B × T is governed by a stress response function

P (X, t) = P̂ (ϕ,F ,Θ,∇0Θ, g,q,∇0q;X, t) . (2.75)

4For example, the balance of angular momentum changes drastically if we consider additional couple
stresses.

5In our case, the spatial metric g is Euclidean and hence neither a function of space nor time.


2.4 Concepts of Material Theory 19

Likewise, we define entropy and heat flux response functions together with the Helmoltz
free-energy function

ψ(X, t) = ψ̂(ϕ,F ,Θ,∇0Θ, g,q,∇0q;X, t) . (2.76)

Latter plays a very important role in material modeling as shown below.

2.4.5. Principle of Objectivity. We consider a rigid body motion superimposed on
the current configuration S yielding a new configuration denoted by S+. The Euclidean
metric tensor defined on this new configuration is denoted by g+ : Tx+S+ → Tx+S+ and
we choose on the new manifold S+ a coordinate system {xα}. The deformation map

ϕ+(X, t, τ) = Ψ(x, τ) ◦ϕ(X, t) with Ψ(x, τ) = c(τ) +R(τ)x (2.77)

maps at time t+ = t + τ , τ > 0 points X ∈ B from the reference configuration to points
x+ ∈ S+ of the new current configuration S+. Here, c denotes a translation vector and
R ∈ SO+(ω) a rotation tensor which is a member of the special orthogonal group

SO+(ω) =
{
R : TxS → Tx+S+ | g+ = R−TgR−1 and detR = 1

}
. (2.78)

Since the reference configuration B is not affected by the rigid body motion (2.77), material
quantities such as the material temperature gradient ∇0Θ ∈ T ∗

XB or the introduced
internal variables q together with their gradient ∇0q are not influenced. In contrast,
two-point tensors such as the deformation gradient transform by their first index slot

F+ = RF or (F+)αA = Rα
a F

a
A . (2.79)

The principle of objectivity now demands, that a scalar valued state of the material such
as energy or entropy must not be changed by the rigid body motion, for example

ψ̂(ϕ+,F+,Θ,∇0Θ, g
+,q,∇0q;X, t+) = ψ̂(ϕ,F ,Θ,∇0Θ, g,q,∇0q;X, t) (2.80)

for all R ∈ SO+(ω) and τ > 0. This condition imposes important restrictions on the
Helmholtz free-energy function, i.e. it rules out the explicit dependence on (i) the de-
formation ϕ and (ii) the time t. The same observation holds for the entropy response
function. For the tensor-valued stress response function (2.75), the principle of objectivity
demands

P̂ (ϕ+,F+,Θ,∇0Θ, g
+,q,∇0q;X, t+) = RP̂ (ϕ,F ,Θ,∇0Θ, g,q,∇0q;X, t) (2.81)

for all R ∈ SO+(ω) and τ > 0, since its image, the first Piola-Kirchhoff stress tensor, is a
two-point tensor that transforms with respect to the first index slot. Finally, the principle
of objectivity stipulates for the heat flux response function

Q̂(ϕ+,F+,Θ,∇0Θ, g
+,q,∇0q;X, t+) = Q̂(ϕ,F ,Θ,∇0Θ, g,q,∇0q;X, t) (2.82)

for all R ∈ SO+(ω) and τ > 0, since its image, the heatflux Q ∈ TXB, is a material
object. Like for the Helmholtz free-energy function, we see from (2.82) that the heat flux
response function cannot explicitely depend on the deformation ϕ and time t.


20 2 Continuum Mechanics and Material Theory

2.4.6. Second Law of Thermodynamics. As a further fundamental restriction, we
must ensure that the dissipation (2.70) is always nonnegative D ≥ 0. This is known as the
second law of thermodynamics. It is common to impose an even stronger condition which
demands nonnegativity on the intrinsic and conductive parts of the dissipation seperately

Dint = gP : Ḟ − Θ̇η̃0 − ψ̇ ≥ 0 and Dcon = −
1

Θ
〈Q,∇0Θ〉 ≥ 0 (2.83)

where (2.83)1 is called the Clausius-Planck inequality and (2.83)2 the Fourier inequality.
Taking the time derivative of the Helmholtz free-energy function (2.76) and considering
the result stemming from the principle of objecticity together with ġ = 0 , we can rewrite
the Clausius-Planck inequality in the form

Dint = [ gP − ∂F ψ̂ ] : Ḟ − [ η̃0 + ∂Θψ̂ ]Θ̇− 〈∂∇0Θψ̂,∇0Θ̇〉 − 〈∂qψ̂, q̇〉 − 〈∂∇0qψ̂,∇0q̇〉 ≥ 0 .
(2.84)

From this inequality, the standard Coleman-Noll argument for unconstrained materials
concludes6

gP = ∂F ψ̂ , η̃0 = −∂Θψ̂ and ∂∇0Θψ̂ = 0 . (2.85)

Hence, the stress and entropy response functions are governed by the Helmholtz free-
energy function. In addition, (2.85)3 excludes the dependency of the Helmholtz free-energy
on the temperature gradient and we finally get the reduced representation

ψ(X, t) = ψ̂(F ,Θ, g,q,∇0q;X) . (2.86)

The inclusion of the spatial metric g as an argument allows to express the Kirchhoff stress
via the formula of Doyle-Ericksen

τ = 2∂gψ̂ or τab = 2
∂ψ̂

∂gab

∣∣∣∣∣
gab=δab

(2.87)

which relies on the spatial representation of the stress power per unit undeformed volume
(2.39) and again on the Coleman-Noll argument for unconstrained materials. Observe,
that according to the principle of objectivity the Helmholtz free-energy function must
fulfill the invariance requirement

ψ̂(RF ,Θ,R−TgR,q,∇0q;X) = ψ̂(F ,Θ, g,q,∇0q;X) (2.88)

for all R ∈ SO(ω). This condition is fulfilled a priori by a representation in terms of the
right Cauchy-Green tensor

ψ(X, t) = ψ̃(C,Θ,q,∇0q;X) . (2.89)

Then, taking into account the material representation of the stress power per unit unde-
formed volume (2.39), the second Piola-Kirchhoff stress tensor follows from the standard
Coleman-Noll argument for unconstrained materials as

S = 2∂Cψ̃ . (2.90)

6Examples for internal material constraints are inextensibility and incompressibility, see Ogden [41].


2.4 Concepts of Material Theory 21

It remains the reduced Clausius-Planck inequality

Dint = −〈∂qψ̂, q̇〉 − 〈∂∇0qψ̂,∇0q̇〉 ≥ 0 (2.91)

which poses a restriction on evolution laws for the internal variables q. To fulfill this
inequality a priori, one often makes use of so called dissipation potential functions. For
rheological models (point problems) this approach is presented in Chapter 3. For a treat-
ment of the inequality (2.91) including the gradient term, we refer to Miehe [30].

Finally, a heat flux response function is constructed such that Fourier’s inequality
(2.83)2 is fulfilled. For the Kirchhoff heat flux vector Jq, we choose a Fourier-type law

Jgq = −κ∇θ or Jδabq
b = −κ θ,a (2.92)

in terms of the heat conduction coefficient κ > 0. Recalling definition (2.43), the material
heat flux response function reads

Q̂(F , g,∇0Θ) = −κC−1∇0Θ or QA = −κ(C−1)AB Θ,B (2.93)

which conforms with the objectivity demand (2.82). Since the right Cauchy-Green tensor
is positive definite, this material law of heat conduction fulfills Fourier’s inequality (2.83)2.

2.4.7. Principle of Material Symmetry. Materials like crystals possess certain
symmetries due to their microstructures. A symmetry is mathematically described by a
symmetry group G which consists of all rotations that leave the microstructure unchanged.
A list of the generators of these groups can be found in Truesdell & Noll [48] for all
crystal systems and transverse isotropy. In this sense, we consider a rotation superimposed
on an infinitesimal subbodyNX ⊂ B of the reference configuration. Hence, tensors defined
on the current configuration are unaffected. It is important to note, that with such a
local concept, a new rotated reference configuration cannot be introduced. We define the
modified tangent map

F ∗ = FR−1 or (F ∗)aA = F a
B (R−1)BA (2.94)

for any rotation tensor R that is a member of the special orthogonal group

SO(ω) = {R : TXB → TXB |G = RTGR and detR = 1} . (2.95)

Then, for a given symmetry group G ⊆ SO(ω), the principle of material symmetry de-
mands7

ψ̂(F ∗,Θ, g,R ⋆ q,R ⋆∇0q;X) = ψ̂(F ,Θ, g,q,∇0q;X) (2.97)

for all R ∈ G. It means that the energetic state must not change if the material element
is rotated by any member of its symmetry group before deformation. According to the

7Rayleigh product. For example, for a covariant tensor T of rank n ≥ 0, the Rayleigh product with
a mixed-variant second-order tensor R ∈ TXB ⊗ T ∗

X
B is defined as

R ⋆ T = TA1...An
(R−TEA1)⊗ ...⊗ (R−TEAn) . (2.96)


22 2 Continuum Mechanics and Material Theory

isotropication theorem, the anisotropic scalar-valued tensor function (2.97) can be made
isotropic by introducing a structural tensor M as additional argument

ψ̂(F ∗,Θ, g,R ⋆ q,R ⋆∇0q;R ⋆M ,X) = ψ̂(F ,Θ, g,q,∇0q;M ,X) (2.98)

for all R ∈ SO(ω), where M models the microstructural symmetry via the invariance
requirement R ⋆M =M for all R ∈ G. A list of invariants forming the basis of isotropic
scalar-valued tensor functions can be found in e.g. Boehler [6].


23

3. Variational Principles for Rheological Models

In this chapter, we consider an elementary characterization of mechanical material
responses. We think about a material element that is loaded within a time interval
T = (0, te) by a stress σ(t) and deforms by a strain ε(t). With the help of experiments,
we can basically distinguish four types of material behavior, see for example Haupt [18]:

1. Elasticity describes rate-independent material behavior without hysteresis.

2. Elasto-plasticity describes rate-independent material behavior with hysteresis.

3. Visco-elasticity describes rate-dependent material behavior and the equilibrium curve,
that is reached by a relaxation process, does not show a hystersis.

4. Visco-elasto-plasticity describes rate-dependent material behavior and the equilib-
rium curve shows a hysteresis.

These four types of material behavior may be described by serial and parallel combinations
of three basic material devices. These are the spring characterizing elastic behavior, the
dashpot characterizing viscous behavior and the friction element which models plastic
behavior, see Figure 4. Especially, the material models of the three devices are

• Hooke’s law σe = E εe

• Newton’s law σv = H ε̇v

• St. Venant’s law

{
|σp| ≤ y0 for ε̇p = 0
|σp| = y0 for ε̇p 6= 0 ,

with Young’s modulus E > 0, viscosity H > 0 and yield stress y0 > 0. Each of these
devices is now analyzed in more detail, where a specific focus is put on variational princi-
ples. Once these are set up, more complicated material models based on these three basic
devices can be formulated quite easily. The basis for this chapter are the lecture notes of
Miehe [29].

3.1. Elementary Formulation of Basic Rheological Devices

We now analyze the single spring, dashpot and friction elements. To do so, we consider
within a time interval (0, t) ⊆ T the work done on the material element

W t
0 =

∫ ε(t)

ε(0)

σ dε̄ =

∫ t

0

P dt̄ with P = σε̇ (3.1)

in terms of the stress power P.

3.1.1. Variational Principles for Elastic Device. A material is called (hyper)-
elastic if the process it undergoes is reversible, i.e. if in a closed strain cycle the work
(3.1) is zero. Hence, for any t ∈ T there must exist a potential ψ(t) such that (3.1) is
expressed by the difference

W t
0 = ψ̂(ε(t))− ψ̂(ε(0)) (3.2)


24 3 Variational Principles for Rheological Models

ε

ε

ε

σ

σ

σ

σ

σ

σ

E

H

y0

a)

b)

c)

Figure 4: Three basic material devices. a) Spring characterizing elastic response via
Hooke’s law, b) dashpot modeling viscous response via Newton’s law and c) friction device
characterizing plastic response via St. Venant’s law.

where ψ̂(ε(·)) is called the stored-energy function. Thus, we have for elasticity the balance

P =
d

dt
ψ̂(ε) (3.3)

or by using the chain rule

σ ε̇−
d

dt
ψ̂(ε) = [ σ − ∂εψ̂ ] ε̇ = 0 . (3.4)

Since this equation has to hold for any strain rate ε̇, we can conclude for the stress

σ(t) = ∂εψ̂(ε(t)) . (3.5)

The stored-energy function must fulfill the conditions of (i) normalization ψ̂(0) = 0 and

∂εψ̂(0) = 0 for a stress-free initial state and (ii) strict convexity

ψ̂(aε1 + (1− a)ε2) < aψ̂(ε1) + (1− a)ψ̂(ε2) for all ε1 6= ε2, a ∈ (0, 1) (3.6)

for a stable (one-to-one) elastic material response. An equivalent statement of the strict
convexity condition (3.6) is

ψ̂(ε2) > ψ̂(ε1) + ∂εψ̂(ε1)(ε2 − ε1) for all ε1 6= ε2 . (3.7)

With such a stored-energy function at hand, we define the potential

π̂e(ε; t) = ψ̂(ε)− σext(t)ε (3.8)

associated with a spring device that is loaded by an external stress σext(t). Then, the
strain at time t follows by minimizing this potential

εt = Arg{ inf
ε
π̂e(ε; t) } . (3.9)


3.1 Elementary Formulation of Basic Rheological Devices 25

From the strict convexity condition (3.7), it follows that a sufficient condition for the
existence and uniqueness of a global minimizer εt, i.e. π̂e(εt) < π̂e(ε) for all ε 6= εt, is a
vanishing first derivate of the potential evaluated at εt. We have the Euler equation

∂επ̂e(εt; t) ≡ ∂εψ̂(εt)− σext(t) = 0 (3.10)

which is just the stress equilibrium in the spring.8

To derive a principle alternative to the primal formulation (3.9), we consider the
Legendre transformation of the stored-energy function

ψ̂∗(σ) = sup
ε
[ σε− ψ̂(ε) ] or ψ̂(ε) = sup

σ
[ σε− ψ̂∗(σ) ] (3.11)

which can be understood as a change of variable by the energetically conjugate quantity.
ψ̂∗(·) is the dual function of the stored energy ψ̂(·) and is refered to as complementary

energy. Due to the strict convexity of ψ̂(·), the complementary energy is also strictly
convex. The optimality conditions of (3.11) define in an inverse format the strain in
terms of the stress and vice versa

σ = ∂εψ̂(ε) and ε = ∂σψ̂
∗(σ) . (3.12)

As an example, we look at the simple quadratic stored-energy function ψ̂(ε) = 1/2Eε2

with Young’s modulus E > 0 and determine the complementary energy via (3.11)1. The
corresponding optimality condition (3.12)1 results in Hooke’s law σ = Eε and solving for

ε yields, when inserting in σε− ψ̂(ε), the complementary energy ψ̂∗(σ) = 1/(2E)σ2.

Recalling the potential (3.8), the Legendre transformation (3.11)2 now motivates the
definition of the mixed potential

π̂∗
e(ε, σ; t) = σε− ψ̂∗(σ)− σext(t)ε (3.13)

in terms of the complementary energy. Then, the strain and the stress at time t are
governed by the mixed saddle point principle

{ εt, σt } = Arg{ inf
ε

sup
σ

π̂∗
e(ε, σ; t) } . (3.14)

We obtain as sufficient optimality conditions the Euler equations

1. ∂επ̂
∗
e(εt, σt; t) ≡ σt − σext(t) = 0

2. ∂σπ̂
∗
e(εt, σt; t) ≡ εt − ∂σψ̂

∗(σt) = 0
(3.15)

which represent the stress equilibrium and the inverse constitutive definition of the stress.
The mixed saddle point principle (3.14) is generally known as the Hellinger-Reissner
principle of elasticity and falls into the category of primal-dual formulations.

8A vanishing first derivative of the potential is always a necessary condition for local unconstrained
optimality as can be shown by doing a Taylor expansion. If the potential is convex, a vanishing first
derivative is a sufficient condition for the existence of (at least) one minimizer.


26 3 Variational Principles for Rheological Models

ε̇ε̇ε̇

φφφ

a) b) c)

Figure 5: Dissipation potential functions. a) Function satisfying conditions (i)–(iii), b)
function not fulfilling the condition (ii) of nonnegativity and violating the dissipation in-
equality (3.24) in the grey region and c) function violating the convexity condition (iii) but
satisfying the dissipation inequality (3.24) everywhere.

To obtain a principle that is dual to the primal formulation (3.9), we assume the order
of the supremum and infimum in (3.14) to be interchangeable9

inf
ε

sup
σ

π̂∗
e(ε, σ; t) = sup

σ
inf
ε
π̂∗
e(ε, σ; t) = sup

σ=σext(t)

[−ψ̂∗(σ) ] , (3.17)

where in the last equality we took note that infε π̂
∗
e(ε, σ; t) = −∞ unless σ = σext(t). The

constrained maximization formulation in (3.17) is known as Castigliano’s principle. So far,
the elastic material device was driven by an externally applied stress σext(t). Alternatively,
we can load the material element by controlling the strain εext(t) and define the potential

π̃∗
e(σ; t) = ψ̂∗(σ)− σεext(t) (3.18)

in terms of the complementary energy. Then, the stress at time t is determined via the
minimization principle

σt = Arg{ inf
σ
π̃∗
e(σ; t) } . (3.19)

The corresponding Euler equation is

∂σπ̃
∗
e(σt; t) ≡ ∂σψ̂

∗(σt)− εext(t) = 0 (3.20)

which defines the stress in an inverse constitutive format.

3.1.2. Variational Principles for Viscous Device. A viscous material is charac-
terized by a fully dissipative response which is rate-dependent. The (dissipative) stress
σ is in a one-to-one relationship with the strain rate ε̇. We consider the work W t

0 done
on the dashpot device as defined in (3.1) which is fully dissipated, that means converted
into heat. This is expressed by the balance

W t
0 = Dt

0 with Dt
0 =

∫ t

0

D dt̄ (3.21)

9Note carefully, that this assumption is not always true. In general, it just holds that

inf
ε
sup
σ

π̂∗
e (ε, σ; t) ≥ sup

σ

inf
ε
π̂∗
e(ε, σ; t) . (3.16)


3.1 Elementary Formulation of Basic Rheological Devices 27

in terms of the dissipation D. Thus, the stress power P is identical to the dissipation

σε̇ = D ≥ 0 (3.22)

which must be nonnegative according to the second law of thermodynamics. Following
Biot [5], Moreau [39] and Halphen & Nguyen [17], we introduce a (viscous) dissi-

pation potential function φ̂(ε̇(·)) that determines the stress in terms of the strain rate

σ(t) = ∂ε̇φ̂(ε̇(t)) . (3.23)

Note, that for viscous material behavior the function φ̂(·) is smooth. The dissipation
inequality (3.22) takes the form

D = ∂ε̇φ̂(ε̇) ε̇ ≥ 0 (3.24)

which is a priori fulfilled for a dissipation potential function that is (i) normalized φ̂(0) = 0,

(ii) nonnegative φ̂(·) ≥ 0 and (iii) convex

φ̂(aε̇1 + (1− a)ε̇2) ≤ aφ̂(ε̇1) + (1− a)φ̂(ε̇2) for all (ε̇1, ε̇2), a ∈ [0, 1] . (3.25)

An equivalent statement of the convexity condition (3.25) is

φ̂(ε̇2) ≥ φ̂(ε̇1) + ∂ε̇φ̂(ε̇1) (ε̇2 − ε̇1) for all (ε̇1, ε̇2) . (3.26)

In addition, we may demand ∂ε̇φ̂(0) = 0 in order to have a stress-free reference state. The
sufficiency of the conditions (i)–(iii) for the fulfillment of the dissipation inequality (3.24)
is seen by setting ε̇2 = 0 in (3.26) yielding

∂ε̇φ̂(ε̇1) ε̇1 ≥ φ̂(ε̇1) ≥ 0 for all ε̇1 . (3.27)

Note that the convexity condition is unnecessarily strong. From the dissipation inequality
(3.24) it just follows that within any connected interval starting from the origin ε̇ = 0,
the dissipation potential function must be monotonic but can be nonconvex, see Figure 5.
In this case, the relation between the stress and strain rate is not one-to-one and globally
as well as locally stable states occur that are global and local minimizers of the potential
(3.28) defined below, see Ericksen [12].

With such a dissipation potential function at hand, we define the rate-type potential

π̂v(ε̇; t) = φ̂(ε̇)− σext(t)ε̇ (3.28)

associated with a dashpot device that is loaded by an external stress σext(t). Then, the
strain rate at time t follows by minimizing this potential

ε̇t ∈ Arg{ inf
ε̇
π̂v(ε̇; t) } . (3.29)

Taking φ̂(·) to be convex, the corresponding Euler equation

∂ε̇π̂v(ε̇t; t) ≡ ∂ε̇φ̂(ε̇t)− σext(t) = 0 (3.30)


28 3 Variational Principles for Rheological Models

is a sufficient condition for ε̇t to be a global minimizer of the rate-type potential (3.29)
and is just the stress equilibrium in the device. From a formal viewpoint, multiple mini-
mizers may exist which is expressed by the element sign in (3.29). However, for a viscous
material response the dissipation is positive since a nonzero evolution of the strain sets in
immediately at loading (in contrast to plasticity, see below). With regard to the notion of
convexity, this leads to the definition of a strictly convex (viscous) dissipation potential
function and the element sign in (3.29) can be replaced by an equality sign.

To derive a mixed principle of Hellinger-Reissner type, we consider the Legendre trans-
formation of the dissipation potential function

φ̂∗(σ) = sup
ε̇

[ σε̇− φ̂(ε̇) ] or φ̂(ε̇) = sup
σ
[ σε̇− φ̂∗(σ) ] (3.31)

and call φ̂∗(·) the dual dissipation potential function. The corresponding optimality con-
ditions define in an inverse format the strain rate in terms of the stress and vice versa

σ = ∂ε̇φ̂(ε̇) and ε̇ = ∂σφ̂
∗(σ) . (3.32)

As an example, we consider the quadratic dissipation potential function φ̂(ε̇) = H/2 ε̇2

with viscosity H > 0. The optimality condition (3.32)1 results in Newton’s law σ = Hε̇

and solving for ε̇ yields, when inserting into σε̇ − φ̂(ε̇), the dual dissipation potential

function φ̂∗(σ) = 1/(2H) σ2.

Recalling the rate-type potential (3.28), the Legendre transformation (3.31)2 now mo-
tivates the definition of the mixed rate-type potential

π̂∗
v(ε̇, σ; t) = σε̇− φ̂∗(σ)− σext(t)ε̇ (3.33)

in terms of the dual dissipation potential function. Then, the strain rate and the stress
at time t are governed by the mixed saddle point principle

{ ε̇t, σt } = Arg{ inf
ε̇

sup
σ

π̂∗
v(ε̇, σ; t) } . (3.34)

The sufficient optimality conditions are the Euler equations

1. ∂ε̇π̂
∗
v(ε̇t, σt; t) ≡ σt − σext(t) = 0

2. ∂σπ̂
∗
v(ε̇t, σt; t) ≡ ε̇t − ∂σφ̂

∗(σt) = 0
(3.35)

which include the stress equilibrium and the evolution equation for the (viscous) strain.

Like for the elastic device, we have a dual formulation based on the potential

π̃∗
v(σ; t) = φ̂∗(σ)− σε̇ext(t) (3.36)

where we control the strain rate of the viscous device. Then, the stress at time t is
governed by the minimization principle

σt = Arg{ inf
σ
π̃∗
v(σ; t) } . (3.37)

The corresponding Euler equation is

∂σπ̃
∗
v(σt; t) ≡ ∂σφ̂

∗(σt)− ε̇ext(t) = 0 (3.38)

which defines the stress in an inverse constitutive format.


3.1 Elementary Formulation of Basic Rheological Devices 29

ε̇ ε̇

φ σ

E

1

y0 y0

−y0
a) b)

Figure 6: Plastic device. a) Dissipation potential function (3.41) that is positively homo-
geneous of degree one and nonsmooth at ε̇ = 0. The set of all slopes of lines constructed
at (0, φ̂(0)) that lie under the graph define the subgradient ∂ε̇φ̂(0) which is called elastic
domain E. b) Set-valued stress law.

3.1.3. Variational Principles for Friction Device. A plastic material behavior is
characterized by a fully dissipative response which is rate-independent. The (dissipative)
stress σ is not in a one-to-one relationship with the strain rate ε̇. As before, we have the
balance

σε̇ = D ≥ 0 . (3.39)

The difference to the viscous device is the choice of a nonsmooth dissipation potential
function that is positively homogeneous of degree one, that means

anφ̂(ε̇) = φ̂(aε̇) with n = 1 for all a > 0 . (3.40)

The classical form of a dissipation potential function representing plastic behavior is

φ̂(ε̇) = y0| ε̇ | (3.41)

where y0 > 0 is the yield stress. This function fulfills the conditions (i)–(iii) from Sec-
tion 3.1.2 and is nonsmooth at ε̇ = 0, see Figure 6a). For nonsmooth functions, the
classical definition of derivatives must be generalized. This leads to the introduction of
the subgradient which at a nonsmooth point is not just a value but a set of values, see
Rockafellar [45]. The subgradient of the dissipation potential function φ̂(ε̇) evaluated
at ε̇ = 0 is the set

E = ∂ε̇φ̂(0) = { σ | σε̇ ≤ φ̂(ε̇) for all ε̇ } . (3.42)

It contains the slopes σ of all lines which go through the point (0, φ̂(0)) and lie below the

graph of φ̂(ε̇), see Figure 6a). The set E is called elastic domain and restricts the range
of the stresses, i.e.

E = ∂ε̇φ̂(0) = {σ | − y0 ≤ σ ≤ y0} . (3.43)

In addition, we define the open range

int(E) = {σ | − y0 < σ < y0} (3.44)


30 3 Variational Principles for Rheological Models

such that E = int(E)∪∂E with ∂E = {−y0, y0}. We get as full derivative of the dissipation
potential function (3.41)

∂ε̇φ̂(ε̇) =





−y0 for ε̇ < 0
E for ε̇ = 0
y0 for ε̇ > 0 .

(3.45)

In contrast to the viscous device discussed in Section 3.1.2 that is characterized by a
smooth evolution of the strain, the stress for a plastic device is defined via a differential
inclusion

σ ∈ ∂ε̇φ̂(ε̇) (3.46)

which is called a set-valued stress law. In particular, we obtain St. Venant’s law

{
σ ∈ E for ε̇ = 0
σ = y0 ε̇/| ε̇ | for ε̇ 6= 0

(3.47)

that is depicted in Figure 6b) and in an inverse form in Figure 7b). Then, the dissipation
in a plastic device results in

D = σε̇ = y0| ε̇ | (3.48)

which is equal to the image of the dissipation potential function. This property does not
hold for the viscous device and is ultimately related to the positive homogeneity of degree
one of the dissipation potential function (3.41). Especially, taking the derivative of (3.40)
with respect to a and evaluating at a = 1 yields10

nφ̂(ε̇) = ∂ε̇φ̂(ε̇) ε̇ with n = 1 . (3.49)

With the dissipation potential function (3.41) at hand, we define like for the viscous device
the rate-type potential

π̂p(ε̇; t) = φ̂(ε̇)− σext(t)ε̇ (3.50)

associated with a friction element that is loaded by an external stress −y0 ≤ σext(t) ≤ y0.
Then, the strain rate at time t is governed by the minimization principle

ε̇t ∈ Arg{ inf
ε̇
π̂p(ε̇; t) } . (3.51)

The sufficient optimality condition is the set-type Euler equation

∂ε̇π̂p(ε̇t; t) ≡ ∂ε̇φ̂(ε̇t)− σext(t) ∋ 0 . (3.52)

The solution ε̇t of the variational principle (3.51) is not unique for |σext(t)| = y0, i.e.
ε̇t ∈ R≥0 for σext(t) = y0 and ε̇t ∈ R≤0 for σext(t) = −y0, respectively. Here, ε̇t = 0
corresponds to unloading from a plastic state. We now want to define rate-independency

10For a viscous device we classically choose a quadratic dissipation potential function, see Section 3.1.2.
Such a function is positively homogeneous of degree two, i.e. set n = 2 in (3.40) and (3.49), and its image
represents half the dissipation.


3.1 Elementary Formulation of Basic Rheological Devices 31

in a precise format. Consider a rescaling of time t 7→ γ(t) where γ(·) is a continuous

and strictly increasing function. We write
◦

ε = dε/dγ and have ε̇ =
◦

ε γ̇. Recalling the
positive homogeneity of degree one (3.40) with the property (3.49), we can rewrite the
Euler equation (3.52) in the form

nφ̂(ε̇) = σext(t)ε̇ ⇐⇒ nγ̇n−1φ̂(
◦

ε) = σext(γ)
◦

ε with n = 1 . (3.53)

Hence, rescaling time does not influence the evolution. This property characterizes rate-
independent systems and is achieved by a dissipation potential function that is positively
homogeneous of degree one.11 For an extensive mathematical treatment on general rate-
independent systems we refer to Mielke & Roub́ıcek [38].

The dual dissipation potential function is obtained by the Legendre transformation of
the nonsmooth function (3.41) and reads

φ̂∗(σ) = sup
σ
[ σε̇− φ̂(ε̇) ] =

{
0 for σ ∈ E

+∞ else .
(3.54)

It is the indicator function of the set (3.43) of admissible stresses and is depicted in
Figure 7a). The subgradient of the dual dissipation potential function is the set

∂σφ̂
∗(σ) =





R≤0 for σ = −y0
0 for σ ∈ int(E)
R≥0 for σ = y0
∅ else

(3.55)

and determines the evolution of the strain via the differential inclusion

ε̇ ∈ ∂σφ̂
∗(σ) (3.56)

which can also be viewed as the constitutive definition of the stress in terms of the strain
rate in an inverse format, see Figure 7b).

The Legendre transformation (3.54) motivates the definition of the mixed rate-type
potential

π̂∗
p(ε̇, σ; t) = σε̇− φ̂∗(σ)− σext(t)ε̇ . (3.57)

Then, the strain rate and stress at time t are governed by the mixed saddle point principle

{ ε̇t, σt } ∈ Arg{ inf
ε̇

sup
σ

π̂∗
p(ε̇, σ; t) } . (3.58)

The corresponding Euler equations read

1. ∂ε̇π̂
∗
p(ε̇t, σt; t) ≡ σt − σext(t) = 0

2. ∂σπ̂
∗
p(ε̇t, σt; t) ≡ ε̇t − ∂σφ̂

∗(σt) ∋ 0 .
(3.59)

Like for the elastic and viscous devices, we have a dual formulation based on the potential

π̃∗
p(σ; t) = φ̂∗(σ)− σε̇ext(t) (3.60)

11For a positively homogeneous dissipation potential function of degree n 6= 1, we clearly observe a
change in the evolution by rescaling time.


32 3 Variational Principles for Rheological Models

where we control the strain rate of the friction element. Then, the stress at time t is
determined by the minimization principle

σt ∈ Arg{ inf
σ
π̃∗
p(σ; t) } . (3.61)

The corresponding Euler equation is

∂σπ̃
∗
p(σt; t) ≡ ∂σφ̂

∗(σt)− ε̇ext(t) ∋ 0 (3.62)

which defines the stress in an inverse format, i.e. σt ∈ int(E) for ε̇ext(t) = 0, σt = y0 for
ε̇ext ≥ 0 and σt = −y0 for ε̇ext ≤ 0.

In plasticity, constitutive modeling is often done by defining yield functions f̂(·) in
stress space which then govern the dissipation potential function. We start with modeling
the elastic domain

E = {σ | f̂(σ) ≤ 0} or E = {σ | χ̂(σ) ≤ y0} (3.63)

where y0 is a threshold value which before was denoted as yield stress. The level set
function χ̂(·) is a gauge, i.e. it is (like the dissipation potential function) (i) normalized,
(ii) nonnegative, (iii) convex and (iv) positively homogeneous of degree one. For our case,
we have

χ̂(σ) = | σ | (3.64)

and the yield function is defined by

f̂(σ) = | σ | − y0 . (3.65)

Then, the dissipation potential function is determined via the principle of maximum
plastic dissipation

φ̂(ε̇) = sup
σ∈E

[ σε̇ ] , (3.66)

which identifies for given strain rate ε̇ the stress σ as the one for which

σε̇ ≥ σ∗ε̇ for all σ∗ ∈ E , (3.67)

see e.g. Lubliner [23]. Cleary, a function φ̂(·) constructed via (3.66) is positively homo-
geneous of degree one. Due to the convexity of χ̂(·), the constrained optimization problem
(3.66) has a unique solution. Using the Lagrange multiplier method, we can rewrite (3.66)
in the form

φ̂(ε̇) = sup
σ

inf
λ≥0

[ σε̇− λ(χ̂(σ)− y0) ] . (3.68)

The optimization with respect to the stress yields the evolution equation

ε̇ = λ ∂σχ̂(σ) (3.69)


3.1 Elementary Formulation of Basic Rheological Devices 33

ε̇φ∗

E

∂σφ̂
∗(−y0) ∂σφ̂

∗(y0)

+∞+∞

σ σy0−y0

a) b)

Figure 7: Friction device. a) Dual dissipation potential function (3.54) which is the indi-
cator function of the elastic domain E and its subgradient. b) Inverse St. Venant’s law.

which is generally known as flow rule. Since ∂σχ̂(σ) = σ/|σ|, the Lagrange parameter
λ characterizes the amount of plastic flow. To obtain the optimality conditions associ-
ated with the minimization with respect to the Lagrange multiplier, consider the Taylor
expansion of the function φ̃(σ, λ) = −λ(χ̂(σ)− y0) for fixed σ, i.e.

φ̃(σ, λ+ δλ) = φ̃(σ, λ) + ∂λφ̃(σ, λ) δλ+ o(|δλ|2). (3.70)

Then, λ is a minimizer iff φ̃(σ, λ+δλ)−φ̃(σ, λ) ≥ 0 for all δλ such that λ+δλ ≥ 0. Letting
|δλ| become very small, we obtain the first-order optimality condition for an unilateral
minimum

∂λφ̃(σ, λ) δλ ≥ 0 (3.71)

for all δλ with λ+ δλ ≥ 0. As a result, we get the classical loading/unloading conditions
in Karush-Kuhn-Tucker form

λ ≥ 0 , χ̂(σ)− y0 ≤ 0 and λ[ χ̂(σ)− y0 ] = 0 . (3.72)

(3.72)1 follows from setting δλ = 0. Assuming that λ = 0, we must have δλ ≥ 0 and

inequality (3.72)2 follows from ∂λφ̃(σ, λ) ≥ 0 according to (3.71). Finally, if λ > 0, δλ

can have any sign and we obtain from (3.71) the condition ∂λφ̃(σ, λ) = 0, or in summary
equality (3.72)3. The Lagrange multiplier is an additional parameter that describes the
nonsmoothness of the formulation in a scalar-valued format. In addition, we mention the
plastic consistency condition

for f̂ = 0 : λ ≥ 0 ,
d

dt
f̂ ≤ 0 , λ

d

dt
f̂ = 0 , (3.73)

see e.g. Simo & Hughes [47]. It can be obtained in the following way: consider at time
t a plastic evolution λt > 0 and the first-order optimality condition (3.71) takes the form

of an equality −f̂ |t δλ = 0 for all δλ. Assuming at time t+ τ , τ ≥ 0 a state λt+τ ≥ 0, we

have −f̂ |t+τ δλ ≥ 0 for all λt+τ + δλ ≥ 0. Subtraction and a Taylor expansion12 yield

1

τ

[
d

dt
f̂ |t τ + o(τ 2)

]
δλ ≤ 0 (3.74)

12We consider a smooth temporal evolution of the yield function f̂ .


34 3 Variational Principles for Rheological Models

for all λt+τ + δλ ≥ 0. Taking the limit τ → 0, we obtain d
dt
f̂ |t = 0 since δλ can have any

sign. On the other hand, for τ small enough we take λt+τ = 0 and we get d
dt
f̂ |t ≤ 0 since

δλ ≥ 0. Collecting the results yields the consistency condition (3.73) which states that
plastic evolution cannot drive the stress state outside of the elastic domain, i.e. σ̇ = 0 for
λ > 0.

Finally, we mention that with the flow rule (3.69), the dissipation follows as

D = σε̇ = λχ̂(σ) ≥ 0 (3.75)

where we made use of the property ∂σχ̂(σ)σ = χ̂(σ) stemming from the positive ho-
mogeneity of degree one of χ̂(·). Refering to (3.75), the image of this function can be
considered as a scalar force that drives the amount λ of plastic flow.

3.2. General Setting for Rheological Models

Rheological models suitable for the description of experimental observations are ob-
tained by a combination of the three basic devices that are discussed in Section 3.1.
Especially, the models are constructed via parallel and serial arrangements of the devices.
The state of the material is described by a minimum number of variables that consist
of the total (external) strain ε and some (strain-like) internal variables {αi}i=1,n. We
summarize these variables in the constitutive state array

c(t) = ( ε(t), α1(t), . . . , αn(t) ) . (3.76)

With this at hand, the strain associated with every single spring, dashpot and friction
device in the rheological model can be found via specific kinematic functions

εi = k̂i(c) for i = 1, . . . , m (3.77)

where m is the number of the devices the model consists of. Next, we define the stored-
energy and dissipation potential functions of the rheological model via the sum of the
single device contributions

ψ̂(c) =
m∑

i=1

ψ̂i(c) and φ̂(ċ) =
m∑

i=1

φ̂i(ċ) , (3.78)

where we assume latter to be a nonsmooth function in all arguments for generality. These
constitutive functions are the basic ingredients for the construction of rate-type variational
principles. We define for a given constitutive state at time t the rate-type potential

π̂(ċ; ct, t) =
d

dt
ψ̂(c) + φ̂(ċ)− σext(t)ε̇ (3.79)

which after applying the chain rule reads

π̂(ċ; ct, t) = ∂cψ̂(ct) · ċ+ φ̂(ċ)− σext(t)ε̇ . (3.80)

Here, the only nonlinearity in ċ occurs due to the dissipation potential function. Then,
the evolution of the constitutive state at time t is governed by the minimization principle

ċt ∈ Arg{ inf
ċ
π̂(ċ; ct, t) } . (3.81)


3.3 Modeling Examples 35

The corresponding Euler equations read

1. ∂ε̇π̂(ċt; ct, t) ≡ ∂εψ̂(ct) + ∂ε̇φ̂(ċt)− σext(t) ∋ 0

2. ∂α̇i
π̂(ċt; ct, t) ≡ ∂αi

ψ̂(ct) + ∂α̇i
φ̂(ċt) ∋ 0

(3.82)

which are the stress equilibrium with the energetic and dissipative stresses

σe = ∂εψ̂(c) and σd ∈ ∂ε̇φ̂(ċ) (3.83)

and Biot’s equation that determines the evolutions of the internal variables {αi}i=1,n.

Often, the dissipative effects are modeled via a dual dissipation potential function that
depends on dissipative driving forces

f(t) = ( σd(t), β1(t), . . . , βn(t) ) (3.84)

which are conjugate to the state c defined in (3.76). In order to get a variational frame-
work which includes these forces, we express the dissipation potential by its dual via the
Legendre transformation

φ̂(ċ) = sup
f

[ f · ċ− φ̂∗(f) ] . (3.85)

This motivates the definition of the mixed rate-type potential

π̂∗(ċ, f; ct, t) =
d

dt
ψ̂(c) + f · ċ− φ̂∗(f)− σext(t)ε̇ (3.86)

for given state at time t. Then, the mixed saddle point principle

{ ċt, ft } ∈ Arg{ inf
ċ

sup
f

π̂∗(ċ, f; ct, t) } (3.87)

determines at time t the rate of the constitutive state as well as the dissipative driving
forces. The Euler equations of this principle are

1. ∂ε̇π̂
∗(ċt, ft; ct, t) ≡ ∂εψ̂(ct) + σd

t − σext(t) = 0

2. ∂α̇i
π̂∗(ċt, ft; ct, t) ≡ ∂αi

ψ̂(ct) + (βi)t = 0

3. ∂σd π̂∗(ċt, ft; ct, t) ≡ ε̇t − ∂σd φ̂∗(ft) ∋ 0

4. ∂βi
π̂∗(ċt, ft; ct, t) ≡ (α̇i)t − ∂βi

φ̂∗(ft) ∋ 0 .

(3.88)

The two differential inclusions are evolution equations for the internal variables or can
alternatively be seen as definitions of the dissipative driving forces in an inverse format.

3.3. Modeling Examples

We now apply the variational framework discussed above to four rheological models
that consist of serial and parallel arrangements of the fundamental spring, dashpot and
friction devices.


36 3 Variational Principles for Rheological Models

ε

α

σext σext
E

H1

H2

1©

2©3©

Figure 8: Visco-elastic material model. An outer dashpot labeled by 1© is parallely con-
nected to a device that consists of a serial arrangement of an inner dashpot labeled by 2©
and an inner spring labeled by 3©.

3.3.1. Visco-Elasticity. As a first example, we consider the material element de-
picted in Figure 8 which describes a (smooth) visco-elastic behavior. Following the recipe
in Section 3.2, we first identify the constitutive state array as

c = (ε, α) (3.89)

which consists of the external strain and the strain associated with the inner dashpot.
The strains in all three basic devices follow as

ε1 = ε , ε2 = α and ε3 = ε− α . (3.90)

Next, we define the stored-energy function associated with the spring device

ψ̂(c) = ψ̂3(ε3) =
1

2
E(ε− α)2 (3.91)

and the dissipation potential function associated with the two dashpot devices

φ̂(ċ) = φ̂1(ε̇1) + φ̂2(ε̇2) =
1

2
H1ε̇

2 +
1

2
H2α̇

2 . (3.92)

These functions, together with the external power term, govern the rate-type potential
(3.80) taking the specific form

π̂(ċ; ct, t) = ∂εψ̂(ct) ε̇+ ∂αψ̂(ct) α̇+
1

2
H1ε̇

2 +
1

2
H2α̇

2 − σext(t)ε̇ (3.93)

for given state at time t. Then, the rates of the external strain and the strain of the inner
dahspot at time t are obtained from the minimization principle

{ ε̇t, α̇t } = Arg{ inf
ε̇,α̇

π̂(ċ; ct, t) } . (3.94)

We have the corresponding Euler equations

1. ∂ε̇π̂(ċt; ct, t) ≡ E(εt − αt) +H1ε̇t − σext(t) = 0

2. ∂α̇π̂(ċt; ct, t) ≡ −E(εt − αt) +H2α̇t = 0
(3.95)


3.3 Modeling Examples 37

ε

α

σextσext

E y0
1© 2©

Figure 9: Perfect elasto-plastic material model. A spring labeled by 1© is connected in
series to a friction device labeled by 2©. The strain associated with the spring is called
elastic strain and is simply ε1 = ε− α.

which are specific forms of (3.82).

As an alternative to this primal formulation, we consider the dual dissipation potential
function associated with the two dashpot devices

φ̂∗(σd, β) =
1

2H1
(σd)2 +

1

2H2
β2 , (3.96)

where the dissipative driving forces f = (σd, β) are conjugate to the strain variables
c = (ε, α). The mixed rate-type potential (3.86) takes the specific form

π̂∗(ċ, f; ct, t) = ∂εψ̂(ct) ε̇+ ∂αψ̂(ct) α̇ + σdε̇+ βα̇−

[
1

2H1

(σd)2 +
1

2H2

β2

]
− σext(t)ε̇ .

(3.97)

Then, the rates of the external strain and the strain of the inner dashpot as well as the
dissipative driving forces at time t are governed by the mixed saddle point principle

{ ε̇t, α̇t, σ
d
t , βt } = Arg{ inf

ε̇,α̇
sup
σd,β

π̂∗(ċ, f; ct, t) } . (3.98)

The Euler equations are

1. ∂ε̇π̂
∗(ċt, ft; ct, t) ≡ E(εt − αt) + σd

t − σext(t) = 0

2. ∂α̇π̂
∗(ċt, ft; ct, t) ≡ −E(εt − αt) + βt = 0

3. ∂σd π̂∗(ċt, ft; ct, t) ≡ ε̇t −
1
H1
σd
t = 0

4. ∂β π̂
∗(ċt, ft; ct, t) ≡ α̇t −

1
H2
βt = 0

(3.99)

and are specific forms of (3.88). In the first relationship, the externally applied stress is
equilibrated with the stresses in the spring and the outer dashpot. The second equation
identifies the dissipative driving force βt with the stress in the spring. The last two
equations are evolution laws for the external strain as well as the strain associated with
the inner dashpot.

3.3.2. Perfect Elasto-Plasticity. As a second example, we take a look at the stan-
dard model of perfect elasto-plasticity that consists of a serial arrangement of a spring
and friction device, see Figure 9. The material response is nonsmooth. As before, we
specify the constitutive state

c = (ε, α) (3.100)


38 3 Variational Principles for Rheological Models

which contains the external strain and the strain associated with the friction device.
Latter will also be referred to as plastic strain. Then, the strains in the single devices are

ε1 = ε− α and ε2 = α , (3.101)

where ε1 is called elastic strain. The stored-energy function of the spring device is

ψ̂(c) = ψ̂1(ε1) =
1

2
E(ε− α)2 (3.102)

and the dissipation potential function associated with the friction element takes the form

φ̂(ċ) = φ̂2(ε̇2) = y0|α̇| , (3.103)

see Section 3.1.3. With these two functions at hand, the rate-type potential (3.80) for a
given state at time t is specified to

π̂(ċ; ct, t) = ∂εψ̂(ct) ε̇+ ∂αψ̂(ct) α̇+ y0|α̇| − σext(t)ε̇ . (3.104)

Note that stress equilibrium

σt = σext(t) with σt = ∂εψ̂(ct) = E(εt − αt) (3.105)

is a necessary condition on the given state. Insertion into (3.104) yields the reduced
rate-type potential

π̂red(α̇; ct) = ∂αψ̂(ct) α̇+ y0|α̇| . (3.106)

Then, the rate of the plastic strain at time t follows from the minimization principle

α̇t ∈ Arg{ inf
α̇
π̂red(α̇; ct) } (3.107)

with the corresponding Euler equation

∂α̇π̂red(α̇t; ct) ≡ −E(εt − αt) + y0





−1 for α̇t < 0
[−1, 1] for α̇t = 0
1 for α̇t > 0

∋ 0 . (3.108)

Recalling the constitutive defintion (3.105)2 of the stress, the solution of (3.107) is α̇t = 0
for |σt| < y0, but is not unique for |σt| = y0, i.e. α̇t ∈ R≥0 for σt = y0 and α̇t ∈ R≤0

for σt = −y0, see also Section 3.1.3. The rate of the total strain at time t follows from
the rate form of stress equilibrium (3.105), i.e. ε̇t = σ̇ext(t)/E for α̇t = 0 and ε̇t = α̇t for
α̇t 6= 0 which means that the elastic strain stays constant during plastic evolution.

We now consider the dissipative force array f = (0, β) conjugate to the constitutive
state c = (ε, α). To construct a mixed formulation, we define the dissipation potential
function via the principle of maximum plastic dissipation

φ̂(α̇) = sup
β∈E

[ βα̇ ] (3.109)

in terms of the elastic domain which is defined in the space of dissipative driving forces

E = { β | χ̂(β) ≤ y0 } with χ̂(β) = |β| . (3.110)


3.3 Modeling Examples 39

PSfrag

ε

α

E1

E3

y0
σextσext

1©

2©

3©

Figure 10: Elasto-plastic material model with kinematic hardening. A spring labeled 1© is
connected in series to a parallel arrangement of a friction device labeled by 2© and another
spring labeled by 3©.

We define for given constitutive state at time t the reduced mixed rate-type potential

π̂∗
λ,red(α̇, β, λ; ct) = ∂αψ̂(ct) α̇ + βα̇− λ (|β| − y0) (3.111)

where the last term should take into account the constraint stemming from the elastic
domain (3.110). Then, the rate of the plastic strain, the dissipative driving force and the
Lagrange multiplier at time t are governed by the mixed saddle point principle

{ α̇t, βt, λt } ∈ Arg{ inf
α̇

sup
β

inf
λ≥0

π̂∗
λ,red(α̇, β, λ; ct) } . (3.112)

The corresponding Euler equations read

1. ∂α̇π̂
∗
λ,red(α̇t, βt, λt; ct) ≡ −E(εt − αt) + βt = 0

2. ∂β π̂
∗
λ,red(α̇t, βt, λt; ct) ≡ α̇t − λt sgn βt = 0

3. ∂λπ̂
∗
λ,red(α̇t, βt, λt; ct) ≡ −|βt|+ y0 ≥ 0, λt ≥ 0, λt(|βt| − y0) = 0

(3.113)

and we identify the dissipative driving force βt associated with the friction device with
the stress (3.105) acting in the spring. In addition, (3.113) contains the evolution law
for the plastic strain together with the classical loading/unloading conditions in Karush-
Kuhn-Tucker form. If |βt| < y0, then λt = 0 and we have no plastic evolution α̇t = 0. On
the other hand, for βt = y0 and βt = −y0 we get the nonunique solutions α̇t = λt ≥ 0 and
α̇t = −λt ≤ 0, respectively. The plastic consistency condition (3.73) specifies to β̇t = 0
at stage of plastic flow and yields a relationship between the Lagrange multiplier and the
total strain rate, i.e. λt = ε̇t sgn σext(t).

3.3.3. Elasto-Plasticity with Kinematic Hardening. As a third example, we
add kinematic hardening to the perfect elasto-plasticity model discussed in Section 3.3.2.
This is achieved by inserting another spring parallel to the friction device, see Figure 10.
The constitutive state c = (ε, α) still contains the external and the plastic strains. We
now have to take into account the strain ε3 of the added spring

ε1 = ε− α , ε2 = α and ε3 = α . (3.114)

The stored-energy function associated with the two springs reads

ψ̂(c) = ψ̂1(ε1) + ψ̂3(ε3) =
1

2
E1(ε− α)2 +

1

2
E3α

2 (3.115)


40 3 Variational Principles for Rheological Models

and the dissipation potential function associated with the friction device remains un-
changed