Structural optimisation of diffusion driven degradation processes

Abstract

In this article, we propose an optimisation framework that can contribute to the prevention of premature failure or damage to building structures and can thereby strengthen their longevity. We concentrate on structures that are contaminated by chemical substances and that have strong destructive effects on the material. The aim of this contribution is a mathematical algorithm that allows the optimisation of a structure exposed to chemical influences and thus the assurance of the static load capacity. To achieve this, we present a coupled mechanical-diffusion-degradation approach embedded in a finite element (FE) framework. Furthermore, we integrate an optimisation algorithm to reduce material degradation. In this paper, we establish shape optimisation of a structure with a gradient based optimisation algorithm.

1 Introduction

Engineering structures are dimensioned according to standards. The maximum existing stresses are evaluated and the material load-bearing capacity is examined, whereby environmental influences are only marginally considered. However, negative influences can additionally change the material composition and thus the mechanical load-bearing capacity over time. Most changes in the internal structure of materials are associated with diffusion processes. One example is the long-term effect of calcium leaching in concrete, where “pure water creates concentration gradients which lead to the diffusion of Ca ions from the pore water and the subsequent degradation of underground concrete” (Choi and Yang 2013).

In the following, an approach is presented which allows the calculation of mechanical and chemical influences on structures. Furthermore, it is possible to calculate an optimal geometry that reduces the damage of long-term effects caused by environmental influences. Therefore, this approach integrates the structural analysis into a structural optimisation algorithm. Within the structural analysis, we couple the diffusion of chemical substances with the mechanical behaviour of the material. We assume permeable structures which allow gradient-based diffusion of concentrations. In addition, we postulate that chemical substances have an influence on the material of the structure, since they can trigger material degradation. To ensure this, we focus on mechanical degradation processes corresponding to negative growth processes using multiplicative decomposition of the deformation gradient. Since this is a space- and time-dependent problem, we embed the material description in a finite element (FE) framework, by using an isoparametric concept for space discretisation and a Newmark-beta approach for temporal discretisation. Furthermore, we develop an algorithm for structural optimisation. In detail, we present a shape optimisation with a gradient-based calculation, which contains information about the sensitivities of the parameters that influence the mechanical behaviour. Using representative examples, we introduce the basic idea of the model and examine practical aspects where components are at risk of being exposed to chemical concentrations and yet their strength must be guaranteed.

In order to determine the approach for optimisation of diffusion driven degradation problems, the following topics are applied and embedded in a short overview of the state of the art. In the context of structural mechanics, the topics coupled problems and growth, respectively degradation processes, are investigated in more detail. For the embedding in optimisation, shape optimisation and sensitivity analysis are applied.

Kuhl (2005) presents a general overview of coupled problems and the numerical implementation. Furthermore, he outlines in detail the coupling of diffusion processes to mechanical behaviour. A coupled theory of fluid permeation and large deformations for elastomeric materials, concentrating on a thermodynamically consistent derivation of the constitutive relations and the resulting partial differential equations, is presented by Chester and Anand (2010).

Menzel and Kuhl (2012) summarise growth and remodelling models for living structures. Mechanical growth and remodelling can be modelled either with a constitutive approach, a kinematic approach or a combination of both. Growth processes can be described by evaluating the time-dependent change in mass, density or volume of a structure. On the one hand, the constitutive approach concentrates on a thermodynamically consistent evaluation of the mass source and the mass flow, which enables the calculation of change in mass or density (Cowin and Hegedus 1976; Harrigan and Hamilton 1992; Epstein and Maugin 2000; Kuhl et al. 2003; Menzel 2005). The kinematic approach, on the other hand, allows the calculation of the variable mass or volume by applying a multiplicative decomposition of the deformation gradient (Rodriguez et al. 1994; Chen and Hoger 2000; Ambrosi and Mollica 2002; Gleason and Humphrey 2005; Menzel 2007). A useful alternative is the combination of the constitutive and kinematic approach, one example for this combination is presented by (Ganghoffer 2010), where the calculation of surface growth in biological tissue is presented. He uses the multiplicative decomposition of the deformation gradient, as originally presented by Rodriguez et al. (1994), and applies a thermodynamically consistent approach to establish an evolutionary law for growth velocity. The multiplicative decomposition of the deformation gradient leads to a “growth tensor describing the local addition of material and an elastic tensor characterizing the reorganization of the body”, see Ganghoffer and Plotnikov (2014). Therefore, the development of the growth tensor, the so-called transplant tensor, is introduced as a state variable in the framework of finite elasticity.

In Gérard et al. (2002), a simplified numerical model of calcium leaching is presented, which concentrates on the chemical state variable calcium and the kinetics of the leaching process. This model considers time-dependent chemical processes and mentions the possibility of a simple coupling in a finite element (FE) algorithm. Kuhl (2005) presents a chemo-mechanical model for the numerical simulation of calcium leaching in concrete. He emphasises the effects of the chemical degradation on the pore volume and the mechanical stability of concrete. A damage function is used for constitutive modelling of the chemo-mechanically degraded material.

An overview on structural optimisation methods applied to discretised linear-elastic structures is given in Christensen and Klarbring (2008). Structural optimisation, which is applied to growth or degradation processes, is usually solved by evolutionary optimisation algorithms (ESO), see Simon (2013) for details. The evolution models are based on the concept of gradually removing inefficient material from a structure. The disadvantage of ESO concepts is that optimisation starts from one defined reference configuration and that the result quickly ends in local minima. Thus, it is not possible to find an optimal solution with a single ESO algorithm, see Vrugt and Robinson (2007). However, the application of shape optimisation to growth or degradation processes could help solve the problem, but has not yet been sufficiently investigated. Similar to the embedding of a growth approach in an optimisation algorithm, shape and topology optimisations are applied to damage models. For example, an optimal damage distribution is computed with a shape optimisation algorithm that embeds an isotropic gradient–enhanced damage model, see Guhr et al. (2020). Suresh et al. (2018) present a topology optimisation method which is connected to a fatigue model. The model enables the computation of optimised topologies, taking fatigue at high cycles as a limiting condition into account. Moreover, Noël et al. (2017) develop a level set–based topology optimisation framework to reduce damage in the context of structural design.

Barthold (2002) presents the basic principles of design changes necessary for structural optimisation and their effects on the structural response on a continuum. He introduces the local convective approach containing local coordinates and derives it from a differentiable manifold. Furthermore, he emphasises the importance of this approach for obtaining information about the kinematic relation required for numerical methods such as the finite element method (FEM) and computer-aided geometric design (CAGD). A detailed discussion of efficient strategies for calculating sensitivity of design parameters is described in Barthold and Stein (1996) where, in addition, the method of variational sensitivity analysis is presented.

The above briefly presented literature motivates the structural analysis which contains a coupled mechanical-diffusion-degradation approach. The special feature of this model is the combination of the coupled problem with a shape optimisation approach so that the failure limit of structures affected by mechanical loads and chemical impacts can be improved. By embedding the degradation model in the optimisation algorithm, this approach represents an effective alternative to evolutionary algorithms.

2 Continuum model

