Startseite The effect of plasticity on damage evolution using a relaxation-based material model
Artikel Open Access

The effect of plasticity on damage evolution using a relaxation-based material model

  • Stephan Schwarz EMAIL logo , Klaus Hackl und Philipp Junker
Veröffentlicht/Copyright: 14. Dezember 2018
Veröffentlichen auch Sie bei De Gruyter Brill

Abstract

As damage occurs in the context of high stresses that are also related to the presence of plastic strains, it is natural to investigate the effect of plasticity on damage evolution and to thus achieve a more realistic model. In this work, the existing and new damage model presented in [Junker P, Schwarz S, Makowski J, Hackl K. Continuum Mech. Therm. 2017, 29 (1), 291–310] is enhanced with plasticity and isotropic hardening. The damage model is based on a relaxation-based approach and does not require additional complex regularization techniques besides considering viscous effects. The benefit of the model are mesh-independent results for the rate-dependent case, even without considering, e.g. gradient terms for mathematical regularization. The enhancement with plasticity and isotropic hardening was investigated for a representative volume element that considerd a damaging matrix material and non-damaging hard precipitates. Two different loading types, pure tension and pure shear, yielded the homogenized stress/strain response for the material at various loading rates. Hereto, several finite discretizations in terms of finite-element meshes were used. The results underline the mesh-independence for physically reasonable loading rates and viscosities.

1 Introduction

Crack and damage processes occur at different length scales and have to be considered for (nearly) all failure investigations of designed materials, i.e. those which are composed of crystallographic phases that possess varying material properties, as well as entire construction parts. Distinct cracks result from accumulating dislocations that ultimately trigger the material to lose its binding on the inter-atomic scale. This results in the separation of the material domains and in an increase in the domain surface. Damage, on the other hand, describes this process in a smeared manner. Due to a local network of cracks, the stiffness of the material is reduced until no load can be carried at the limit of propagating crack length and number. Damage modeling focuses exactly on the local evolution of decreasing stiffness without considering the underlying reasons in detail. Regardless of whether crack or damage processes are investigated, it is nearly impossible to observe cracks without the evolution of dislocations in metallic materials. Furthermore, the evolution of the movement of dislocations is associated with irreversible deformations that are usually referred to as plastic deformations. Consequently, for the modeling and simulation of damage processes, it is essential to take plastic deformations into account as well.

In this paper, we present the implementation of plasticity with isotropic hardening in the damage model published in [1], [2]. This model provides a relaxation-based approach to damage modeling without the requirement of any distinct regularization procedures. Such as the use of common field functions, inclusion of gradients or complex integration techniques. In order to achieve this advantage, slight modifications were applied to the relaxed (condensed) energy of the undamaged and damaged material and viscous effects were added to the model. As no additional quantity is needed to consider the spatial behavior, which is something common for models based on, e.g. gradient-enhanced regularization strategies (cf. [3], [4]), the number of unknowns is reduced, which leads to several numerical advantages including a considerably reduced computational effort and an improved convergence behavior. However, as this paper focuses on the coupling of damage and plasticity, the interaction phenomena of both material characteristics are of particular importance and, furthermore, the performance and efficiency of this new coupled damage-plasticity model will be presented by fairly realistic numerical calculations.

First, it is important to provide some information about softening effects and the reasons for the need of regularization techniques. Common material models that take account of softening effects due to damage encounter the problem of ill-posed boundary value problems (in the case of no regularization). This condition leads to a nonunique solution for the resulting algebraic system, and the numerical results exhibit a strong mesh dependence in a finite-element context. For this reason, regularization techniques can be applied to prevent this problem by taking account of the nonlocal behavior of the damage. This can be achieved, e.g. by using a field function to couple the local damage parameter to a nonlocal level in which differences between the local and nonlocal parameter as well as the gradient of the nonlocal parameter are penalized [4]. Due to the additional field function there are more unknowns at the nodes, and consequently, the calculations are more expensive. In contrast, the new damage model [1], [2] is a two-phase model in which the transition between the undamaged and damaged energy states is carefully chosen: the energy is defined between the Reuß bound, a linear transition between the two energy states which provides a smooth transition between undamaged and damaged material, and the Voigt bound, a sudden transition in the crossing point of the two energy states. The Reuß bound provides a convex energy so that the mathematical problem becomes well posed and the results are mesh-independent; however, it does not allow a stress decrease, which is what we would expect in the context of damage problems. This is why we modified the energy slightly towards the Voigt bound, which allows a stress decrease and reductions of the mesh-independence are prevented by adding viscous effects to the model.

In order to extend this damage model [1], [2] to further applications, this paper is devoted to the implementation of plasticity phenomena. It is natural to assume plasticity in a damage model because damage occurs in the context of high stresses that are also related to the development of plastic strains. In reality, plastic strains continue to evolve until the material fails, which in turn can be described by a completely damaged configuration. As damage and plasticity phenomena often occur together, they should be considered by means of a coupled model. Therefore, isotropic plasticity is taken into account and coupled with the existing model. This is achieved by adding a further internal variable that captures the plastic strains εp and by extending the energy by an additional term 1/2H^α2 that describes the isotropic hardening. Essentially, all other equations result from both extensions and can be determined by applying the principle of the minimum of the dissipation potential [5], [6], and [7]. The merging of a pure damage model and a pure plastic model to produce a coupled one is schematically plotted in terms of a force vs. displacement diagram in Figure 1.

Figure 1: Schematic plot of a force vs. displacement diagram for damage (left), plasticity and isotropic hardening (center) and both (right).
Figure 1:

Schematic plot of a force vs. displacement diagram for damage (left), plasticity and isotropic hardening (center) and both (right).

Although this is the first paper that extends the new damage model [1], [2], there are several similar approaches based on other damage models. Basically, regularized damage models can be separated into three groups: the gradient type, the integral type and the viscous regularization to which our model belongs.