We present the continuum model including coupling between diffusion and mechanical behaviour. First we introduce the extended kinematic framework. We then formulate the partial differential equations which are determined by balance equations and constitutive equations together with boundary and initial conditions. The constitutive equations, including mechanical stresses, particle diffusion and chemical reactions, are determined by a thermodynamically consistent framework based on the evaluation of the entropy inequality.

2.1 Kinematics

A graphical illustration of the applied kinematic framework is presented in Fig. 1, which provides a multiplicative decomposition of the deformation gradient into an elastic and a growth part. We assume that the temporal evolution of the mass density takes place in the reference configuration. An infinitesimal volume element is represented by dV in the reference configuration, by dv_d in the degradation space and by dv in the actual configuration. We introduce positions in the reference configuration X, in the actual configuration x and in the degradation space x_d. Furthermore, φ_t describes the mapping of the reference particles X onto their spatial position x = φ_t(X,t) with t representing time. We introduce the convective tangent vectors G_i,g_i and h_i. The contravariant basis vectors that determine the dual basis are defined by $\textbf {G}_{\text {i}}\cdot \textbf {G}^{\text {j}}=\delta ^{\text {j}}_{\text {i}}$, $\textbf {G}_{\text {i}}\cdot \textbf {G}^{\text {j}}=\delta ^{\text {j}}_{\text {i}}$ and $\textbf {h}_{\text {i}}\cdot \textbf {h}^{\text {j}}=\delta ^{\text {j}}_{\text {i}}$. With this, we introduce the deformation tensors which map infinitesimal line elements represented in different configurations, i.e. F, F^d and F^e

$$ \begin{array}{lllc} \text{deformation gradient} : & \textbf{F} &=& {\textbf{g}}_{\text{i}} \otimes {\textbf{G}}^{\text{i}}\\ \text{degradation gradient} : & \textbf{F}^{\text{ d}} &=& \textbf{h}_{\text{i}} \otimes {\textbf{G}}^{\text{i}}\\ \text{elastic deformation gradient} : & \textbf{F}^{\text{e}} &=& {\textbf{g}}_{\text{i}} \otimes \textbf{h}^{\text{i}} . \end{array} $$

(1)

Fig. 1

Graphical illustration of the kinematic concept

These gradients are important two-point tensors that allow transformations between objects in relation to the respective configurations. We apply a multiplicative decomposition of the deformation gradient F into an elastic part F^e and a degradation part F^d as introduced in Himpel et al. (2005), Kuhl et al. (2004), and Menzel and Kuhl (2012). According to Lubarda and Hoger (2002), we use an isotropic approach for the degradation fraction of the deformation gradient with the stretch ratio ν, to be specific

$$ \begin{array}{@{}rcl@{}} \begin{array}{lllcl} \textbf{F} &=& \textbf{F}^{\mathrm{e}} \textbf{F}^{\mathrm{d}} \qquad \qquad &\text{with}& \qquad \\ \textbf{F}^{\mathrm{d}} &=& \nu \textbf{1} & \text{and}& \qquad \nu = \sqrt[3]{\frac{\rho_{0}}{\rho_{0}^{*}}} , \end{array} \end{array} $$

(2)

whereby 1 denotes the second order identity tensor and the mass densities ρ₀ and $\rho _{0}^{*}$ respectively, being introduced as this work proceeds. Using the elastic part of the deformation gradient, we can introduce the elastic right Cauchy Green tensor as

$$ \begin{array}{llll} \textbf{C}^{\mathrm{e}} = (\textbf{F}^{\text{e}})^{\mathrm{T}} \textbf{F}^{\text{e}} = {\text{g}}_{\text{ij}} \textbf{h}^{\text{i}} \otimes \textbf{h}^{\text{j}} . \end{array} $$

(3)

Moreover, the material time derivative of the elastic right Cauchy Green tensor is obtained by

$$ \begin{array}{llll} \dot{\textbf{C}}^{\mathrm{e}}&=& (\dot{\textbf{F}}^{\text{e}})^{\mathrm{T}} \textbf{F}^{\text{e}} +(\textbf{F}^{\text{e}})^{\mathrm{T}} \dot{\textbf{F}}^{\text{e}} \\ &=& (\textbf{F}^{\mathrm{e}})^{\mathrm{T}} \textbf{l}^{\mathrm{T}} \ \textbf{F}^{\mathrm{e}} + (\textbf{F}^{\mathrm{e}})^{\mathrm{T}} \textbf{l} \ \textbf{F}^{\mathrm{e}} \\ &=& 2 (\textbf{F}^{\mathrm{e}})^{\mathrm{T}}\textbf{d} \textbf{F}^{\mathrm{e}} , \end{array} $$

(4)

with

$$ \dot{\textbf{F}}^{\mathrm{e}}= \textbf{l} \textbf{F}^{\mathrm{e}} \quad \text{and} \quad \textbf{d}= \frac{1}{2} (\textbf{l}^{\mathrm{T}} +\textbf{l}) , $$

(5)

wherein $\textbf {l} = \text {grad} \dot {\textbf {x}}$ is the spatial velocity gradient.

Figure 1 introduces the definition of the mass densities referred to the respective configurations. In the reference configuration, $\rho _{0}^{*}$ represents the initial and ρ₀ the referential mass density, which arises as a result of material degradation. The mass density ρ_d is referred to the degradation space and ρ_t to the actual configuration. Mass exchange is realised by a mass sink term R₀ per unit volume in the reference configuration which triggers the degradation from the initial mass contribution dM and results in the mass contribution dm, i.e.

$$ \begin{array}{@{}rcl@{}} \text{dm} = \text{dM} + \int\limits_{\text{t}_{0}}^{\mathrm{t}} \text{R}_{0} \text{dt} \text{dV} . \end{array} $$

(6)

This leads to time-dependent update of the initial mass density $\rho _{0}^{*}$ to the resulting referential mass density ρ₀, to be specific

$$ \begin{array}{@{}rcl@{}} \rho_{0} = \rho_{0}^{*} + \int\limits_{\text{t}_{0}}^{\mathrm{t}} \text{R}_{0} \text{dt} . \end{array} $$

(7)

Within the kinematic concept, we assume that the initial referential mass density and the mass density referred to the degradation space coincide, i.e.

$$ \begin{array}{@{}rcl@{}} \rho_{0}^{*} = \rho_{\mathrm{d}} . \end{array} $$

(8)

The mapping between the configurations, see Fig. 1, leads to a change in volume, whereas mass is not influenced by such mappings, i.e. using (8) results in

$$ \begin{array}{@{}rcl@{}} \rho_{0} = \rho_{0}^{*} \text{J}^{\mathrm{d}} = \rho_{\mathrm{d}} \text{J}^{\mathrm{d}} = \rho_{\mathrm{t}} \text{J} ; \rho_{\mathrm{d}} = \rho_{\mathrm{t}} \text{J}^{\mathrm{e}} . \end{array} $$

(9)