The integral type introduces nonlocal variables in order to average the local variables of the defined surroundings of a specific point. Due to the definition of the surroundings, the material length is a prescribed parameter. In the literature, there is a comprehensive discussion of several integral-type damage approaches with and without consideration of the plasticity [8]. Furthermore, in [9], a coupled damage-plasticity model has been turned into a nonlocal model by applying integral formulations to energy terms that are associated with the damage processes. Moreover, in [10], a nonlocal damage model has been coupled with the plasticity, in which the damage is linked to the evolution of plastic strains and was regularized by a weighted spatial averaging of the damage-driving variable. Another example of coupling plasticity and nonlocal integral-type damage is [11], in which plastic-damage dissipation consists of a plastic-independent part and a coupled part. A recent paper [12], presents a nonlocal coupled damage-plasticity model based on integral regularization techniques with a focus on the determination of model parameters with a novel calibration procedure related to experimental data.

The gradient-type models apply a gradient formulation on a nonlocal quantity in order to influence or determine its evolution, which in turn is connected to the local quantity for regularization purposes. A typical example for gradient-type regularization models with plasticity can be found in the literature in [4]. This model is enhanced by a nonlocal field function in the energy terms to regularize the damage parameter, whereby the same procedure is applied to the plasticity. In [13] the plasticity is considered within the framework of nonlocal damage by means of the dissipation function and two different yield criteria. In contrast, the work of [14] regularizes by introducing the second-order gradient of the plastic multiplier in the yield function to achieve the gradient-dependent softening plasticity. There is also the idea in [15] to build up an over-nonlocal enhancement of a gradient-type plastic-damage model by considering the weighted sum of local and nonlocal parameters, which corresponds to the idea of integral-type regularizations. To introduce a crack propagation focused work with implicit formulation, one can mention [16], in which the long-range interactions of an integral-type model were combined with the computational efficiency of a gradient formulation by using a nonlocal effective strain measure. In the context of crack prorogation, it is also interesting to mention the work in [17], in which an adaptive finite-element algorithm was implemented on a gradient-based, regularized damage-plasticity model. Borst et al. [18] present several combinations of coupling scalar and gradient damage with gradient and hardening plasticity. A focus on algorithmic aspects is given in [19], which presents gradient-enhanced damage and plasticity approaches with the ability to model localization phenomena in quasi-brittle and frictional materials. In order to consider the strain delay effect, a rate-dependent model for quasi-brittle materials that includes damage and multi-variable plasticity under dynamic loading was recently published in [20]. In [21], a coupled damage-plasticity model, which is implemented within the stress-based variational formulation providing a very robust behavior, is presented. Furthermore, a robust and consistent model can be found in [22] for finite-strain elastoplasticity with and without nonlocal damage by using remeshing with a focus on the transfer of history data from one mesh to another. The concept of effective stress and strain equivalence and the stress triaxiality for considering ductile failure is shown in the nonlocal damage-plasticity framework in [23]. In the literature one can also find [24] in order to simulate the failure behavior of concrete, in which a thermodynamically consistent nonlocal gradient and fracture energy-based plasticity model is used that also realistically predicts the transition from brittle to ductile post-peak response by the use of two characteristic lengths. Last but not least, a plastic nonlocal damage model for cementitious materials based on a regularization technique with use of a damage Laplacian and limit function and on a Drucker-Prager plastic limit function that considers isotropic hardening is presented in [25].

Viscous models use the fact, that rate-dependent effects delay the evolution of internal variables, which in the present case of damage modeling can be interpreted as a damage limiter for a certain finite element. This forces the surrounding elements to evolve damage as well to account for the same amount for dissipated energy as in an unregularized, mesh-dependent simulation. Consequently, the delaying effect of viscosity acts as a regularization strategy. In contrast to integral- or gradient-type models, there is no spatial internal length scale but there is an internal length scale in time which corresponds to a spatial description. Viscous regularization in the context of damage models originates from the work of Needleman [26] who first showed that material rate dependence indirectly introduces a length scale and pathological mesh sensitivity can be eliminated for quasi-static and dynamic problems. He also encountered numerical instability problems for approaching the rate-independent limit. Viscous regularization for localization problems was also used in [27] and applied to a cracked medium. It was also shown that the introduced implicit internal length scale smears the crack and the width of the localization zone becomes independent from the mesh for high strain rates. In [28], a rate-independent damage model was transformed into a rate-dependent one in order to catch more effects of the loading rate for concrete. It was furthermore possible to show the well posedness and mesh objectivity of numerical investigations. The transformation into a rate-dependent model was inspired from the Perzyna viscoplasticity model and they were able to show the same responses with the rate-dependent model as for the rate-independent one for very small loading rates. Viscoplasticity was also used for regularization in [29] where further and deeper studies for this kind of regularization were performed with a focus on ductile materials and the use of two different hardening models. Also plasticity models can be regularized by a viscous term shown in [30] where the relation between the material length scale and the width of shear bands was examined. In order to limit the rate dependency for viscous regularization, a new approach was introduced in [31], which was applied to interface debonding models whereas a damage model including a delay effect was investigated in [32].

For classification purposes, our new coupled model basically belongs to the third group: damage models with viscous regularization. Advantages are that such models do not rely on complex regularization procedures and are therefore rather easy to implement. Consequently, the presented model directly provides mesh-independent results and, due to its good numerical performance, it also does not need any adaptive mesh techniques such as those mentioned in the literature. One result of being free of gradients is that the model does not possess an internal length scale, due to the added viscous effects it has rather an internal length scale in time. A drawback of the model is that such models include a conceptual rate-dependence. More details on the original model for pure damage without considering plastic effects can be found in [1], [2] whereas a detailed investigation regarding numerical aspects and further model development can be found in [33].

The following section describes the material model and presents the coupling of plasticity and damage. The finite-element approximation is given in Section 3 and the numerical results with several examples are subsequently given in Section 4. Finally, the coupling of damage and plasticity as well as the numerical results are discussed in Section 5, which also includes an outlook for possible future work.

2 Coupled material model

2.1 Energies

The damage approach in this model, which is based on [1], [2], is basically a two-phase model with an undamaged part and a damaged part. For this reason, we define two Helmholtz free energies, Ψ0 for the undamaged material and Ψ1 for the damaged material:

(1)Ψ0(ε):=12(εεp):E0:(εεp)+12H^α2
(2)Ψ1(ε):=12(εεp):E1:(εεp)+12H^α2

where ε is the total strain, εp is the plastic strain, α is the hardening variable, and H^ is the hardening modulus. The energies differ only in the elastic constants E0 and E1, which relate the material to the undamaged case “0” and the damaged case “1”. The stiffness of the damaged material is naturally much softer and can be described with a numerical scalar κ when assuming isotropic damage by

(3)E1=κE0with:κ1.

As already implied in the introduction, the key idea is now to establish a mixture rule for both energies Ψ0 and Ψ1 that is between the Reuß energy, which is convex but not appropriate for damage problems due to the missing drop of stresses, and the Voigt energy, which allows a drop of stresses but is no longer convex. An “appropriate” choice between these two limit cases and adding viscous effects to the model allow a stress decrease in conjunction with mesh-independent results. We chose a mixture in the following way, see [1] for more details:

(4)Ψβ=infεi{(1α^ddβ)Ψ0(ε0)+(α^d+dβ)Ψ1(ε1)|(1α^ddβ)ε0+(α^d+dβ)ε1=ε}.

Obviously, the damage variable d[0,1] is actually used as a volume fraction and describes the amount of undamaged material (d=0) and damaged material (d=1). The very small parameter α^>0 serves as a numerical quantity to ensure nonzero driving forces for ε0 and d=0, whereas β includes the nonlinear mixture and allows us to switch between the limit cases of the energies: setting α^=0 and β=1 leads to the convex Reuß energy, and setting α^=0 and β leads to the non-convex Voigt energy.

Combining both energies by applying the described mixture rule results in an energy Ψβ containing an effective elastic tensor Eeff:

(5)Ψβ=12(εεp):Eeff:(εεp)+12H^α2.

The effective elastic tensor is accordingly defined as

(6)Eeff=[(1α^ddβ)E01+(α^d+dβ)E11]1=:f(d)E0.

This in turn allows the identification of a damage function f(d)

(7)f(d)=[1+(α^d+dβ)(1κ1)]1

as well as to calculate its derivative

(8)f(d)=(α^+βdβ1)(1κ1)(1+(α^d+dβ)(1κ1))2

which will be used for the calculation of the thermodynamic driving forces. The combined Helmholtz free energy can also be written as

(9)Ψβ=12(εεp):f(d)E0:(εεp)+12H^α2.

The Gibbs energy G is given as

(10)G:=ΩΨdVΩbudVΩtudAminu

where Ω is the body’s volume and Ω its surface, b denotes the body forces (due to, e.g. gravity), and t the external tractions. The energy is minimized with respect to the displacements u, thus the balance of momentum is achieved.

The Gibbs energy can finally be set up by inserting the previously established Helmholtz free energy, reading

(11)G=Ω12(εεp):f(d)E0:(εεp)dV+Ω12H^α2dVΩbudVΩtudAminu

which is minimized with respect to the unknown displacement field u=u(x), where x is the spatial coordinate. The necessary condition for G being minimal is the stationarity condition δG=0δu, where δG denotes the Gatêaux derivative given by

(12)δG=Ωδε:f(d)E0:(εεp)dVΩbδudVΩtδudA=0δu.

Identifying the stresses to be σ=f(d)E0:(εεp), the stationarity condition results, as expected, in the weak form of the balance of linear momentum

(13)Ωδε:σdVΩbδudVΩtδudA=0δu

which can be solved by a finite-element approach in a standard manner. However, the evolution of plastic strains, the damage variable and the hardening variable remain undetermined until here. Therefore, we apply a variational concept for the derivation of the associated evolution equations in the next subsection.

2.2 Evolution of internal variables

In order to determine evolution equations for the internal variables, the principle of the minimum of the dissipation potential (see e.g. [5], [6] and [7]) is applied. Therefore, it is necessary to set up the dissipation potential L first. It consists of three parts and has to be minimized with respect to the rates of all internal variables d˙,ε˙p,α˙

(14)L:=Ψ˙β+D+consmind˙,ε˙p,α˙.

The first part is the time derivative of the Helmholtz free energy Ψ˙β, which is calculated using the chain rule to be

(15)Ψ˙β=Ψε:ε˙+Ψdd˙+Ψεpε˙p+Ψαα˙

The second part is the dissipation function D. At this point, a rate-dependent, elastoviscoplastic ansatz is used for the damage and a rate-independent ansatz for the plasticity. This demands D to be homogeneous of order one and two in d˙ and to be homogeneous of order one in ε˙p, viz.

(16)D=r1|d˙|+r22d˙2+r3|ε˙p|.

The third part considers the constraints cons, which must be fulfilled. For the present model, three constraints have to be accounted for: firstly, the damage parameter d has to be bounded, i.e. d[0,1], and may only increase. Secondly, the plastic deformations evolve in a volume-preserving way. Thirdly, the evolution of the hardening variable is coupled to the evolution of plastic strains via α˙=|ε˙p|.

The first constraint is considered by a Kuhn-Tucker parameter γ , defined by

(17)γ={γ¯:d˙<0{d˙>0d=1}0:else.

The second constraint is taken into account by a first Lagrange parameter λ1, which ensures ε˙p:I=0. Finally, the third constraint is included by a second Lagrange parameter λ2. Combination of the rate of Ψ, the dissipation function D , and the constraints gives the final dissipation potential

(18)L=Ψε:ε˙+Ψdd˙+Ψεpε˙p+Ψαα˙+r1|d˙|+r22d˙2+r3|ε˙p|+γd˙+λ1ε˙p:I+λ2(α˙+|ε˙p|)mind˙,ε˙p,α˙.

The minimization conditions for Equation (18) read

(19)Ld˙=0=Ψβd+r1d˙|d˙|+r2d˙+γ
(20)Lε˙p=0=Ψβεp+r3ε˙p|ε˙p|+λ1I+λ2ε˙p|ε˙p|
(21)Lα˙=0=Ψβαλ2.

Introducing pd=Ψβ/d=12(εεp):f(d)E0:(εεp) allows us to reformulate Equation (19) to

(22)d˙=1r2[|pd|r1]+,

whereas [x]+:=(x+|x|)/2. Details can be found in [1]. The evolution equation for the damage parameter depends on the elastic strains and thus directly depends on the present inelastic deformation state described in terms of the plastic strains.

The evolution of the plastic strains is given by Equation (20), which first has to be solved for the unknown first Lagrange parameter λ1. For this purpose, we double-contract Equation (20) with the identity tensor I, which yields under consideration of I:ε˙p=0 and I:I=3

(23)I:σ+3λ1=0λ1=13I:σ

as Ψβ/εp=f(d)E0(εεp)=σ. Consequently, Equation (20) transforms to

(24)ε˙p=|ε˙p|r3+λ2[σ13I:σ]=|ε˙p|r3+λ2devσ

with the stress deviator devσ. Equation (24) remains very cumbersome due to the presence of the absolute value of the plastic strains’ rate. It is thus convenient to perform a Legendre transformation of the dissipation function Dϵ=r3|ε˙p|, resulting in

(25)Dε=supε˙p{ε˙p:devσDελ2(α˙+|ε˙p|)}=supε˙p{|ε˙p|r3+λ2devσ:devσr3|ε˙p|λ2|ε˙p|)}+λ2α˙=supε˙p{|ε˙p|r3+λ2(|devσ|2(r3+λ2)2=:Φ~)}+λ2α˙.