Moreover, $\text {J}^{\mathrm {e}}=\det \textbf {F}^{\mathrm {e}}>0$ denotes the determinant of the elastic part of the deformation gradient and $\text {J}=\det \textbf {F}>\text {0}$ as well as $\text {J}^{\mathrm {d}}=\det \textbf {F}^{\mathrm {d}}=\rho _{0} / \rho _{\mathrm {0}}^{*}>0$, cf. (2), follow by analogy.

2.2 Degradation approach

The degradation model is implemented by a combination of a kinematic and a constitutive approach. Thus, in Section 2.1, the multiplicative decomposition of the deformation gradient is performed and the degradation space is introduced. Within the constitutive approach, a mass exchange is established which leads to the change of referential mass density ρ₀. In return, the evolving mass density is used to calculate the degradation part of the deformation gradient, see (2). We assume that the mass degradation is caused due to chemical concentrations c_γ, which are defined by the corresponding molar density ρ_γ and the molar mass M_γ, i.e.

$$ \text{c}_{\gamma} = \frac{\rho_{\gamma}}{\text{M}_{\gamma}} . $$

(10)

Thus, we apply the constitutive approach of the mass sink term as follows

$$ \begin{array}{cll} \text{R}_{0} = -\dot{\rho}_{\gamma} . \end{array} $$

(11)

Assuming only one concentration to be related to the chemical interaction, we can summarise the degradation impact of the concentrations as

$$ \begin{array}{@{}rcl@{}} \rho_{0} = \rho_{0}^{*} - \text{c}_{\gamma} \text{M}_{\gamma} . \end{array} $$

(12)

This paper does not address chemical processes in detail. For further details on a model with more complex chemical processes, the reader is referred to, e.g. Gérard et al. (2002). In the following, the concentrations are introduced as an additional degree of freedom in the continuum.

2.3 Balance equations

In order to implement the above-introduced problem numerically, the following section presents the necessary balance approaches to solve the coupled multi-field problem. On the one hand, we need the balance of mass to calculate the concentrations, and on the other hand, the momentum balance to calculate the displacement. The concentration distribution depends on spatial conditions and the stresses on mass density which, in turn, is influenced by concentration-related mass reduction. For this reason, this is a strongly coupled problem. Furthermore, the balance of energy and entropy is needed to determine the constitutive material model.

2.3.1 Balance of mass

We evaluate the balance of mass for the macroscopic body and for the chemical concentrations. The reduction of the mass in the reference configuration leads to a change of mass with $\text {d}\dot {\mathrm {m}}= \text {R}_{0} \text {dV}$ and $\dot {\{{\bullet }\}}$ denoting the material time derivative. Using the transformations between the configurations, we can represent the balance of mass on the actual configuration as

$$ \begin{array}{cll} \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \left( \dot{\rho}_{\text{t}} + \rho_{\text{t}} \text{div} \dot{\textbf{x}}\right) \text{dv} = {\int\limits_{\text{B}_{\text{t}}}} \frac{\text{R}_{0}}{\text{J}} \text{dv} , \end{array} $$

(13)

wherein ‘div’ denotes the spatial divergence operator. The balance of mass of the chemical concentrations contains the material time derivative of the concentrations $\dot {\text {c}}_{\gamma }$ and the flux of the concentrations j_γ, i.e.

$$ \begin{array}{cll} \displaystyle{\int\limits_{\text{B}_{\text{t}}}} (\dot{\text{c}}_{\gamma} + \text{div} \textbf{j}_{\gamma}) \text{dv} = 0 . \end{array} $$

(14)

2.3.2 Balance of linear momentum

The balance of linear momentum is shown below for the actual configuration and contains the Cauchy stress tensor T. Although we apply mass exchange, we can neglect effects of the impulse resulting from the degrading mass, since the process of degradation is assumed to be very slow. This results in

$$ \begin{array}{cll} \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \text{div} \textbf{T} \text{dv} = \textbf{0} , \end{array} $$

(15)

wherein additional volume forces as well as acceleration are neglected. From the balance equation of moment of momentum, we can derive that the Cauchy stress tensor is symmetric, so that T^T = T.

2.3.3 Balance of energy

The balance of energy consists of the temporal derivative of internal energy $\dot {\text {E}}$, and kinetic energy $\dot {\text {K}}$, which correspond to mechanical energy change dW, thermal energy change dQ and chemical flow on the surface $\text {E}_{\text {c}_{\gamma }}$, resulting in the following equations

$$ \begin{array}{lllllll} \dot{\text{E}} + \dot{\text{K}}&=& \text{dW} + \text{dQ} + \text{E}_{\text{c}_{\gamma}} \\ \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \dot{\text{e}} \text{dv} &=& \displaystyle{\int\limits_{\text{B}_{\text{t}}}} (\textbf{d} : \textbf{T} + \text{r} - \text{div} \textbf{q}) \text{dv}\\ &&- \displaystyle{\int\limits_{\partial\text{B}_{\text{t}}}} (\mu_{\gamma} \textbf{j}_{\gamma}) \cdot \text{d}\textbf{a} , \end{array} $$

(16)

wherein $\dot {\text {e}}$ is the material time derivative of the volume specific internal energy, r is the volume specific heat source, q is the heat flux density, μ_γ is the chemical potential and da = n da with spatial outward unit surface vector n. The time derivation of the kinetic energy has no impact on the balance of energy as we neglect acceleration and consider only slow progression of degradation. The surface part can be reformulated as follows,

$$ \begin{array}{llll} \displaystyle{\int\limits_{\partial\text{B}_{\text{t}}}} (\mu_{\gamma} \textbf{j}_{\gamma}) \cdot \text{d}\textbf{a} &=& \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \text{div} (\mu_{\gamma} \textbf{j}_{\gamma}) \text{dv} \end{array} $$

(17)

with

$$ \text{div} (\mu_{\gamma} \textbf{j}_{\gamma}) = \textbf{j}_{\gamma} \cdot \text{grad} \mu_{\gamma} + \mu_{\gamma} \text{div} \textbf{j}_{\gamma} .\\ $$

(18)

Using the mass balance of chemical concentrations from (14), the local form of the balance of energy reads

$$ -\text{r} + \text{div} \textbf{q} = - \dot{\mathrm{e}} + \textbf{d} : \textbf{T} - \textbf{j}_{\gamma} \cdot \text{grad} \mu_{\gamma} + \mu_{\gamma} \dot{\text{c}}_{\gamma} . $$

(19)

2.3.4 Balance of entropy

Entropy is summarised in the second law of thermodynamics and states that entropy never decreases in a closed system. Accordingly, entropy always increases or remains constant. In this context, the entropy inequality is adapted as

$$ \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \dot{\text{s}} \text{dv} \geq \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \frac{1}{\Theta} \text{r} \text{dv} - \displaystyle{\int\limits_{\partial\text{B}_{\text{t}}}} \frac{1}{\Theta} \textbf{q}\cdot \text{d}\textbf{a} , $$

(20)

where $\dot {\text {s}}$ denotes the material time derivative of the volume-specific entropy and where Θ is the absolute temperature. The surface part can be reformulated as follows,