Note that the constraint on the relation between the rate of the hardening variable and the plastic strain now enters with a negative sign due to the supremum we are seeking (instead of the minimum of the dissipation potential). The result of the Legendre transformation is the function Φ~0, which indicates whether plastic strains evolve or not. This enables us to reformulate Φ~ as

(26)Φ:=|devσ|(r3+λ2)!0

which is interpreted as yield function. The function Φ obviously possesses the same roots as Φ~, but has a form that is more usual in the context of plasticity. Finally, Equation (21) allows us to calculate the second Lagrange parameter, which simply reads

(27)λ2=Ψβα=H^α.

Defining the consistency parameter ρ:=|ε˙p|/(r3+λ2) closes the system of material point equations

(28)d˙=1r2[12(εεp):f(d)E0:(εεp)r1]+
(29)ε˙p=ρdevσ
(30)α˙=|ε˙p|=ρ(r3+H^α)

with the yield function Φ

(31)Φ:=|devσ|(r3+H^α)!0

and the Kuhn-Tucker conditions

(32)ρ0,Φ0,ρΦ=0.

The parameter r1 is identified to give an energetic threshold for damage evolution, which takes place in a rate-dependent way with viscosity r2. The parameter r3 is the deviator-type norm of the yield stress, whereas H^ is the hardening modulus.

3 Numerical treatment

3.1 Finite-element method

We employ a finite-element approximation to simulate finite bodies. The displacements are approximated by the usual shape functions N, whereas the strains are calculated using the strain-displacement operator B. Determination of the variations is based on the same approximation. Thus:

(33)u=Nu^,ε=Bu^,u=Nu^,
(34)δu=Nδu^,δε=Bδu^,δu=Nδu^.

Applying the finite-element approximation to the variation of the Gibbs energy δG, Equation (12), leads to the residual R

(35)R=ΩBTσdVΩbNdVΩtNdA=0

with the stress written in the usual Voigt notation as

(36)σ=f(d)E0(εεp).

The roots of R are found by employing a Newton-Raphson iteration scheme

(37)Ri+1=Ri+dRdu^Δu^=!0Δu^=[dRdu^]1Ri.

This requires the tangent operator dR/du^. It can be derived to be

(38)dRdu^=ΩBTdσdu^dV=ΩBTdσdεBdV.

The notation is converted with respect to the loading steps:

(39)dσdεdσn+1dεn+1.

For the calculation of the tangent operator, let us linearize the stresses σ=σ(ε,d,εp) by

(40)σn+1=σn+σε|nΔε+σd|nΔd+σεp|nΔεp=σn+σε|n(εn+1εn)+σd|n(dn+1dn)+σεp|n(εpn+1εpn).

The differences with respect to the current and the previous load step are denoted by the increments, and the respective derivatives are calculated as

(41)σε|n=f(dn)E0
(42)σd|n=f(dn)E0(εnεpn)
(43)σεp|n=f(dn)E0.

The tangent operator follows straight forwardly from Equation (40)

(44)dσn+1dεn+1=σε|n+σd|ndn+1εn+1+σεp|nεpn+1εn+1.

It can be specified only if the evolution equations for the internal variables according to Equations (28)–(30) are discretized in time.

3.2 Tangent related to damage

To calculate the tangent operator, the damage part requires determination of the derivative of the evolution of damage with respect to the strains dn+1/εn+1. Thus, we have to choose between, e.g. an explicit and implicit Euler integration scheme. Whereas the later one is unconditionally stable, the first one demands a simpler numerical implementation and faster computation for each time step. The onset and evolution of damage is accompanied by distinct and sharp changes of the stress/strain curve. An accurate numerical resolution of this non-monotonic curve demands the choice of rather small time increments anyways for which the explicit and implicit Euler integration schemes converge. We consequently choose an explicit Euler integration scheme and apply it to Equation (28) for simplicity, which yields