$$ \begin{array}{llll} \displaystyle{\int\limits_{\partial\text{B}_{\text{t}}}} \frac{1}{\Theta} \textbf{q} \cdot \text{d}\textbf{a} &=& \displaystyle{\int\limits_{\text{B}_{\text{t}}}} \text{div} \left( \frac{1}{\Theta} \textbf{q}\right) \text{dv} \end{array} $$

(21)

with

$$ \begin{array}{llll} \text{div} \left( \frac{1}{\Theta} \textbf{q}\right) &=& \textbf{q} \cdot \text{grad}\left( \frac{1}{\Theta}\right) + \frac{1}{\Theta} \text{div} \textbf{q} \\ &=& - \frac{1}{{\Theta}^{2}} \textbf{q} \cdot \text{grad} {\Theta} + \frac{1}{\Theta} \text{div} \textbf{q} . \end{array} $$

(22)

Inserting the reformulation into the entropy inequality (20), we obtain the local form

$$ 0 \leq {\Theta} \dot{\text{s}} - \text{r} - \frac{1}{\Theta} \textbf{q} \cdot \text{grad} {\Theta} + \text{div} \textbf{q} . $$

(23)

Under isothermal conditions, i.e. q = 0 as well as $\dot {\Theta }=0$, and with the local form of the energy balance, see (19), we end up with the restriction

$$ 0 \leq {\Theta} \dot{\text{s}} - \dot{\mathrm{e}} + \textbf{d} : \textbf{T} \\ - \textbf{j}_{\gamma} \cdot \text{grad} \mu_{\gamma} + \mu_{\gamma} \dot{\text{c}}_{\gamma} . $$

(24)

2.4 Constitutive formulation

We perform a Legendre(-Fenchel) transformation between the extensive thermodynamic entropy s and the conjugate quantity, the temperature Θ. This results in the so-called Helmholtz energy ψ, i.e.

$$ \begin{array}{lllllll} \dot{\psi} = \dot{\mathrm{e}} - {\Theta} \dot{\text{s}} . \end{array} $$

(25)

By introducing the Helmholtz free energy into the entropy inequality (24) and by selecting the elastic right Cauchy Green tensor and the concentrations to derive the energy process with $\psi (\textbf {C}^{\text {e}},\text {c}_{\gamma })$, the following results are obtained

$$ \begin{array}{cccc} \underbrace{- \frac{\partial \psi}{\partial \textbf{C}^{\text{e}}} : \dot{\textbf{C}}^{\text{e}} + \textbf{d} : \textbf{T}}_{\text{reversible}} &+& \underbrace{ - \frac{\partial \psi}{\partial \text{c}_{\gamma}} \dot{\text{c}}_{\gamma} + \mu_{\gamma} \dot{\text{c}}_{\gamma} }_{\text{reversible}} \\ \underbrace{- \textbf{j}_{\gamma} \cdot \text{grad} \mu_{\gamma}}_{\text{ireversible}} & & \geq 0 . \end{array} $$

(26)

We can separate the approach for the Helmholtz function into a mechanical ψ^M and a chemical ψ^C part

$$ \psi(\textbf{C}^{\text{e}},\text{c}_{\gamma})= \psi^{\text{M}}(\textbf{C}^{\text{e}})+\psi^{\text{C}}(\text{c}_{\gamma}) . $$

(27)

In order to ensure the entropy inequality, the reversible contributions are evaluated and the Cauchy stress T and the chemical potential μ_γ are determined, i.e.

$$ \begin{array}{lll} \textbf{T} &=& 2 \textbf{F}^{\text{e}} \frac{\partial \psi^{\text{M}}}{\partial\textbf{C}^{\text{e}}} (\textbf{F}^{\text{e}})^{\mathrm{T}} \\ \mu_{\gamma} &=& \frac{\partial \psi_{\gamma}^{\text{C}}}{\partial \text{c}_{\gamma}} , \end{array} $$

(28)

wherein (4) is used. From the irreversible part of the entropy inequality we motivate the flux of the concentrations

$$ \begin{array}{lll} \textbf{j}_{\gamma} &=& -\text{D} \text{grad} \text{c}_{\gamma} , \end{array} $$

(29)

where D ≥ 0 is introduced as the diffusion coefficient.

2.4.1 Specifications of energy contributions

The mechanical energy is described by adopting a hyperelastic Neo-Hooke material, ψ^Neo, which contains material parameters μ and λ, to be specific

$$ \begin{array}{lll} \psi^{\text{M}} &=& \rho_{\text{t}} \psi^{\text{Neo}} \\ \psi^{\text{Neo}} &=&\frac{1}{\rho_{0}*} \left[\frac{1}{2} \lambda (\text{J}^{\mathrm{e}} - 1)^{2} - \mu \ln\text{J}^{\mathrm{e}} \right.\\ & & \left.+ \frac{1}{2} \mu (\text{I}_{\mathrm{C}}^{\text{e}} - 3)\right] , \end{array} $$

(30)

with invariants $\text {I}_{\mathrm {C}^{\text {e}}} = \text {tr} \textbf {C}^{\text {e}}$ and $\text {J}^{\text {e}} = \sqrt {\det \textbf {C}^{\text {e}}} = \sqrt {\text {III}_{\mathrm {C}^{\text {e}}}}$, which depend on the elastic right Cauchy Green deformation tensor C^e. In addition, we specify the chemical part ψ^C as

$$ \psi^{\text{C}} = \text{c}_{\gamma} \mu_{\gamma}^{0} + \text{R} {\Theta} (-\text{c}_{\gamma} + \text{c}_{\gamma} \ln \text{c}_{\gamma}) , $$

(31)

where $\mu _{\gamma }^{0}$ is the constant standard potential and where R is the gas constant. Inserting the mechanical and chemical energy contributions into (28) results in

$$ \begin{array}{lll} \textbf{T} &=& \frac{\rho_{\text{t}}}{\rho_{0}*} \textbf{F}^{\text{e}} \left[\lambda (\text{J}^{\mathrm{e}} -1)\text{J}^{\mathrm{e}} (\textbf{C}^{\text{e}})^{\mathrm{-1}} \right.\\ && -\left. \mu (\textbf{C}^{\text{e}})^{\mathrm{-1}} + \mu \textbf{1}\right](\textbf{F}^{\text{e}})^{\text{T}} \\ &=& \frac{\rho_{\text{t}}}{\rho_{0}*} \left[\lambda (\text{J}^{\mathrm{e}} -1)\text{J}^{\mathrm{e}} \textbf{1} + \text{2} \mu \textbf{K}^{\mathrm{e}}\right] , \end{array} $$

(32)

with the spatial Karni-Reiner strain tensor K^e, i.e.

$$ \textbf{K}^{\mathrm{e}}=\frac{1}{2} [\textbf{F}^{\text{e}}(\textbf{F}^{\text{e}})^{\text{ T}} - \textbf{1}] . $$

(33)

Furthermore, the chemical potential follows as

$$ \mu_{\gamma} = \mu_{\gamma}^{0} + \text{R} {\Theta} \ln \text{c}_{\gamma} . $$

(34)

2.4.2 Illustration of stress state

In the evaluations of the examples in Sections 4 and 6 we use the principal Cauchy stresses T₁ and T₂ in the two-dimensional space considered, namely

$$ \text{T}_{1,2} = 1 \pm \sqrt{- \frac{\text{I}_{\mathrm{T}}}{2} - \text{II}_{\mathrm{T}}} , $$

(35)

with the invariants of the Cauchy stress I_T and II_T, i.e.

$$ \begin{array}{lll} \text{I}_{\mathrm{T}} &=& \text{tr} \textbf{T} \\ \text{II}_{\mathrm{T}} &=& \frac{1}{2} ((\text{tr} \textbf{T})^{\mathrm{2}} - \text{tr} \textbf{T}^{\mathrm{2}}) . \end{array} $$

(36)

3 Finite element formulation

The previously presented continuum framework for the coupled mechanical-diffusion-degradation model leads to a set of coupled differential equations, which are time-dependent and highly nonlinear. For further procedure, a numerical approximation, the FEM, is used. We apply the Galerkin method, whereby the balance equations are represented in their weak form and weighted by independent test functions. For the fully discrete problem, we solve for the displacements u = x −X and the concentrations c_γ on the basis of the balance of linear momentum and the balance of mass for the concentrations.

3.1 Weak form

In the following we present the weak form of the balance of momentum weighted by the independent test functions for the displacements δu and the weak form of the balance of mass of the concentrations weighted by the independent test functions for the concentrations δ c_γ. The equations are posed in the reference configuration. On the one hand, the weak form of the balance of linear momentum results in

$$ \begin{array}{lllll} \displaystyle{\int\limits_{\text{B}_{0}}}\textbf{P} : \text{Grad} \delta\textbf{u} \text{dV} = {\int\limits_{\partial\text{B}_{0}}} \textbf{P} \textbf{N} \cdot \delta\textbf{u} \text{dA} , \end{array} $$

(37)

wherein the mapping of the Cauchy stress T to the first Piola-Kirchhoff stress tensor P is applied with the transformation P = JTF^-T and wherein N is the material normal unit vector on the Neumann boundary. On the other hand, the weak form of the balance of mass is given as

$$ \begin{array}{lllll} \displaystyle{\int\limits_{\text{B}_{0}}} (\dot{\text{c}}_{\gamma} \delta \text{c}_{\gamma} - \textbf{j}_{\gamma}\cdot \text{grad} \delta \text{c}_{\gamma}) \text{J} \text{dV} \\ = \displaystyle{\int\limits_{\partial\text{B}_{0}}} \textbf{j}_{\gamma} \cdot \delta \text{c}_{\gamma} \text{J} \textbf{F}^{\text{-T}} \textbf{N} \text{dA} , \end{array} $$

(38)

with JF^-TN dA = n da. We describe the Dirichlet boundary conditions as follows

$$ \begin{array}{llllllll} \textbf{u} &=& \textbf{u}^{*} && \forall && \textbf{X} \in &{\partial\text{B}_{0}^{\text{u}}} \\ \text{c}_{\gamma} &=& \text{c}_{\gamma}^{*} && \forall & &\textbf{X} \in &{\partial\text{B}_{0}^{\text{c}_{\gamma}}} , \end{array} $$

(39)

and the Neumann boundary conditions as

$$ \begin{array}{llllllll} \textbf{P} \textbf{N} &=& \textbf{t}^{*} && \forall && \textbf{X} \in &{\partial\text{B}_{0}^{\text{t}}}\\ \textbf{J}_{\gamma} \cdot \textbf{N} &=& \text{J}_{\gamma}^{*} && \forall & &\textbf{X} \in & {\partial\text{B}_{0}^{\text{J}_{\gamma}}} , \end{array} $$

(40)

with $\textbf {J}_{\gamma }= \text {J} \textbf {j}_{\gamma } \textbf {F}^{\text {-T}}$. The initial conditions are given with u(t₀) = u₀ and $\text {c}_{\gamma }(\text {t}_{0}) = \text {c}_{\gamma }^{0}$.

3.2 Discretisation in space

For space discretisation, we apply the isoparametric concept which is based on approximating geometry, displacement and concentrations by the same set of ansatz functions h^I(ξ). The discrete form of the test functions for the displacement δu^h and for the concentrations $\delta {\text {c}_{\gamma }^{\text {h}}}$ results in

$$ \begin{array}{llllllll} \delta\textbf{u}^{\mathrm{h}} &=& {\sum}_{\text{I}=1}^{\text{NN}} \text{h}^{\mathrm{I}} (\boldsymbol{\xi}) \delta\textbf{u}^{\mathrm{I}} \\ \delta\text{c}_{\gamma}^{\text{h}} &=& {\sum}_{\text{I}=1}^{\text{NN}} \text{h}^{\mathrm{I}} (\boldsymbol{\xi}) \delta{\text{c}_{\gamma}^{\text{I}}} , \end{array} $$

(41)

wherein NN denotes the number of nodes per element and where ξ represents the local coordinates. In this work, we apply an eight-noded element description, i.e. two-dimensional Serendipity elements under plane strain conditions. Moreover, approximations of the degrees of freedom as well as of all related gradient operations follow straightforwardly.

3.3 Discretisation in time

The simulation of the diffusion of the concentrations requires a discretisation in time for which we apply the Newmark-beta method. Within the considered time interval, we approximate a constant average acceleration of the concentrations $\ddot {\text {c}}_{\beta }$, i.e.

$$ \begin{array}{llllllll} \ddot{\text{c}}_{\beta} = \frac{1}{2} (\ddot{\text{c}}_{\gamma}(\text{t}) + \ddot{\mathrm{c}}_{\gamma}(\text{t} + {\Delta}\text{t})) , \end{array} $$

(42)

with the previous acceleration of the concentrations $\ddot {\mathrm {c}}_{\gamma }(\text {t})$ and the acceleration to be approximated in the present time step $\ddot {\mathrm {c}}_{\gamma }(\text {t} + {\Delta }\text {t})$. Based on this, one obtains

$$ \begin{array}{lll} \dot{\mathrm{c}}_{\gamma}(\text{t} + {\Delta}\text{t}) &=& \dot{\mathrm{c}}_{\gamma}(\text{t}) + {\Delta}\text{t} \ddot{\text{c}}_{\beta}\\ \text{c}_{\gamma}(\text{t} + {\Delta}\text{t}) &=& \text{c}_{\gamma}(\text{t}) + {\Delta}\text{t} \dot{\mathrm{c}}_{\gamma}(\text{t}) + \frac{1}{2} {\Delta}\text{t}^{2} \ddot{\text{c}}_{\beta} . \end{array} $$

(43)

Finally, we gain the approximation of the velocity for the concentrations $\dot {\mathrm {c}}_{\gamma }(\text {t} + {\Delta }\text {t})$ in the present time step, i.e.

$$ \begin{array}{lllllll} \dot{\mathrm{c}}_{\gamma}(\text{t} + {\Delta}\text{t}) &=& \frac{\mathrm{2}}{\Delta\text{t}} (\text{c}_{\gamma}(\text{t} + {\Delta}\text{t}) - \text{c}_{\gamma}(\text{t})) - \dot{\mathrm{c}}_{\gamma}(\text{t}) . \end{array} $$