(45)dn+1=dn+{Δtr2[pd(εn+1,dn,εpn)r1]forpd>r10else

with the discretized driving force

(46)pd(εn+1,dn,εpn)=12(εn+1εpn):f(dn)E0:(εn+1εpn).

Taking the derivative simply leads to

(47)dn+1εn+1=Δtr2[pd(εn+1,dn,εpn)r1]εn+1=Δtr2f(dn)E0:εn+1.

Thus, the damage part in the tangent operator is calculated to be

(48)σd|ndn+1εn+1=f(dn)E0:(εnεpn)Δtr2f(dn)E0:εn+1=Δtr2f2(dn)(E0:(εnεpn))(E0:εn+1)

for Δd0. For Δd=0, it holds dn+1/εn+1=0.

3.3 Tangent related to plasticity

To calculate the tangent operator, the plasticity part requires determination of the derivative of the evolution of plastic strains with respect to the strains εpn+1/εn+1. Thus, we again first apply an explicit Euler integration scheme to the Equations (29) and (30), which yields

(49)εpn+1=εpn+{ρdev[σ(εn+1,dn,εpn)]for Φ=00else

and

(50)αn+1=αn+ρ(r3+H^αn)

where we make use of the discretized stress

(51)σ(εn+1,dn,εpn)=f(dn)E0:(εn+1εpn).

Taking the derivative of the updated plastic strains with respect to the current total strains leads to the equation

(52)εpn+1εn+1=ρεn+1dev[σ(εn+1,dn,εpn)]+ρdevσ(εn+1,dn,εpn)εn+1=ρεn+1dev[σ(εn+1,dn,εpn)]+ρdev[f(dn)E0].

As the plastic part is based on a rate-independent ansatz, it is not possible to directly achieve the necessary derivative εpn+1/εn+1 owing to the dependence on the consistency parameter ρ, or rather its derivative with respect to the strains ρ/εn+1. To solve this equation, investigations of the yield function Φ=0 are required

(53)Φ=|dev[σ(εn+1,dn,εpn+1)]|(r3+H^αn+1)=0

for which holds Φ=Φ(εn+1,ρ(εn+1)) since εpn+1=εpn+1(εn+1,ρ(εn+1)). Thus, the consistency parameter ρ is given implicitly by Equation (53). Both sides of Equation (53) may be differentiated with respect to the current total strains, which yields

(54)dΦdεn+1=Φεn+1+Φρρεn+1=0.

Consequently, the missing derivative of ρ can be found by

(55)ρεn+1=(Φρ)1Φεn+1.

To solve this equation, two further derivatives must be determined: Φ/ρ and Φ/εn+1. The first one, Φ/ρ, is calculated as follows:

(56)Φρ=dev[σ(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]ρH^αn+1ρ.

Hereto, the middle part

(57)dev[pεp(εn+1,dn,εpn+1)]ρ=dev[f(dn)E0:dev[σ(εn+1,dn,εpn)]]

and the last part is calculated separately to be

(58)αn+1ρ=r3+H^αn+1

yielding

(59)Φρ=dev[pεp(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]|dev[f(dn)E0dev[pεp(εn+1,dn,εpn)]]r3+H^αn+1.

The second one, Φ/εn+1, is calculated as follows:

(60)Φεn+1=dev[σ(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]εn+1

with

(61)dev[σ(εn+1,dn,εpn+1)]εn+1=dev[f(dn)E0f(dn)E0ρdev[f(dn)E0]].

Finally, the derivative of the yield function with respect to the current trains reads

(62)Φεn+1=dev[σ(εn+1,dn,εpn+1)]|dev[σ(εn+1,dn,εpn+1)]|dev[f(dn)E0f(dn)E0ρdev[f(dn)E0]]

which, in combination with Equations (59) and (55), closes the necessary equations for the computation of the tangent operator.

4 Numerical results

This section provides numerical examples that are used to investigate the coupled damage-plasticity model. The main focus in this work concentrates on the influence of plasticity; investigations regarding the damage model are given in [1], [2].

We present two different loading types, tension (Section 4.1) and shear (Section 4.2), and use four different finite-element meshes with 1024, 3136, 6400 and 10,816 elements, respectively, see Figure 2. A linear refinement is employed such that the coarsest and the finest mesh differ by approximately one decade, which allows to evaluate the quality of the mesh-objectivity. Furthermore, numerical results with different strain rates are presented in order to show the influences of the rate-dependent formulation of the damage evolution (Section 4.3).

Figure 2: 2D: Meshes with 1024, 3136, 6400 and 10,816 elements (from left to right), inclusion highlighted in dark gray.
Figure 2:

2D: Meshes with 1024, 3136, 6400 and 10,816 elements (from left to right), inclusion highlighted in dark gray.

In all cases, periodic boundary conditions with fixed corner points are used and the prescribed strains are applied homogeneously. The presented specimen consists of two materials: surrounding the centered round inclusion (“inc”) (dark gray) there is a matrix (“mtx”) material (light gray) that is able to undergo damage and evolve plastic strains, whereas no damage or plastic effects are occurring in the inclusion because it behaves purely elastically. Furthermore, the matrix material is much softer than the inclusion material and has a slightly larger Poisson’s ratio. The exact material parameters (i.e. [34]) are shown in Table 1. Here, the viscosity r2mtx is paramount because of the applied viscous regularization. The chosen viscosity (2 × 106 MPa s) is in a realistic order of magnitude; it is in a range of the solids bitumen (up to 105 MPa s) and asphalt (up to 107 MPa s) and much smaller than for granite (about 1013 MPa s) or glass (up to 1015 MPa s). As demonstrated by the numerical results in the following, this physically sound value for the viscosity allows for the computation of mesh-independent results even for rather slow loading velocities. Consequently, no numerically motivated viscosity has to be chosen, which is usually unrealistically high.

Table 1:

Parameters for the respective model.

Einc104.5 GPar2mtx2 × 106 MPa s
νinc0.31 [–]r1mtx5 GPa
Emtx50.0 GPaαmtx1 × 10−5 [–]
νmtx0.36 [–]βmtx2.2 [–]
κmtx1 × 10−7 [–]r3mtx20 GPa
Δtmtx1 sH^mtx0.03 × Emtx

4.1 Tension test

The first numerical example is a tension test in the horizontal direction with a loading rate of ε˙11tension=5×1061/s (all other components of averaged ε are zero) that is applied in 270 time steps. The global answer is plotted in terms of the stresses, which are averaged over the total geometry, see Figure 3. The coarsest mesh with 1024 elements is represented by the blue line, followed by the finer meshes with 3136 and 6400 elements in green and red and finally the finest mesh with 10,816 elements is represented by the black line. Specific time steps at which several contour plots for internal variables are shown subsequently, are highlighted in red.

Figure 3: Tension: global stress response for 1024, 3136, 6400 and 10,816 elements (blue, green, red and black), highlighted in red at steps 200, 245, 255 and 270 (from left to right).
Figure 3:

Tension: global stress response for 1024, 3136, 6400 and 10,816 elements (blue, green, red and black), highlighted in red at steps 200, 245, 255 and 270 (from left to right).

It is obvious that the results of all meshes are quasi-identical and the global answer of all geometries provide mesh-independent results. The plot also shows the three characteristics of the model: firstly, there is an elastically evolving part represented by a linear curve. Afterwards, there is a kink in the curve, which then continues linearly, this is where the plastic strains start to evolve and isotropic hardening is taking place. And finally, there is a very strong drop in the stresses, which is caused by the evolving damage. The smooth characteristic of the curve along with the good convergence of the simulations show that the use of an explicit Euler integration scheme is justified.

The results in Figure 4Figure 6 compare the damage variable d, damage function f(d), and the norm of the plastic strains |εp| of the different meshes for the case of tension loading.

Figure 4: Tension: distribution of damage $d$d for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).
Figure 4:

Tension: distribution of damage d for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).

Figure 5: Tension: distribution of damage function $f(d)$f(d) for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).
Figure 5:

Tension: distribution of damage function f(d) for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).

Figure 6: Tension: distribution of plastic strains $\boldsymbol{\varepsilon}_\text{p}$εp for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).
Figure 6:

Tension: distribution of plastic strains εp for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 200, 245, 255 and 270 (downwards).

Figure 4 shows the evolution of damage during loading steps 200, 245, 255 and 270 (downwards). The scale goes from 0 (blue, no damage) to 4.5 × 10−3 (red, damage). The results show damage occurring very locally on the left and right side of the inclusion within the matrix material. This narrow damaged area describes a localized material failure initiating and propagating at the boundary of the inclusion. However, the local loss of stiffness, described by the damage function f(d), shows a much more smeared distribution, see the following paragraph. All meshes provide almost the same distributions and similar absolute values. Only the coarse mesh seems to be too inaccurate for a reasonable resolution of the distribution of the damage variable. Naturally, with a finer mesh the localization can be characterized in greater detail, but absolute values of the damage counterbalance its refinement so that the global response stays the same and the results are therefore mesh-independent. This is confirmed by the smooth transitions between undamaged and damaged material independently of the element borders.

The results in Figure 5 show the damage function and its evolution during the same loading steps. The scale goes from 1 (blue, intact material) to 0 (red, failed material). It is obvious that even already very small values of d lead to a strong loss of stiffness, which is described by the damage function. That is why the distribution of failed material is significantly larger than one would expect from the distribution of damage. As these areas of failed material now possess almost no stiffness, we can assume the properties of ultimately damaged material here. In addition, the main direction of propagation is oriented horizontally, while there is also a smaller propagation in the diagonal directions. Again, the results are mesh-independent and there are very nice smooth transition zones between intact and failed material.

The plastic strains are finally presented in Figure 6. At this point, the same loading steps (200, 245, 255 and 270) are chosen. As the yield stress constitutes a lower energetic barrier for the evolution of plastic effects than the energetic damage barrier r1, the evolution of εp starts earlier. The scale for these results goes from 0 (blue, no plastic strains) to 4.3 × 10−3 (red, plastic strains). The distribution of the plastic strains is a little more homogeneous than the distribution of the damage and is not only limited to the edges of the inclusion. The plastic strains are evolving much more into the matrix material, especially in the diagonal directions, while there is almost no evolution in the horizontal direction where most of the damage is taking place. It is also interesting that the high plastic strains in diagonal directions obviously have an effect on the evolution of damage: while damage is mainly evolving in the horizontal direction, the evolution of high plastic strains in diagonal directions consequently leads to a reduction in the stiffness in the same area. Thus, the plasticity influences the damage behavior. Once again, the results show a mesh-independent behavior with smooth transition areas.

4.2 Shear test

The second numerical example is a shear test in the x/y direction with a loading rate of ε˙12shear=5×1051/s (all other components of averaged ε are zero) that is applied in 290 time steps. Again, the global answer is plotted against the stresses, which are averaged over the total geometry, see Figure 7. Also again, the meshes with 1024, 3136, 6400 and 10,816 elements are represented by the blue, green, red and black line, while specific time steps, at which several contour plots for internal variables are shown subsequently, are highlighted in red.

Figure 7: Shear: global stress response for 1024, 3136, 6400 and 10,816 elements (blue, green, red and black), highlighted in red at steps 150, 255, 270 and 290 (from left to right).
Figure 7:

Shear: global stress response for 1024, 3136, 6400 and 10,816 elements (blue, green, red and black), highlighted in red at steps 150, 255, 270 and 290 (from left to right).

In the case of shear loading, the results of all meshes are quasi-identical as well and thus the global answer of all geometries again provides mesh-independent results. The characteristic of the curve is a little different from the case of tension loading: the middle part, where plastic strains are evolving and isotropic hardening is taking place, shows a lower slope as in the tension case. The drop is also slightly weaker at this point (less steep descent as in tension), but besides that, the characteristic is basically the same: it starts with an elastic part, followed by the plastic part, and finally the drop due to the damage.

Likewise, the results in Figure 8Figure 10 compare the damage variable d, damage function f(d), and the norm of the plastic strains εp of the different meshes now for the case of shear loading.

Figure 8: Shear: distribution of damage $d$d for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270, and 290 (downwards).
Figure 8:

Shear: distribution of damage d for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270, and 290 (downwards).

Figure 9: Shear: distribution of damage function $f(d)$f(d) for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270, and 290 (downwards).
Figure 9:

Shear: distribution of damage function f(d) for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270, and 290 (downwards).

Figure 10: Shear: distribution of plastic strains $\boldsymbol{\varepsilon}_\text{p}$εp for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270 and 290 (downwards).
Figure 10:

Shear: distribution of plastic strains εp for 1024, 3136, 6400 and 10,816 elements (from left to right) for steps 150, 255, 270 and 290 (downwards).

Firstly, the damage evolution during loading steps 150, 255, 270 and 290 (downwards) is presented in Figure 8. In this case the scale goes from 0 (blue, no damage) to 4.0 × 10−3 (red, damage), which is only slightly smaller than in tension. Obviously, there is not as much damage evolved as in the tension case, but instead it is not only arising at the left and right side, but also at the top and bottom of the inclusion, naturally within the matrix material. Again, the damage is occurring very locally and is now propagating from the described locations until it completely circles the inclusion. Besides the damage around the inclusion, there are very homogeneous damaged areas in diagonal directions, which are a result of homogeneous stress conditions due to plastic strains shown afterwards. Again, the material failure must be proven by the damage function. Moreover, mesh-independent results are achieved again; all meshes provide the same distributions and almost the same absolute values and exhibit smooth transitions between undamaged and damaged material.