(44)

4 Numerical analysis examples

In this section, we discuss the properties of the proposed diffusion controlled degradation model in the context of representative examples focusing on the basic coupling mechanism between diffusion impact and structural response.

4.1 Material behaviour

In the following, the correlations between mechanical and degradation processes are investigated on the basis of an analysis including one finite element and two different sets of boundary conditions as shown in Fig. 2. Furthermore, we apply the material and loading parameters shown in Table 1.

Fig. 2

a Dirichlet boundary conditions Δu. b Neumann boundary conditions F

Table 1 Material and loading parameters

The degradation process is identically established for both examples, at this point independently of chemical concentrations, with the change of the referential mass density ρ₀, as shown in Fig. 3. With this at hand, the degradation part of the deformation gradient F^d can be evaluated with (2).

Fig. 3

Approach for the time-dependent decrease of density with $\rho _{0}/\rho _{0}^{*}=1-{5\times 10^{-4}} \text {t}^{2} [\text {s}]^{-2}$, wherein t is time

For mechanical impact, the first example is uniaxially loaded with a displacement Δu = 0.5cm in y-direction, see Fig. 2a, which results in the deformation gradient

$$ \begin{array}{@{}rcl@{}} \textbf{F} = \textbf{1} + \varepsilon \textbf{e}_{\mathrm{y}} \otimes \textbf{e}_{\mathrm{y}} , \end{array} $$

(45)

wherein ε is the strain in y-direction. For comparison, the second example is loaded with a constant traction with resulting force F = 0.015MN in y-direction, see Fig. 2b.

The determinant of the degradation contribution to the deformation gradient, J^d, follows directly from the time-dependent approach for the change of referential density, see Fig. 4, so that we observe the same referential mass density evaluations for both boundary conditions. However, the total deformation changes differently in the two examples. For the boundary condition displayed in Fig. 2a, the total deformation represented by J is constant, because the structure does not allow shrinkage. The elastic deformation changes inversely so that, i.e. J^e = [J^d]^− 1J. The stress T_yy increases, because the total deformation is fixed, whereas the degradation part decreases. In comparison to the first example with fixed Dirichlet boundary conditions, the example illustrated in Fig. 2b results in a degradation process and the total deformation represented by J decreases over time. In return, the elastic deformation J^e and the stress T_yy are constant. The related states of deformation are additionally illustrated in Fig. 5.

Fig. 4

Evaluation of the boundary condition Fig. 2a in red and of the boundary condition Fig. 2b in blue with a time slot of 29 s

Fig. 5

The grey colour shows the initial state, the black colour the deformation state and the blue colour the degraded structure: displacement loading (left, a) and force loading (right, b) by analogy with Fig. 2

4.2 Diffusion driven degradation

In this example, we discuss the coupled effects of diffusion leading to a degradation of the material. This example is motivated by the idea to analyse chemical influence on hollow concrete blocks, since concrete is a porous medium that allows the inflow of concentrations and is susceptible to chemical degradation. The boundary value problem is illustrated in Fig. 7a, which shows a structure with a hole. The hole of the structure is defined by the parameters (s₁,s₂), which represent the axes of an ellipse. We consider a concentration inflow from the left side of the structure. The choice of parameters, as stated in Table 2, is based on values common for concrete and the orders of magnitude for chemical diffusion processes, although no specific chemical process is described here. The material degradation triggered by molecular processes with concentrations occurs very slowly. Therefore, the calculation is accelerated by considering the diffusion rate per day (d) and a time-dependent increase in concentrations per day, as shown in Fig. 6. In total, the computation considers a time period of 4 d with a time step size of 1 d. A stable mechanical environment is enabled by fixed displacements as shown in Fig. 7a. The concentration inflow on the left side leads to the contour plot given in Fig. 7b. Figure 8a illustrates the degradation induced by the concentrations and in Fig. 8b the impact of the deformation on the first principal stress T₁ is evaluated. Since material degradation is triggered by a high concentration of chemicals, the main material reduction takes place close to the inlet area.

Table 2 Material parameters for the structure with the hole

Fig. 6

Time dependent Dirichlet boundary conditions for the concentrations $\text {c}_{\gamma }^{*}(\text {t})$

Fig. 7

a Dirichlet boundary conditions of a structure with a hole. b Contour plot of the concentrations after the first time step

Fig. 8

a Deformed mesh after 4 d printed in red vs. the initial area outlined in black. b First principal stress T₁ after 4 d, with a maximum stress of 1.7 × 10^− 1MNcm^− 2

The maximum of the first principal stress is located in this area close to the hole. Overall, this example clearly shows the coupling between diffusion and deformation.

5 Structural optimisation framework

This section outlines the connection between structural analysis and optimisation framework. Furthermore, the information required for the optimisation problem, such as the objective function, the constraints and design parameters, are presented. The main programme runs in the numerical computing environment of MATLAB which contains a link to the open software gmsh to create a mesh on the one hand and an interface to a Fortran code using MEX-file interfaces on the other. The linearisation of the weak form, in discrete form, is implemented in a Fortran based FE code and the MEX-file transfers the data to the workspace. The assembly and calculation are implemented in MATLAB.

5.1 Sketch of algorithmic framework

An overview of the algorithm is shown in Fig. 9. The algorithm can be divided into three main sections, the boundary value problem, structural analysis and structural optimisation. First, the boundary value problem is relevant for the definition of the model problem by means of parametric geometry description. The information about the boundary value problem provides the input for the structural analysis. In addition, the solution of the structural analysis computes a structure deformed by chemical and mechanical loads, which is used as initial design in the structural optimisation. In this paper, continuum mechanical quantities such as the evaluation of the stress restriction and geometry parameters are selected for the objective function and the constraints. These variables are calculated within the framework of structural analysis and directly passed on to the mathematical solver. The structural optimisation minimises the defined objective function while maintaining the given constraints and provides new design parameters. The optimisation task is solved with the help of a MATLAB toolbox.

Fig. 9

Illustration of the algorithmic framework

Furthermore, the technical implementation can be briefly summarised. The global programme runs in MATLAB, while the information at element level is implemented in a Fortran code and embedded via MEX interfaces. Information at element level retrieves the coupled equations for the description of the problem presented, the variation formulations with the discrete weak forms and the gradient information for the FEM solution. In addition, an interface to the open software gmsh, a finite element mesh generator, is generated during the framework of the boundary value problem. In this paper, geometrical values are chosen as design parameters; therefore, a parametric mesh design is implemented. The functions introduced as part of the optimisation process are described in the typographic style ‘italics’, e.g. objective function J, constraints g, design parameters s.

5.2 Optimisation problem

The structural optimisation problem is solved by using the nonlinear optimisation function fmincon provided by MATLAB toolbox, see (MathWorks 2019). This solver finds the minimum of a restricted nonlinear and multivariable function. The general optimisation problem follows with an objective function J(v,s), nonlinear inequality constraints g(v,s), upper and lower limit values s^u and s^l for a set of design parameters, i.e.