Figure 9 shows the evolution of the damage function during the same loading steps. The scale goes from 1 (blue, intact material) to 0 (red, failed material). Basically, the same expectations arise as in the case of tension loading. The results are mesh-independent. Furthermore, the damage function illustrates a total loss of stiffness completely around the inclusion, which is corresponding to material failure initiating at the sides as well as at the top and bottom of the inclusion. These failure zones are propagating until the failed material is separating the matrix material from the inclusion. Consequently, both areas are no longer connected and can no longer transfer displacements or cause stresses.

The results in Figure 10 finally present the plastic strains. The same loading steps (150, 255, 270, and 290) are shown. Again, the evolution of plastic deformations initiates earlier than the evolution of damage. The scale now goes from 0 (blue, no plastic strains) to 2.5 × 10−2 (red, plastic strains). In the case of shear loading, the distribution of plastic strains is even more homogeneous and propagated than in the previous tension case and the amount of plastic strains is also significantly higher. Essentially, all matrix material is influenced by plastic strains, especially a specific zone around the inclusion with a concentration in the diagonal directions and even higher values in smaller regions next to it. Also here, the plasticity considerably influences damage evolution: after high plastic strains are accomplished in diagonal directions, a reduction in stiffness occurs in the same area. Again, the results are mesh-independent and have smooth transitions.

4.3 Dependence on the strain rates

As the damage part of this coupled material model is rate-dependent, we present for a full analysis the dependence of the material behavior on the applied strain rates in this section. Hereto, different strain rates were applied to the tension as well as the shear loading calculations and were compared by the global stress responses.

For a better comparison the same maximum strain loading was applied to all calculations; the maximum applied strain for the tension case was εmax=0.00135 and for the shear case it was εmax=0.0145. The results from the tension loading case are shown in Figure 11 with strain rates increasing from 0.3 to 0.5 and finally to 0.7 × 10−5 1/s (black, blue and green curves) in correspondingly 450, 270 and 193 time steps. Accordingly, the results from the shear loading case are shown subsequently in Figure 12 with strain rates from 3 to 5 and finally to 7 × 10−5 1/s (black, blue and green curves) in correspondingly 484, 290 and 208 time steps. In both loading cases, the calculations were performed with the previously mentioned meshes consisting of 1024, 3136, 6400 and 10,816 elements. Basically, only the coarsest mesh with 1024 elements differs slightly in the tension loading case and significantly in the shear loading case, for this reason the related curve is presented in a dashed way.

Figure 11: Tension: global stress response for a fast, medium and slow strain rate (green, blue and black), applied to 1024 (dashed), 3136, 6400 and 10,816 elements.
Figure 11:

Tension: global stress response for a fast, medium and slow strain rate (green, blue and black), applied to 1024 (dashed), 3136, 6400 and 10,816 elements.

Figure 12: Shear: global stress response for a fast, medium and slow strain rate (green, blue and black), applied to 1024 (dashed), 3136, 6400 and 10,816 elements.
Figure 12:

Shear: global stress response for a fast, medium and slow strain rate (green, blue and black), applied to 1024 (dashed), 3136, 6400 and 10,816 elements.

Figure 11 as well as Figure 12 basically show the same influences due to the different strain rates: the faster the loading is applied the later the drop of stresses occurs, i.e. the material strength is higher at higher loading velocities. This is caused by the viscous effects of the damage formulation, which delay the evolution processes. In contrast, the slower the loading is applied the lower the material strength is and the drop of stresses occurs earlier. Consequently, there is a direct influence of the strain rate to the onset of damage evolution. Interestingly, the curves run in parallel during damage evolution. Thus, the different strain rates rather affect the initialization of damage than its evolution.

In addition, the presentation of different meshes also allows to investigate the influence of the loading velocity on the mesh-objectivity. In the case of tension loading presented in Figure 11, the strain rates do not influence the mesh-objectivity anyhow. It can be assumed that the coarsest mesh is sufficiently converged and the strain rate is fast enough that viscous effects regularize the mathematical problem. In contrast, the case of shear loading presented in Figure 12, exhibits a different behavior regarding the coarsest mesh (dashed curve). Obviously, the combination of a very coarse mesh with a slower strain rate forces the model to its limit: the results are no longer mesh-independent, which is well known for damage models based on viscous regularization. However, this only applies to the very coarse mesh as the meshes with 3136, 6400 and 10,816 elements still have a very good agreement. For the two faster strain rates, all meshes supply convincing mesh-independent results.

5 Conclusions and outlook

This paper presented the plasticity enhancement of the new damage model published in [1], [2]. The result is a more comprehensive damage model that is characterized by a relaxation-based approach without the need for complex regularization procedures besides considering viscous effects. The efficient model provides mesh-independent results and not only softening, but also irreversible as well as hardening effects can be described. In particular, the model provides mesh-independent numerical results for physically reasonable viscosity.

The detailed implementation and enhancement with plasticity was presented in Section 2 and 3. The numerical examples were presented in Section 4, in which two loading types acting on a representative volume element consisting of a damaging matrix material and an undamaging, stiff inclusion material were applied. The results were always mesh-independent for four different finite-element meshes, which were discretized with element sizes differing by a factor of 10. Smooth transitions between undamaged and damaged material were observed, independently of any element borders. Damage evolved rather locally and the evolution of plastic strains was in agreement with practical experience from material science: firstly, plastic effects are occurring and then damage starts to evolve as loading increases. These characteristics were especially evident in the global stress response diagrams. For a complete investigation, several loading rates were applied. The numerical results demonstrated that rather the onset of damage than its evolution is influenced by the respective strain rate and mesh-independence remains for most cases.

Altogether, the numerical results underline the functionality and efficiency of the new coupled damage-plasticity model. The results are physically sound and mesh-independent. In the future, a comparison with experimental data for verification purposes as well as further modifications or enhancements are possible.

References

[1] Junker P, Schwarz S, Makowski J, Hackl K. Continuum Mech. Therm. 2017, 29 (1), 291–310.10.1007/s00161-016-0528-8Suche in Google Scholar

[2] Schwarz S, Junker P, Hackl K. Proc. Appl. Math. Mech. 2016, 16 (1), 173–174.10.1002/pamm.201610075Suche in Google Scholar

[3] Dimitrijevic BJ, Hackl K. Tech. Mech. 2008, 28 (1), 43–52.Suche in Google Scholar

[4] Dimitrijevic BJ, Hackl K. Int. J. Numer. Method. Biomed. Eng. 2011, 27 (8), 1199–1210.10.1002/cnm.1350Suche in Google Scholar

[5] Carstensen C, Hackl K, Mielke A. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2002, 458 (2018), 299–317.10.1098/rspa.2001.0864Suche in Google Scholar

[6] Hackl K, Fischer FD. Proc. R. Soc. A Math. Phys. Eng. Sci. 2008, 464 (2089), 117–132.Suche in Google Scholar

[7] Junker P, Jerzy M, Hackl K, Makowski J, Hackl K. Continuum Mech. Therm. 2014, 26 (3), 259–268.10.1007/s00161-013-0299-4Suche in Google Scholar

[8] Bazant ZP, Jirásek M. J. Eng. Mech. 2002, 128 (11), 1119–1149.10.1061/(ASCE)0733-9399(2002)128:11(1119)Suche in Google Scholar

[9] Nguyen GD, Houlsby GT. Int. J. Numer. Anal. Methods Geomech. 2008, 32, 39–413.Suche in Google Scholar

[10] Grassl P, Jirásek M. Int. J. Numer. Anal. Methods Geomech. 2006, 30, 71–90.10.1002/nag.479Suche in Google Scholar

[11] Sciarra FM. Int. J. Plasticity 2012, 34, 114–138.10.1016/j.ijplas.2012.01.009Suche in Google Scholar

[12] Nguyen GD, Korsunsky AM, Belnoue JPH. Int. J. Plasticity 2015, 64, 56–75.10.1016/j.ijplas.2014.08.001Suche in Google Scholar

[13] Nedjar B. Int. J. Solids Struct. 2001, 38, 5421–5451.10.1016/S0020-7683(00)00358-9Suche in Google Scholar

[14] Comi C, Perego U. Int. J. Numer. Method Eng. 1996, 39, 3731–3755.10.1002/(SICI)1097-0207(19961115)39:21<3731::AID-NME24>3.0.CO;2-ZSuche in Google Scholar

[15] Poh LH, Swaddiwudhipong S. Int. J. Solids Struct. 2009, 46, 4369–4378.10.1016/j.ijsolstr.2009.08.025Suche in Google Scholar

[16] Engelen RAB, Geers MGD, Baaijens FPT. Int. J. Plasticity 2003, 19, 403–433.10.1016/S0749-6419(01)00042-0Suche in Google Scholar

[17] Svedberg T, Runesson K. Int. J. Solids Struct. 2000, 37, 7481–7499.10.1016/S0020-7683(00)00208-0Suche in Google Scholar

[18] Borst R, Pamin J, Geers MGD. Eur. J. Mech. A/Solids 1999, 18, 939–962.10.1016/S0997-7538(99)00114-XSuche in Google Scholar

[19] Borst R, Pamin J, Peerlings RHJ, Sluys LJ. Comput. Mech. 1995, 17, 130–141.10.1007/BF00356485Suche in Google Scholar

[20] Ren X, Zeng S, Li J. Comput. Mech. 2015, 55, 267–285.10.1007/s00466-014-1100-7Suche in Google Scholar

[21] Ibrahimbegovic A, Jehel P, Davenne L. Comput. Mech. 2008, 42, 1–11.10.1007/s00466-007-0230-6Suche in Google Scholar

[22] Javani HR, Peerlings RHJ, Geers MGD. Comput. Mech. 2014, 53, 625–639.10.1007/s00466-013-0922-zSuche in Google Scholar

[23] Mediavilla J, Peerlings RHJ, Geers MGD. Comput. Methods Appl. Mech. Eng. 2006, 195, 4617–4634.10.1016/j.cma.2005.10.001Suche in Google Scholar

[24] Vrech SM, Etse G. Comput. Methods Appl. Mech. Eng. 2009, 199, 136–147.10.1016/j.cma.2009.09.025Suche in Google Scholar

[25] Addessi D, Marfia S, Sacco E. Comput. Methods Appl. Mech. Eng. 2002, 191, 1291–1310.10.1016/S0045-7825(01)00325-5Suche in Google Scholar

[26] Needleman A. Comput. Methods Appl. Mech. Eng. 1988, 67 (1), 69–85.10.1016/0045-7825(88)90069-2Suche in Google Scholar

[27] Sluys LJ, de Borst R. Int. J. Solids Struct. 1992, 29 (23), 2945–2958.10.1016/0020-7683(92)90151-ISuche in Google Scholar

[28] Dube JF, Cabot GP, La Borderie C. J. Eng. Mech. 1996, 122 (10), 939–947.10.1061/(ASCE)0733-9399(1996)122:10(939)Suche in Google Scholar

[29] Niazi MS, Wisselink HH, Meinders T. Comput. Mech. 2013, 51 (2), 203–216.10.1007/s00466-012-0717-7Suche in Google Scholar

[30] Wang WM, Sluys LJ, de Borst R. Eur. J. Mech. A 1996, 15 (3), 447–464.Suche in Google Scholar

[31] Chaboche JL, Feyel F, Monerie Y. Int. J. Solids Struct. 2001, 38 (18), 3127–3160.10.1016/S0020-7683(00)00053-6Suche in Google Scholar

[32] Suffis A, Lubrecht TAA, Combescure A. Int. J. Solids Struct. 2003, 40 (13–14), 3463–3476.10.1016/S0020-7683(03)00153-7Suche in Google Scholar

[33] Langenfeld K, Junker P, Mosler J. Continuum Mech. Therm. 2018, 30 (5), 1125–1144.10.1007/s00161-018-0669-zSuche in Google Scholar

[34] Shankar S, Mayuram MM. Int. J. Solids Struct. 2008, 45, 3009–3020.10.1016/j.ijsolstr.2008.01.017Suche in Google Scholar

Published Online: 2018-12-14

©2018 Walter de Gruyter GmbH, Berlin/Boston

This article is distributed under the terms of the Creative Commons Attribution Non-Commercial License, which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Heruntergeladen am 19.9.2025 von https://www.degruyterbrill.com/document/doi/10.1515/jmbm-2018-2001/html
Button zum nach oben scrollen