$$ \begin{array}{ccccll} \min J(\textbf{v}, \boldsymbol{s})&: & \boldsymbol{g}(\textbf{v}, \boldsymbol{s}) &\le& 0 & \text{constraints}\\ & & \boldsymbol{s}^{\mathrm{l}} &\le& \boldsymbol{s} \le \boldsymbol{s}^{\mathrm{u}}& \text{limit} \text{values} , \end{array} $$

(46)

wherein v ∈{u,c_γ} are the field variables and s geometric design parameters. In this paper, we deal with the two shape optimisation problems as follows:

• 1. Minimisation of the maximum first principal stress $\text {T}_{1}^{\max \limits }$ by changing geometrical parameters, with the constraint of a maximum loss of area.

• 2. Minimisation of the area by changing geometrical parameters, with the constraint of an upper limit for the maximum first principal stress $\text {T}_{1}^{\max \limits }$.

In both optimisation setups, structural analysis is applied to calculate the objective and constraint functions, since there is a dependency on the field variables. Further specifications on the objective function and constraint are provided in Sections 6.1 and 6.2. Fmincon uses the ‘sqp-legacy’ algorithm to solve the optimisation task, wherein the gradients of the objective function and constraints are obtained numerically by finite differences. The Hessian matrix is iteratively integrated using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) procedure. Further details can be found in (MathWorks 2019). The application of the finite difference method is a precise but time-consuming method. In the outlook, we address a more efficient approach for future work.

6 Numerical optimisation examples

In the following, two examples are presented for shape optimisation of structures that are loaded by chemical concentrations. The first example is taken from Section 4.2 and shows a hollow concrete block which is loaded by chemical substances. We calculate the optimised shape of the brick in such a way that the maximum stress caused by the concentrations is reduced while still retaining material. The second example is inspired by a mechanically loaded bridge, additionally loaded with chemical substances, which can be caused by, for example calcium leaching. The goal of this example is a structure with minimised material and limited stresses.

6.1 Optimisation of a structure with a hole

Using the example of the structural analysis in Section 4.2 with the material parameters from Table 2, an optimised shape is discussed below. We use the sum of the Gaussian point values of the first principal stress in the maximum loaded element $\text {T}_{1}^{\max \limits }$ as the objective function J(v,s) and the change of the area as inequality constraints g(v,s) with

$$ \begin{array}{lll} J(\textbf{v},\boldsymbol{s}) &=& {\sum}_{1}^{\text{NG}} \text{T}_{1}^{\max} \\ g(\textbf{v},\boldsymbol{s}) &=& \left|\frac{\text{A}^{\text{ini}} - \text{A}}{\text{A}^{\text{ini}}}\right| - \bar{\mathrm{A}} , \end{array} $$

(47)

wherein NG is the number of Gaussian points, Aⁱⁿⁱ is the initial area, A is the actual area and $\bar {\mathrm {A}}$ is a threshold. The optimisation problem follows as

$$ \begin{array}{ccclll} \mathrm{ min} J (\textbf{v}, \boldsymbol{s})&: & g(\textbf{v},\boldsymbol{s}) &\le& 0 & \text{constraints} \\ & & \boldsymbol{s}^{\mathrm{l}} &\le& \boldsymbol{s} \le \boldsymbol{s}^{\mathrm{u}} & \text{limit} \text{values}. \end{array} $$

(48)

The axes of the ellipse (s₁,s₂) are the design parameters. The optimisation algorithm applies the input parameters listed in Table 3.

Table 3 Input parameters for the optimisation algorithm

The following diagramme, see Fig. 10, shows the iteration course of the optimisation to minimise the objective function while fulfilling the constraints. The solver requires a total of 9 iterations until an optimal shape is found to minimise the average first principal stress in the Gaussian points to 0.083MNcm^− 2.

Fig. 10

Iteration of the optimisation solver ‘sqp-legacy’, which shows the decrease of the objective function, i.e. the first principal stress

After 18,468 s of computing time, the optimisation leads to the following new design parameters, i.e. (Fig. 12).

$$ \boldsymbol{s} = \begin{bmatrix} 5; 5.525 \end{bmatrix}, $$

(49)

where the length of the axes of the hole change to reduce the first principal stress.

The contour plot of the first principal stress in Fig. 11 shows that the change of the hole parameters, in other words the design parameters, leads to an overall decrease of the first principal stress. Thereby, the constraint is fulfilled, i.e.

$$ g(\textbf{v},\boldsymbol{s}) = \left|\frac{\text{A}^{\text{ini}} - \text{A}}{\text{A}^{\text{ini}}}\right| - \bar{\mathrm{A}} = -2.7\times10^{-6} \le 0, $$

(50)

Fig. 11

Evaluation of the first principal stress with new design after 4 days

This example illustrates how shape has a major influence on the effects of degrading concentrations and that damage can be minimised by changing the shape. This observation is supported by Fig. 13 which shows the temporal course of the deformation based on the determinant of the degradation gradient, J^d, at the maximum loaded point. The curve is compared on the one hand with the initial shape, and on the other hand with the optimised shape, showing that the optimised form leads to a lower degradation.

Fig. 12

Evaluation of the new design. a Initial structure vs. b optimal structure after 9 optimisation steps

Fig. 13

Evaluation of the determinant of the degradation gradient J^d in the marked point within 4 days, the blue line refers to the initial and the red line to the optimised design

6.1.1 Remarks on numerical investigation

In this paper, the mathematical optimisation is solved utilising the MATLAB function fmincon, which calculates the minimum of a constrained nonlinear multivariable function, see MathWorks (2019). The user can chooses between different algorithms to solve the task. In this section, the influence of the solution algorithm on the optimised result is evaluated. For this purpose, the result of the example ‘Optimisation of a structure with a hole’ is evaluated using three different algorithms. The ‘sqp-legacy’ and ‘active-set’ algorithm are based on sequential quadratic programming (SQP) method. Thereby, SQP is derived using the Newton method and taking into account inequality constraints. The two algorithms differ in their implementation, e.g. they apply different definitions for the strict feasibility with respect to bounds or the choice of the solution algorithm for the subproblems. In contrast, the third algorithm ‘interior-point’ combines two different approaches to solve the optimisation task. The algorithm uses Newton steps or conjugated gradient steps depending on the solution of each iteration step.

Table 4 compares the efficiency of the algorithms and shows that algorithm ‘sqp-legacy’ runs most efficiently, because the solution is obtained after only 7 iterations, respectively after 18,468 seconds. The solutions are identical for all approaches.

Table 4 Comparison of different solution algorithms

Furthermore, Fig. 14 shows the optimised design taking into account different mesh sizes. The identical optimised parameters with s= [5;5.525] are determined, but the computing time increases with the refinement of the mesh. In detail, the finer mesh shown in Fig. 14 needs 3.3 times more calculation time than the coarse mesh.

Fig. 14

Optimisation results under consideration of different mesh

6.2 Optimisation of a bridge-like structure

This example is inspired by a bridge with mechanical loads under environmental influences, such as calcium leaching. This scenario can occur, for example during long-term exposure to pure water, which triggers the diffusion of calcium ions. To illustrate whether environmental influences are taken into account or not, the calculation is listed in the first step without and in the second step with the influence of chemical concentrations. In order to reduce the material costs, the objective J(v,s) of the problem is to minimise the area within a plane strain setting while holding a threshold for the first principal stress, i.e. g(v,s). The optimisation problem is obtained as follows

$$ \begin{array}{lll} J(\textbf{v},\boldsymbol{s}) &=& {\sum}_{\text{e}=1}^{\text{NE}}\displaystyle{\int\limits_{\text{B}_{0}^{\text{e}}}} \text{dV} \\ g(\textbf{v},\boldsymbol{s}) &=& \mathrm{T}_{1} - \text{T}_{1}^{\max} , \end{array} $$

(51)

wherein the area is calculated by the sum of the total element volumes. We apply a threshold for the maximum first principal stress $\text {T}_{1}^{\max \limits }$, thus the total optimisation problem follows with, i.e.

$$ \begin{array}{ccrcll} \min J(\textbf{v},\boldsymbol{s})&:& \boldsymbol{g}(\textbf{v},\boldsymbol{s}) &\le& 0 & \text{constraints}\\ & & \boldsymbol{s}^{\mathrm{l}} \le \boldsymbol{s} &\le& \boldsymbol{s}^{\mathrm{u}} & \text{limit} \text{values}. \end{array} $$

(52)

To save computing time, the symmetry of the structure is utilised and the calculation is performed on half of the system using symmetry boundary conditions on the right side, see Fig. 15. Table 5 presents the applied material parameters. The applied forces with F = 10kN lead to the first principal stress distribution, which is displayed in the contour plot in Fig. 18b. The formation of a tensile area in red and compression area in blue becomes visible.

Fig. 15

Mechanical boundary conditions of the structure

Table 5 Material parameters for the bridge-like structure

The lower edge of the structure is defined by a B-spline function with four control points, where the design parameters are the control points that allow vertical displacement, see Fig. 16, i.e.

$$ \boldsymbol{s} = \begin{bmatrix} \text{x}_{\mathrm{1}};\text{x}_{\mathrm{2}};\text{x}_{\mathrm{3}};\text{x}_{\mathrm{4}} \end{bmatrix}. $$

(53)

Utilisation of symmetry must also be taken into account for the calculation of design parameters during optimisation. For this purpose, the condition x₃ = x₄ is introduced.

Fig. 16

Design parameters for the structural optimisation

The parameters for the optimisation algorithm are shown in Table 6, where s^u and s^l are the matrices containing the minimum and maximum allowable change of design parameters. Figure 17 shows the iteration of the optimisation solver with the decrease of the objective function, i.e. the area of the structure. After 6 iterations the objective function converges, the side condition is fulfilled and a local minimum is given with the following design parameters

$$ \boldsymbol{s} = \begin{bmatrix} 1;0.893;-0.1372;-0.1372 \end{bmatrix}. $$

(54)

Table 6 Input parameters for the optimisation algorithm

Fig. 17

Iteration of the optimisation solver ‘sqp-legacy’, which records the decrease of the objective function, i.e. the area

The design parameters reveal that the largest saving occurs in the area of the least stress, i.e. in the area where neither compressive nor tensile stress is present, whereby more material is required in the middle of the structure to ensure load-bearing capacity. In total, the optimised design saves 20.6 % of the area. The contour plot in Fig. 18 represents the first principal stress in the optimised and initial design.

Fig. 18

Evaluation of the first principal stress induced by mechanical load. a Optimal structure with $\text {T}_{1}^{\max \limits }= 120 \text {kN} \text {m}^{-2}$ vs. b initial structure with $\text {T}_{1}^{\max \limits }= 112 \text {kN} \text {m}^{-2}$

In a second step, the optimised bridge is additionally loaded by chemical concentrations. In this example, we focus on the general effect of any chemical concentrations that trigger material degradation. Additionally, to the mechanical force, two concentration inflows are located on the top of the structure and are provided by a concentration increase of 0.33molm^− 3 per time step, see Fig. 19.

Fig. 19

Mechanical and concentrations boundary conditions of the structure

As already mentioned in the example in Section 4.2, the chemically induced degradation process is accelerated by running the simulation in time steps of days. The contour plots in Fig. 20 show the distribution of the concentrations and the resulting influence on the determinant of growth after 9 days. Figure 21 shows the time progression of the first principal stress in one nodal point, P₁, where the maximum first principal stress occurs. The first principal stress increases due to the increased concentrations. The previously defined maximum stress of $\text {T}_{1}^{\max \limits }= 120 \text {kN} \text {m}^{-2}$ can no longer be maintained.

Fig. 20

a Evaluation of the concentrations and b the impact on the determinant of the degradation part of deformation after 9 days

Fig. 21

Increase of maximum first principal stress over 9 days with chemical impact

The aim of the optimisation is to save as much material as possible while still ensuring the load capacity. For this reason, smaller deviations in stress due to diffusion processes are also relevant, since otherwise the load-bearing capacity can no longer be guaranteed over a long load duration.

By analogy with the previous example, the points of the polygon chain at the lower edge of the structure are design parameters. Using the same parameters for the optimisation algorithm as provided in Table 6 and starting from the initial design as illustrated in Fig. 15, the optimisation yields the following optimal design parameters

$$ \boldsymbol{s} = \begin{bmatrix} 0.933;0.793;-0.196;-0.196 \end{bmatrix}. $$

(55)

The iteration process illustrated in Fig. 22 shows how the design parameters reduce the objective function by 17.4 % while maintaining the constraints. The increased maximum first principal stress due to the degradation process changes the optimal design and allows for less material saving, as shown in Fig. 23. It is clear that the consideration of environmental conditions, such as calcium leaching, are necessary to predict long-term performance.

Fig. 22

Iteration of the optimisation solver ‘sqp-legacy’ for the example coupled to chemical impact, which shows the decrease of the objective function, i.e. the area

Fig. 23

Comparison of the results: the green coloured area shows the optimal result when chemical influence is considered and the red area shows the result without degradation processes

7 Summary

In this paper, a coupled mechanical-diffusion-degradation model is presented. The degradation process is derived by a multiplicative split of the deformation gradient and a constitutive approach for the development of growth, which assumes chemical concentrations as a trigger for degradation. The numerical FE framework is briefly presented. Furthermore, the embedding of a structural optimisation framework is outlined. The applicability of the model is presented for a structure with a hole and a beam for which a practical reference is outlined. The main focus of the examples lies in the analysis of the influence of long-term acting chemical concentrations, which can influence the mechanical stress and can keep certain upper limits of the stress by constructive changes.

For future work, an alternative to the numerical finite difference method will be applied, namely the variational sensitivity analysis as outlined in Barthold and Stein (1996), which accelerates the simulation time.

Overall, it is shown that the model can provide an optimal design taking into account long-term effects from concentrations that damage the material, while still maintaining certain limits for the load-bearing capacity. The algorithm offers the possibility to integrate different chemical processes, to calculate the interaction with the mechanical behaviour and to solve optimisation problems.