Optimization of fiber geometry for fiber reinforced composites considering damage

Optimization of fiber geometry for fiber reinforced composites considering damage

ARTICLE IN PRESS Finite Elements in Analysis and Design 46 (2010) 401–415 Contents lists available at ScienceDirect Finite Elements in Analysis and ...

966KB Sizes 5 Downloads 140 Views

ARTICLE IN PRESS Finite Elements in Analysis and Design 46 (2010) 401–415

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Optimization of fiber geometry for fiber reinforced composites considering damage J. Kato , E. Ramm Institute of Structural Mechanics, University of Stuttgart, Pfaffenwaldring 7, D-70550 Stuttgart, Germany

a r t i c l e in fo

abstract

Article history: Received 19 May 2009 Accepted 4 January 2010 Available online 18 February 2010

The present contribution deals with an optimization strategy of fiber reinforced composites. Although the methodical concept is very general we concentrate on Fiber Reinforced Concrete with a complex failure mechanism resulting from material brittleness of both constituents matrix and fibers. Because of these unfavorable characteristics the interface between fiber and matrix plays a particularly important role in the structural response. A prominent objective for this kind of composite is the improvement of ductility. The influential factors on the entire structural response of this composite are (i) material parameters involved in the interface, (ii) the material layout at the small scale level, and (iii) the fiber geometry on the macroscopic structural level. The purpose of the present paper is to improve the structural ductility of the fiber reinforced composites applying an optimization method with respect to the geometrical layout of continuous long textile fibers. The method proposed is achieved by applying a so-called embedded reinforcement formulation. This methodology is extended to a damage formulation in order to represent a realistic structural behavior. For the optimization problem a gradient-based optimization scheme is assumed. An optimality criteria method is applied because of its numerically high efficiency and robustness. The performance of the method is demonstrated by a series of numerical examples; it is verified that the ductility can be substantially improved. & 2010 Elsevier B.V. All rights reserved.

Keywords: Shape optimization Damage Fiber reinforced composites

1. Introduction 1.1. Overview The most widely used light-weight fiber reinforced composites are fiber reinforced polymers (FRP) where often long glass, carbon or aramid fibers are placed into a polymer matrix leading to a good-natured ductile material. The present study addresses a new composite material, namely Fiber Reinforced Concrete (FRC), sometimes also called Textile Reinforced Concrete. It differs from FRP that the fibers are placed in a fine grained concrete or mortar matrix (Fig. 1), often as a reinforcement mesh with a relatively low fiber content, making it economically attractive. Unlike conventional steel reinforcement, this kind of textile fiber is corrosion free; this property allows to manufacture light-weight thin-walled composite structures. However, FRC structures show very complex failure mechanisms resulting on the one hand from the material brittleness of

 Corresponding author. Tel.: + 49 711 685 66123; fax: +49 711 685 66130.

E-mail addresses: [email protected], [email protected] (J. Kato). 0168-874X/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2010.01.001

both constituents, fibers and matrix, and on the other hand from their interface behavior introducing the necessary ductility. The specific characteristic of FRC is an ideal target for optimization applying the overall structural ductility as objective which ought to be maximized for a prescribed fiber volume. In this context the ‘structural ductility’ means ‘energy absorption capacity’ which is measured by the internal energy summed over the entire structure up to a prescribed displacement of a dominant control point. For this task geometrical parameters on the small scale level like fiber location, orientation, size, length, spacing as well as combinations of different fiber materials can be considered as the representative design variables, see Kato et al. [15]. As an example for FRC Karihaloo and Lange-Kornbak [14] introduce a multi-objective optimization problem to maximize tensile strength, deformability and toughness of high performance fiber reinforced concrete mixes (HPFRC) in which chopped ‘short’ fibers are utilized. In their investigation the fiber diameter, length, fracture toughness and fiber volume fraction are chosen as the design variables and an analytical solution instead of the finite element method is applied for the optimization problem. Kato et al. [16] introduce a multiphase material optimization to improve the ductility of FRC with respect to fiber size, length and the combination of different fiber materials using the finite

ARTICLE IN PRESS 402

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

Fig. 1. FRC structures (Refs. [6,13]).

element method for the discretization. In that study the concept of ‘volume fraction’ of fiber material(s) is applied describing the design variables mentioned. This approach is considered as a material distribution problem derived from conventional topology optimization. An optimal fiber orientation of composites has been investigated in several contributions. Most of them focus on the optimal fiber orientation using laminated FRP structures. In those cases the fiber angle in individual plies is chosen as design variable, see Stegmann and Lund [34] and Stolpe and Stegmann [35], just to mention only a few references. The design variables in both mentioned approaches, i.e. optimizing material or fiber angle distribution, are defined locally within certain plies or even in a patch of finite element; this restriction limits to achieve final optimal fiber layouts. For example in the investigation of Kato et al. [16] fiber materials are defined only in prescribed design elements, furthermore only straight fibers are allowed. In Stegmann and Lund [34] and Stolpe and Stegmann [35] the fibers are discontinuous between adjacent elements, leading to a so-called discrete fiber distribution. The purpose of the present study is to introduce a methodology improving the structural ductility of FRC with respect to the fiber geometry which is independent of the fixed finite element mesh. The mechanical model of FRC is briefly summarized in the next section. The optimization problem is solved by a gradientbased optimization scheme. An optimality criteria method, see Patnaik et al. [28], is used because of its numerically high efficiency and robustness. For the sensitivity analysis a variational semi-analytical direct method is used. The present method is not restricted to FRC but may as well be applied to other fiber reinforced composites, for example fiber reinforced glass (FRG).

1.2. Present numerical model We introduce so-called embedded reinforcement elements where the fiber geometry is globally defined. These embedded reinforcement elements have been originally introduced by Phillips and Zienkiewicz [32]. Chang et al. [8] modified the concept allowing for straight reinforcement segments to be placed at any angle with respect to the local axes of isoparametric concrete elements. Balakrishnan and Murray [1] introduce an embedded formulation with bond–slip relation between concrete and reinforcement. Further improvements by Elwi and Hrudey [10] allow for a general curved reinforcement formulation in the embedded element. Hofstetter and Mang [11] apply the embedded reinforcement formulation for a thin-walled prestress concrete shell structure where geometry of curved tendons are introduced by an analytical expression. The extension to a threedimensional formulation is discussed by Barzegar and Maddipudi [2]. Recently Huber [12] applies the bond–slip relation by

Balakrishnan and Murray [1] for a 3D model with straight reinforcement bars considering nonlinear material models for both concrete and steel reinforcement. In the present study we apply the embedded reinforcement element including the bond–slip relation by Balakrishnan and Murray [1] for a two-dimensional model in which a curved fiber geometry is allowed. The fiber geometry is defined globally by Be´zier-splines. The materials for both concrete and fibers are modeled by a gradient-enhanced damage formulation, see Peerlings et al. [29,30], Peerlings [31]. For the interface between ¨ concrete and fiber a discrete bond model is applied, see Kruger et al. [18–20], Xu et al. [37]. The assumption introduced by Balakrishnan and Murray [1] is adopted for the interfacial kinematical relation. 1.3. Structure of paper The paper is organized as follows. Firstly the applied material models are briefly presented in Section 2. Secondly the concept of the globally defined fiber geometry is introduced in Section 3. Then the embedded reinforcement element is described in Section 4. The finite element formulation of FRC is introduced in Section 5, which is composed of three individual material formulations, namely the gradient enhanced damage for both concrete and fibers and the interface model. Some details are shifted to the Appendix, e.g. the transformation matrices relevant for the fiber orientations and the linearization of the model. Finally the optimization problems and the derivation of the sensitivity analysis are discussed in Sections 6 and 7. Note that in the present paper the following superscripts ðÞc , ðÞf , and ðÞi denote the terms for ‘concrete’, ‘fiber’, and ‘interface, respectively. In some cases we utilize a compact expression, e.g. ðÞc þ f ¼ ðÞc þðÞf . A subscript ðÞL or ðÞG indicates that the value ðÞ is measured in the local or global coordinate system, respectively. However the notion ðÞG is introduced only when we need to emphasize it, otherwise we skip it for simplicity.

2. Applied material models 2.1. Material model for constituents concrete and fiber In this study the nonlinear material behavior of both concrete and fiber is described by an isotropic continuum damage model. Firstly an equivalent strain measure is defined. For the concrete matrix de Vree’s definition (de Vree et al. [36]) of equivalent strains ecv is adopted as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k1 1 ðk1Þ2 2 12k c ð1Þ I1 þ ev ðI1 ; J2 Þ ¼ I  J ; 2kð12nÞ 2k ð12nÞ2 1 ð1 þ nÞ2 2

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

403

where I1 denotes the first invariant of the strain tensor and J2 the second invariant of the deviatoric strain tensor. k indicates the ratio of compression relative to the tension strength and n is Poisson’s ratio. For the fiber we follow Mazars’s definition (Mazars and Pijaudier-Cabot [24]) since the fiber is assumed to be a onedimensional model qffiffiffiffiffiffiffiffiffiffiffiffiffiffi efv ¼ /efL S2 ; ð2Þ where efL is the fiber strain and /S denotes the Macauley bracket /xS ¼ ðx þjxjÞ=2. For the damage evolution of both concrete and fiber we use an exponential damage law introduced by Mazars and Pijaudier-Cabot [24] as DðkÞ ¼ 1

k0 ð1a þ aebðkk0 Þ Þ if k Z k0 ; k

ð3Þ

where D stands for the damage parameter ð0 r Dr 1Þ, a defines the final softening stage and b governs the rate of damage growth. k0 is a threshold variable which determines damage initiation and k represents the most severe deformation the material has experienced during loading. In a conventional local damage model, k is related to the local equivalent strain ev and the history variable k is defined by the Kuhn–Tucker relations, i.e. k_ Z0, ev k r 0, k_ ðev kÞ ¼ 0. For non-local damage, k is related to a weighted volume average of the local equivalent strain ev , denoted as non-local equivalent strain e~ v . In the gradientenhanced damage model (Peerlings et al. [29,30], Peerlings [31]) e~ v is approximated implicitly as follows:

e~ v cr2 e~ v ¼ ev ;

ð4Þ

2

where r denotes the Laplacean operator and c is a positive parameter of the dimension length squared regularizing the localization of the deformation. Thus in the Kuhn–Tucker equations ev is replaced by the non-local equivalent strain e~ v which is discretized in the finite element sense. Elastic unloading is included in the traditional way.

Fig. 2. Discrete bond model.

with 2

c ¼ 1 þ tanh4ar

sR 0:1f c

af nes 1

r2 ðr þ hÞ2

!1 3 5;

ð7Þ

where c denotes an additional parameter ð1 o c o 2Þ which considers the influence of the kind of fiber material, the loading condition and the stresses perpendicular to a fiber direction. sm;0 and sf;0 denote the initial adhesion strength and sliding friction strength, respectively. r describes a fiber (roving) radius, n is Poisson’s ratio of a fiber and h is the surface roughness of a fiber. ar and af are constants assuming the lateral deformation of a fiber. These properties are given depending on the kind of fiber material used. f c is the uniaxial compressive strength of concrete, es the uniaxial strain and sR defines the stress perpendicular to a fiber. For a detailed description of this model it is referred to ¨ Kruger et al. [18–20]. In this model loading and unloading conditions are also considered. This one-dimensional interface model is originally formulated for a fiber in a three-dimensional setting. If this model is utilized in a two-dimensional space, the interface has to be approximated to hold the original total interface area.

3. Global layout of fiber geometry 2.2. Material model for interface 3.1. Concept and design variables In this study nonlinear interfacial behavior between fiber and ¨ matrix is expressed by a discrete bond model, see Kruger et al. [18]. This model was obtained by experiments for different textile fiber materials and leads to a realistic interface response of FRC. The significant factors governing interfacial response are the bond strength and the debonding behavior. The influence of material properties at a small scale level and the stresses perpendicular to the fiber direction are included in the material formulation as important parameters. The bond stress–slip ðsiL uiL Þ relation is expressed as (  1=R ) 1  s0 for uiL r w1 ; siL ¼ w~  b þ ð1bÞ  ð5Þ ~R 1þw ~ ¼ uiL =w0 denotes the normalized slip. uiL is the slip length where w which will be introduced in the following section. w0 is a factor defined by the initial tangent k1 . k2 is the tangent at slip w1 (see Fig. 2) where the bond stress achieves the maximum bond strength. b ¼ k2 =k1 and s0 ¼ k1  w0 are parameters to calculate the stresses and R defines the radius of curvature at slip w1 . The stress–slip relation for the range uiL 4 w1 is simply described by the adhesion strength sm and the friction bond strength sf (see Fig. 2),

The geometry of a continuous long fiber is defined in the global coordinate system. One of the big advantages of FRC structures is that they do not need thick concrete covers. Furthermore hooks of textile fibers are not used. Due to these characteristics the layout of textile fibers in FRC can be rather simple often parallel fibers or a mesh of straight fibers are used, see Fig. 1. Slightly curved fibers are advantageous if an optimal structural response is looked for. We approximate the fiber geometry by Be´zier-splines, defined as parametric curves by control points. A quadratic Be´zier-spline and its mathematical formulation are introduced in Fig. 3, where r stands for a position vector of the spline; W ð0 r W r1Þ is the local coordinate system of the spline. pj indicates the j-th control point. Of course we could apply another parameterization allowing more general geometrical definitions such as a level set function (Fig. 4). The fiber approximated by Be´zier-splines is embedded in the structure and the control points are moved in order to obtain the optimal fiber layout. The entire domain of structure is defined in a parametric space s ð0 r s r1Þ, see Fig. 3. Thus the normalized coordinates of control points turn out to be the design variables defining the global fiber layout in the physical space. According to this the j-th position vector of control point pj can be expressed as follows

sm ¼ sm;0 c;

^ yÞ ^ þ ðsxj Lx ; syj Ly Þ; rj ðsxj ; syj Þ ¼ Oðx;

sf ¼ sf;0 c

ð6Þ

ð8Þ

ARTICLE IN PRESS 404

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

Fig. 3. (a) Quadratic Be´zier-spline and (b) concept of global layout of fiber geometry.

Fig. 4. Element patch describing intersections and related Newton–Raphson algorithm.

^ y^ where O stands for the coordinate origin of the structure; x, are the corresponding global coordinates of O. L denotes the contour lengths of the structure and the superscripts x, y on L as well as s indicate the direction. Inserting Eq. (8) into the general mathematical formulation of Be´zier-splines leads to the geometric formulation of a fiber including the design variables s as follows rðW; sx ; sy Þ ¼

nb X j¼0

Fj ðWÞrj ðsxj ; syj Þ with Fj ¼

nb ! Wj ð1WÞnb j ; ðnb jÞ!j!

fiber and fixed finite element mesh can be calculated. It is necessary to determine the global coordinates for intersections of fibers and mesh in order to establish the stiffness matrix and afterwards the internal forces of embedded fiber elements. This procedure is described in the sequel.

3.2. Determination of intersections between finite element mesh and fiber ð9Þ

where nb is the order of the Be´zier-spline. Note that the coefficients F are independent of the design variables s. Once the fiber geometry is defined by Eq. (9) the intersections between

The basic procedure determining the intersections between the mesh and a fiber is explained in Box 1. This procedure is continued until the end of the fiber. Note that the local coordinate system W and additional parameter $ have been introduced

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

405

only for the determination of intersections and are no longer necessary thereafter. Step 1: (1.1) start with control point p0 and obtain its position vector r0(W), where W =0 (1.2) calculate the descent along the fiber at p0, i.e. rW rðWÞ at W = 0 (1.3) find intersection of element boundary with descent line Step 2: (2.1) parameterize element boundary by normalized parameter $ð0 r $ r 1Þ using position vectors r^ a and r^ b of points a and b (2.2) determine global coordinates of intersection (x, y)p by Newton – Raphson scheme, see box in Fig 4, where R is the residual vector between points rðWÞ and rð$Þ Step 3: (3.1) if converged at (2.2), store the global coordinates of intersection p; ðx; yÞp ¼ rðW^ Þ (3.2) otherwise choose other element boundary of element and restart from (2.1) Step 4: (4.1) replace reference point p0 by the determined intersection p and repeat the same procedure from (1.2) with the new descent rW rðW^ Þ simplicity approximated by a straight line leading to a polygonal layout as already indicated in Fig. 5.

3.3. Inverse mapping for local coordinate of fiber In the previous section we introduced the procedure to determine the global coordinates of the intersections ðx; yÞp . It is also necessary to determine the associated parametric (natural) coordinates ðx; ZÞp of the corresponding finite element in order to perform the integration for the internal virtual work of the fiber. This so-called nonlinear inverse mapping is described in detail in Elwi and Hrudey [10] and Barzegar and Maddipudi [2]. For the general isoparametric mapping the global coordinate ðx; yÞ of an arbitrary point p in an element is expressed by the shape function N and the global coordinates of the element nodes ðx^ k ; y^ k Þ, " " # Nðx; ZÞ x ¼ y 0

0 Nðx; ZÞ

p

# " p

# x^ k ; y^ k

ð10Þ

where the expression ðÞ ^ emphasizes a known value. k is the number of nodes of an element. In case of the inverse mapping Eq. (10) is transformed to the following equation such that the residual vector function R vanishes Rðx; ZÞ ¼

"

" # x^ y^

 p

Nðx; ZÞ

0

0

Nðx; ZÞ

# " p

x^ k y^ k

#

  0 ¼ : 0

ð11Þ

Finally, this nonlinear equation is solved by a Newton– Raphson scheme as shown in the box of Fig. 5 obtaining the associated natural coordinates ðx; ZÞp of point p. Once these coordinates of the intersections are determined the integration of the embedded reinforcement can be performed. For further processing of the fiber mechanics the curved fiber is for

4. Embedded reinforcement element 4.1. Kinematical assumption for interface between concrete and fiber The embedded reinforcement element applied in this study considers a bond–slip relation between concrete and fiber, see Fig. 6(a). We introduce the kinematical relation of the interface based on the assumption by Balakrishnan and Murray [1]. In the kinematical assumption the slip at an arbitrary point is considered as the relative displacement between concrete and fiber measured along the axis of the fiber. The components of the displacements can be written as ufL ¼ ucL þ uiL ;

ð12Þ

uiL

where is the slip length or relative displacement introduced in Eq. (5). ufL and ucL are the displacements of fiber and concrete at the considered point, respectively, see Fig. 6(b). The slips of the original curved fiber between two adjacent elements have to be equal; this is not the case for the polygonal geometry assumed above. In order to satisfy the compatibility at least in an average sense the slip length uiL is projected onto the global x-axis d ¼ cosy  uiL -uiL ¼ td

with t ¼ ðcosyÞ1 ;

ð13Þ

where y is the angle between fiber axis and x-axis, see Fig. 6. Thus the compatibility of the slip-length is enforced for d. From uiL the

Fig. 5. Determination of natural coordinates of intersections.

ARTICLE IN PRESS 406

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

Fig. 6. (a) Embedded reinforcement element patch and (b) notion for displacements of slip.

local bond strain eiL is obtained which in turn leads to the local fiber strain efL

efL ¼ ecL þ eiL :

ð14Þ

|{z} Te1 ecG

e

Matrix T transforms the global strain eG of a two-dimensional continuum into the local one eL under plane stress condition (see Appendix A). Te1 represents the first row of Te extracting the local strain ecL in fiber direction from the global concrete strain ecG .

5. Finite element formulation of FRC

Since we apply a gradient enhanced damage model for both concrete and fibers and use a nonlinear interface model between fibers and matrix the virtual work is decomposed into

dW ¼ dWint dWext ¼ dWcint þ dWfint þ dWiint dWext ¼ 0;

ð15Þ

where dWcint , dWfint , and dWiint stand for the internal virtual work of concrete, fibers and interfaces, respectively, and dWext denotes the external virtual work. The gradient enhanced damage model leads to a two-field formulation ðO ¼ Oc [ Of Þ and its formulation at the actual time t þ 1 is Z Z Z dWu ðu; duÞ ¼ de : r dO du  b^ dO du  t^ dG ¼ 0; ð16Þ

dWe ðe~ v ; de~ v Þ ¼

O

Z

dre~ v  s dO þ O

G

Z

Of

dWfe ¼

Z O

f

ð18Þ

Of

dre~ fv;L tfL dOf þ

Z Of

de~ fv;L ðe~ fv;L efv;L ðefL ÞÞdOf :

ð19Þ

According to Eq. (14) the fiber strain efL can be decomposed into a contribution of the local concrete strain ecL and that of the interface eiL , as indicated in Eq. (18). The work of the latter part together with the virtual work inside the interface due to the slip length uiL Z duiL siL dOi Oi

5.1. Virtual work

O

stress and strain field in the fibers, Z Z dWfu;int ¼ defL sfL dOf ¼ ðdecL þ deiL ÞsfL dOf ;

de~ v ðe~ v ev ðeÞÞ dO ¼ 0

ð17Þ

O

where ‘time’ t does not mean the ‘real time’ but simply the ‘loading step number’ for a nonlinear static problem. du and de~ v are the virtual displacement and non-local equivalent strain fields, respectively. r is the stress tensor, b^ and t^ are the prescribed body and traction forces. Eq. (16) is the usual virtual work expression, whereas Eq. (17) defines the weak form of an additional equilibrium equation for the non-local equivalent strain where s ¼ cre~ v is a work equivalent stress vector (see Peerlings et al. [29,30], Peerlings [31]). Without loss of generality we omit the body force in this study. Both equations of the virtual work can be split up into the parts of concrete and fibers, dWcu=e , dWfu=e , the latter one being reduced to a one-dimensional expression referring to the local

defines the total virtual work of the interface slip: Z Z dWiint ¼ deiL sfL dOf þ duiL siL dOi ¼ 0 8 duiL : Of

Oi

ð20Þ

5.2. Discretization The virtual work expressions (16)/(17) and (18)/(19) together with (20) contain three independent variables, namely the displacement field in the concrete element u, the non-local equivalent strain e~ v and the slip length uiL which are discretized in the finite element sense. In the present study a two-dimensional 8-node quadratic plane stress element is applied for the concrete matrix. The non-local strain is discretized by bilinear shape functions within this element. Please note, that the interface slip is discretized as a 3-node quadratic one-dimensional element ðni ¼ 3Þ, whereas the non-local strain enhancement of the fibers are only linearly interpolated based on the two values at the fiber beginning and end obtained from the non-local strain values in the concrete element nc X



k

Nk d

or

u ¼ Nd;

ð21Þ

k¼1

e~ v ¼

ne X

~ k ek N

or

~ e~ v ¼ Ne;

ð22Þ

k¼1

uiL ¼

ni X k¼1

Nik ðuiL Þk ¼

ni X k¼1

Nik ðtdÞk ¼

ni X

Nkd

k

or

uiL ¼ Nd;

ð23Þ

k¼1

where d is a vector with eight nodal displacements, e with four nodal values and d contains three nodal slip values.

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

The two nodal values of the projected slip lengths are summed in the vector d d ¼ ½d

1

d

2

3

d T :

ð24Þ

N contain the shape function for the interface defined in the global coordinate system. Analogously the local bond strain eiL of Eq. (14) in one element can be expressed as

eiL ¼

ni X

Bik ðtdÞk ¼ tBi d ¼ Bd;

ð25Þ

k¼1

where Bi and B stand for the so-called B-operator matrices for the interface defined in local and global coordinate systems, respectively. The local fiber strain efL can be written according to Eq. (14)

efL ¼ Te1 ecG þ eiL ¼ Te1 Bf d þBd:

ð26Þ

407

e and d leads after assembly to the following stiffness expression 3n 3n 2 2 c f cþf 2 3n þ 1 Kc þ f Kde Kfdd þ f int;u f ext f 7 Dd 7 6 dd 6 int;u c f 7 6 Kc þ f Kc þ f 6 6 7 0 7 ¼ 6 f int;e þf int;e 7 ; ð32Þ 7 4 De 5 6 ed ee 5 5 4 4 f i i D d Kdd 0 Kdd f int;i |fflfflfflffl{zfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Du KT

R

where KT , Du and R stand for the tangential stiffness matrix, the incremental displacement/strain vector and the residual force vector, respectively. The superscripts n and n þ 1 on the matrix and vectors indicate the iteration number in the increment. For the derivation of the corresponding stiffness matrices in KT and the forces in R it is referred to Appendices B to D.

6. Structural optimization of FRC

Introducing Eqs. (21)–(23) into the virtual work expressions leads to

6.1. Optimization problem

dWu ¼ dWcu;int þ dWfu;int dWext

In general an optimization problem is defined by an objective fðsÞ, equality constraints hðsÞ and inequality constraints gðsÞ. In this study the objective is to maximize the structural ductility for a prescribed fiber volume. As the ductility is defined by the internal energy summed up over the entire structure with a prescribed nodal displacement u^ (Maute [22] and Maute et al. [23]), the mathematical formulation of the optimization problem of FRC can be written as follows:

8dd

2 ¼

n ele [ e¼1

3

6Z 6

7 Z Z 7 T T T Bc rc dOc þ Bf ðTe1 ÞT sfL dOf lt þ 1 Nc t0 dG 7 7 ¼ 0; f O O G 4|fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}5

ddT 6 6

c

c

f ext

f

f int;u

f int;u

ð27Þ

dWe ¼ dWce þ dWfe 2

8de

[

e¼1

c f int;e

nfele

[

e¼1

subject to hðsÞ ¼ 3 sL r si rsU

6Z 6

7 Z 7 f ~ f ÞT ðe~ f ef ÞdOf 7 ¼ 0; ðB~ ÞT ðTd1 ÞT tfL dOf þ ðN v;L v;L 7 f f O O 4|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 5

deT 6 6

f f int;e

ð28Þ 3

2

7 6Z Z [ T6 T f T i f i7 7¼0 dWiint ¼ dd 6 B s d O þ N s d O L L 7 6 f O Oi 5 4|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} e¼1

Oc e^

rc dec dOc þ c

Z Z

sfL defL dOf þ f

Of e^ L

Z Z

#

siL duiL dOi i

Oi u^ L

ð33Þ

7 Z 7 c ~ c ÞT ðe~ c ec ÞdOc 7 ðB~ ÞT sc dOc þ ðN v v 7 c c O O 5 4|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

deT 6 6 2

þ

3

6Z 6

ncele

¼

"Z Z

minimize fðsÞ ¼ 

Z Of

ltf dOf V^ ¼ 0

i ¼ 1; . . . ; ns

ð34Þ

ð35Þ

where V^ denotes the prescribed fiber volume, sL and sU the lower and upper bounds of the design variables, and ns the number of design variables. tf represents the thickness of a fiber, which is assumed to be constant along the entire fiber. l is the length of a single fiber within an embedded reinforcement element and depends on the design variables si .

niele

8dd:

ð29Þ

6.2. Equilibrium conditions and total derivative of design function

f iint;i

c

Bc is the usual kinematic operator matrix; B~ is derived from the gradient of the non-local equivalent strain c

re~ cv ¼ B~ e;

ð30Þ

f

and B~ that of the corresponding part in the fiber f

re~ fv;L ¼ Td1 re~ fv;G ¼ Td1 B~ e;

ð31Þ

where Td1 is the first row of a rotation matrix Td , see Appendix A. l inserted in Eq. (27) denotes the load factor with respect to a reference traction load t0 .

The sensitivities of the design functions (objective, constraints etc.) depend on the gradients of the state variables uðd; e; dÞ. These are derived from the three equilibrium conditions (16), (17), and (20) at position n þ 1. The total derivative of the design functions with respect to the design variables can be decomposed into an explicit and an implicit part. As indicated the design functions depend on the structural response which in turn is implicitly related to the optimization variables, for example the objective f ¼ f ðs; uðd; e; dÞÞ. This leads to ex rs ðÞ ¼ rex s ðÞ þ ru ðÞr s u ¼ r s ðÞ þ r d ðÞrs d þ re ðÞr s e þ rd ðÞr s d;

ð36Þ 5.3. Element matrices Introducing damage and interface models into the virtual work expressions and linearizing with respect to the primary variables d,

ex s ðÞ

where r describes the explicit derivative with respect to the design variables. An optimality criteria method (see Patnaik et al. [28]) is applied to solve the optimization problem. For the sensitivity analysis a variational semi-analytical direct method is adopted and described in the sequel.

ARTICLE IN PRESS 408

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415 f

7.1. Overview of sensitivity analysis The main effort of sensitivity analysis is the calculation of implicit part rs u. For a direct sensitivity analysis this part is obtained by exploiting the stiffness expression containing the tangent stiffness matrix and the so-called pseudo-load vector, see Eq. (64). The accuracy of the sensitivities strongly depends on that of the pseudo-load vector. This pseudo-load vector is obtained through the derivatives of the equilibrium conditions, Eqs. (16), (17), and (20) with respect to the design variables and by assembling the individual pseudo-load vectors for each equilibrium condition. For a linear elastic analysis the pseudo-load vector can be easily obtained formulating a discrete sensitivity approach. It is more demanding for materially or geometrically nonlinear problems. Therefore the derivation of the pseudo-load vector for the damage formulation used in the present study is detailed in Sections 7.4 to 7.6. For this the gradients of constitutive equations and also the explicit part of the derivative of objective function are described first in the next two sections.

f f f ¼ rs ðTd1 ÞB~ e þ Td1 rs ðB~ Þe þTd1 B~ rs e : |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} explicit

0 1 0 1 c c c ~c @ r @ e @ r @ e c B ex im C ex c im r s rc ¼ c þ c v ¼ Cced @rs ec þ rs ec A þE @rs e~ v þ r s e~ cv A @e @s |fflffl{zfflffl} |fflfflffl{zfflfflffl} @e~ v @s

¼0

r e ¼r e

r e : e ¼e

i L

e ¼e

and

i L ðBðsÞ; dÞ;

with im c e f e f e f c rs ecL ¼ rex s eL þ r s eL ¼ r s ðT1 ÞB d þT1 r s ðB Þd þ T1 B rs d ;

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

|fflfflfflfflfflffl{zfflfflfflfflfflffl}

explicit

implicit

ð39Þ im i i rs eiL ¼ rex s eL þ r s eL ¼ r s ðBÞd þ B r s d:

ð40Þ

 local slip uiL im i i rs uiL ðNðsÞ; dÞ ¼ rex ð41Þ s uL þ r s uL ¼ rs ðNÞd þN r s d:  local and non-local equivalent strains ev , e~ v and gradient of non-local equivalent strain re~ v for

3 concrete

rs siL ¼

@siL @uiL ex im ¼ kL ðrs uiL þ rs uiL Þ; @uiL @s

ð50Þ

ð51Þ

@rc @e~ cv

E 

and

f

EL 

@sfL @e~ fv;L

:

ð52Þ

f

E and E L are detailed in Eqs. (80) and (89) in Appendices B and C. Ced is the so-called elasto-damage ‘secant’ material tensor, see Eq. (83). kL in Eq. (51) denotes the tangent modulus of the interface which is explicitly obtained from Eq. (5). In the present case the design variables s are defined on the global level and not related directly to variables on the element level needed in the present case. Therefore a semi-analytical approach is most appropriate calculating the above introduced derivatives by a finite difference method. 7.3. Sensitivity for explicit part of objective function The explicit part of sensitivity of the objective function given in Eq. (33) is expressed as follows, f i ex c rex s f ¼ r s ðf þ f þf Þ

ð53Þ

with

ecv ¼ ecv ðI1 ðec Þ; J2 ðec ÞÞ; c e vð

ex c þ s

r e ¼r e r e

ð42Þ

im c Þ¼ s

r e

c

c e vð

ex c ~v þ s

r e ¼r e

im c ~v s

r e rs ðB Þ d þ B rs dÞ; |fflfflffl{zfflfflffl}

c

c rex s f ¼ 0;

c

¼0

c s ~v

@sfL @efL @sf @e~ þ f L v;L f @eL @s @e~ v;L @s

c

f e c L ðB ðsÞ; T1 ðsÞ; dÞ

ð38Þ

c s v

ð49Þ

f

rs sfL ¼

c

ð37Þ

|fflfflffl{zfflfflffl}

c L

~c e þ E rim s ev;

with the abbreviations

im c c c c rs ec ðdÞ ¼ rex s e þ r s e ¼ r s ðB Þ d þ B r s d;

i s L

¼0

¼0

c

im c

¼ Cced rs

f

 concrete and fiber strains ecL , efL ,

c s Lþ

ð48Þ

implicit

ex im ex im ¼ Cfed;L ðrs efL þ rs efL Þ þE L ðrs e~ fv;L þ rs e~ fv;L Þ;

The gradients of the variables in the constitutive equations at position n þ 1 are calculated for the total sensitivity analyses. The derivatives with respect to a design variable s are decomposed into explicit and implicit parts:

ð47Þ

Note that the explicit terms for concrete in Eqs. (37), (43)–(45) ~ c , and B~ c do not depend on the design are zero because Nc , Bc , N variables s; in other words all explicit terms of derivatives for concrete vanish. Utilizing the above equations the stress derivatives of concrete matrix rs rc , fiber rs sfL , and interface rs siL at position n þ 1 are

7.2. Gradients of constitutive equations

f s L

f

im f ~ ~ ~f ~ rs e~ fv;L ¼ rex s e v;L þ r s e v;L ¼ rs ðN Þe þ N r s e; ex im f f f rs ðre~ v;L Þ ¼ rs ðre~ v;L Þ þ rs ðre~ v;L Þ

7. Sensitivity analysis

c

r e ¼ rs ðN~ Þ e þ N~ rs e; |fflfflfflffl{zfflfflfflffl}

f rex s f ¼

ð54Þ

Z Z Of e^ fL

ð55Þ

ð44Þ

Z Z Z Z i i i ex i i i rex ðrex siL duiL rs jJi jdOix : s f ¼ s ðsL ÞduL þ sL r s duL ÞdO  Oi u^ iL

c

c

|fflfflffl{zfflfflffl}

Ofx e^ fL

ð43Þ

¼0

im ~ ~ ~c ~c rs ðre~ cv Þ ¼ rex s ðr e v Þ þ r s ðr e v Þ ¼ r s ðB Þ e þ B rs e;

Z Z ex f f f f f ðrex sfL defL rs jJf jdOfx ; s ðsL ÞdeL þ sL r s deL ÞdO 

Oix u^ iL

ð56Þ ð45Þ

¼0

3 fibers im f f rs efv;L ðefL Þ ¼ re efv;L ðrex s eL þ r s eL Þ with Eqs: ð38Þ2ð40Þ

ð46Þ

The second terms are integrated in the parametric space x. Eq. (54) is satisfied because the functions of concrete, e.g. shape functions and B-operator, are independent of the design variables as mentioned above. The determinants of Jacobian matrices jJf j and jJi j for fiber and interface elements map the parametric element domains onto their real space. The stress

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415 ex

ex

derivatives rs sfL and rs siL are related to the explicit parts as shown in Eqs. (50) and (51). In the following sections the implicit part of sensitivity of the objective function is discussed.

Inserting Eqs. (43)–(48) into the obtained derivative of the equilibrium condition Eq. (17) and arranging it as in the previous section yields Z Z c c ~ c ÞT N ~ c ÞT F c Bc dOc rs d ~ c dOc rs e ½cðB~ ÞT B~ þ ðN ðN c c O O |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

O

Z

|fflfflfflffl{zfflfflfflffl}

Oc

¼0

Z

þ

Z

O

f

Ofx

B

fT

ðTe ÞT sf r 1

L

f

s jJ

f

jdOx rs lt þ 1

Z

O

Gx

~ Gx ¼ 0; N t0 jJjd

Z

T

Kcdd

þ

fT

fT

Kfdd

e T

Cfed;L BdOf

B ðT1 Þ rs d Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Z þ

f Kdd

Z

with the abbreviations c

F 

f

FL 

@efv;L

ð60Þ

@efL

detailed in Eqs. (81) and (90) in Appendices B and C. The virtual non-local equivalent strain de~ v in (17) is assumed to be arbitrary, thus its derivative rs de~ v vanishes.

Similarly the derivative of equilibrium condition Eq. (20) is obtained considering Eq. (29) Z Z Z T T rs ðBÞT sfL dOf þ B rs sfL dOf þ B sfL rs jJf jdOfx

Kcde

Of

Kfde

Of

rs ðNÞT siL dOi þ i

O

Ofx

Z O

T

i

N rs siL dOi þ

Z

T

i

Ox

N siL rs jJi jdOix ¼ 0:

ð61Þ Again the derivative rs duiL vanishes because duiL in (20) is an arbitrary test function. The pseudo-load vector is derived by inserting Eqs. (50) and (51) into Eq. (61). However the bond–slip relation Eq. (5) does not include any term related to the non-local equivalent strain opposite to Eq. (50). Excluding the non-local term from Eq. (50), i.e.

d P~

2 1 Z Z T f f f f fT e T f ~  Bf ðTe1 ÞT E L rex e d O  B r ðT s 1 Þ sL dO s v;L Of Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} d P~ 4

Ofx

Z þ

¼ rs lt þ 1 P Z Z T f f  rs ðBf ÞT ðTe1 ÞT sfL dOf  Bf ðTe1 ÞT Cfed;L rex s e L dO Of Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

d P~

and

7.6. Sensitivity for the third equilibrium equation

~ c dOc rs e þ B E N Oc |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}

T f f ~ dOf rs e Bf ðTe1 ÞT E L N |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

3 Z T  Bf ðTe1 ÞT sfL rs jJf jdOfx :

@ecv @ec

c

cT

Of

d P~

ð59Þ

Of

x |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} e P~ 3

Bc Cced Bc dOc rs d þ B Cfed Bf dOf rs d Oc Of |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Z

2 Z f ~ f ÞT ðe~ f ef Þrs jJf jdOf ;  ½cðB~ ÞT ðTd1 ÞT re~ fv;L þðN v;L x v;L

ð57Þ

where the virtual displacement field du in (16) is assumed to be arbitrary so that its derivative rs du vanishes. jJc j is the determinant of Jacobian matrix of the concrete element. The ~ maps a line differential on the boundary. metric jJj Substituting Eqs. (49) and (50) into Eq. (57) results in Z

Kfed

e P~

Bf ðTe1 ÞT rs ðsfL ÞdOf

cT

f

e P~

T

f

f T

1 Z ~ f ÞT ðe~ f ef Þ þðN ~ f ÞT ðrex e~ f F f rex ef ÞdOf  ½rs ðN v;L v;L s v;L L s v;L Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

¼0

½rs ðBf ÞT ðTe1 ÞT þ Bf rs ðTe1 ÞT sfL dOf þ

f

Kf

T

Z

~f T

ee Z f f f f ¼  c½rs ðB~ ÞT B~ þ ðB~ ÞT rs ðB~ ÞedOf Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

Bc rc rs jJc j dOcx |fflffl{zfflffl} Ocx

T

þ Z

T

Bc rs ðrc ÞdOc þ

Kc

Z ed f ~ ~ f ÞT F f Te Bf dOf rs d ~ ~ þ ½cðB Þ B þðN Þ N dO rs e ðN L 1 Of Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

The derivative of equilibrium condition Eq. (16) with respect to a design variable s is obtained considering Eq. (27) as follows

rs ðBc ÞT rc dOc þ c

Kcee

Z

7.4. Sensitivity for first equilibrium equation

Z

409

ð58Þ

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} d P~ 5

In Eq. (58) all implicit and explicit terms are assembled to the left and right hand side, respectively. The right hand d side of Eq. (58) leads to the pseudo-load vector P~ . Note that the stiffness matrices on the left hand side of Eq. (58) correspond to those in the tangent stiffness matrix KT introduced in (32). The same procedure is adopted for the derivatives of the second and third equilibrium conditions.

im f f rs sfL ¼ Cfed;L ðrex s eL þ r s eL Þ;

ð62Þ

and substituting Eqs. (51) and (62) into Eq. (61) results in the following expression Z  Z T T B Cfed;L BdOf þ N kL NdOi rs d Of Oi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} Ki

dd Z T B Cfed;L Te1 Bf dOf rs d ¼  rs ðBÞT sfL dOf f O Of |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

Z þ

Kf

d P~

dd

1 Z Z Z T f T f ex f f f  B Ced;L rs eL dO  B sL rs jJ jdOfx  rs ðNÞT siL dOi Of Ofx Oi |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} d P~ 2

7.5. Sensitivity for the second equilibrium equation Analogously the derivative of equilibrium condition Eq. (17) with respect to a design variable s is obtained considering Eq. (28).

d P~

3 Z Z T T i ex i i  N kL rs uL dO  N sL rs jJi jdOix : Oi Oix |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} d P~ 5

d P~ 6

d P~ 4

ð63Þ

ARTICLE IN PRESS 410

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

7.7. Total sensitivity analysis

where u j and ðrs u^ j Þpseudo are the j-th component of vectors u and ^ pseudo expressed as ðrs uÞ

Assembling Eqs. (58), (59) and (63) leads to the following compact matrix expression 3 2 5 X d ~ P l 7 32 2 2 3 6 3 6l¼1 7 cþf 7 6 Kcdeþ f Kfdd Kdd P r d 7 6 s 3 76 6 6 7 6X 7 e 7 6 7 7 6 Kc þ f Kc þ f 6 7 ~ 7: 6 0 0 r e ð64Þ  ¼ r l P 6 7 7 6 ed 6 7 s s t þ 1 ee l 7 54 4 4 5 6 5 6l¼1 7 0 rs d 7 6 Kfdd 0 Kidd 7 6X 4 6 ~d 5 Pl

^ u ¼ K1 T lt þ 1 P;

ð67Þ

~ ^ pseudo ¼ K1 ðrs uÞ T P pseudo :

ð68Þ

l¼1

Substituting u j and ðrs u^ j Þpseudo into (66) yields

rs lt þ 1 ¼ 

ðrs u^ j Þpseudo lt þ 1 : uj

ð69Þ

According to the above equations the derivative of the total nodal displacement vector rs u^ is

r s lt þ 1

which has the format of the typical stiffness equation adding up all terms on the right hand side to a new pseudo-load vector

rs u^ ¼ u

Ppseudo :

Finally, the total sensitivity of the objective function can be obtained by inserting Eq. (70) into Eq. (36) and accumulating each sensitivity over the load increment step number nstep as

KT rs u^ ¼ Ppseudo ¼ rs lt þ 1 P^ þ P~ pseudo :

ð65Þ

KT denotes the tangent stiffness matrix at the time step t þ1 as mentioned before. The next question is how to deal with the derivative of load factor rs l. Note that the derivatives based on a load-controlled algorithm differ from those based on a displacement-controlled algorithm controlling a certain nodal displacement ‘component’ uj ¼ u^ j of the structure which is more suitable for the optimization of ductility. For a detailed description it is referred to Schwarz et al. [33], Lipka et al. [21]. A load-controlled algorithm renders rs l ¼ 0 while for a displacement controlled algorithm only the sensitivity of the nodal displacement for the controlled degree of freedom u^ j is equal to zero. The sensitivity of the load factor based on a discretized formulation is derived subsequently, u rs u^ j ¼ rs lt þ 1 j þ ðrs u^ j Þpseudo ¼ 0; ð66Þ

lt þ 1

rs f ¼

lt þ 1

n step X

^ pseudo : þ ðrs uÞ

ex

ðrs f t þ ru^ f t rs u^ t Þ

ð70Þ

ð71Þ

t¼1

where f t indicates the ductility increment in the t-th load increment. Note that in the damage model the stress r is an explicit function of the total strain e even if un-/re-loading situations occur. Thus an incremental iterative procedure for the update of the stress r, mandatory for example in plasticity, is not necessary. This allows to use directly the derivative of the total strain rs e or displacement rs u at position n þ 1 in the sensitivity analysis, i.e. both derivatives do not have to be accumulated over load increments. This also avoids that errors of the sensitivities are accumulated during load incrementation.

ν κ

Fig 7. Deep beam.

ν κ

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

411

0.250

0.125

0 0

0.005

12 8 4 0 0

0.2

0.4

Fig 8. Results of optimization for (a) a linear elastic and (b) a materially nonlinear response.

As said above this approach is restricted to material models in which the stress r is a function of the total strain e. Details of this direct sensitivity approach using a total strain or displacement is described in Bugeda et al. [7].

8. Numerical examples 8.1. Optimization for a deep beam In the first numerical simulation a FRC beam reinforced with four carbon fibers is investigated as displayed in Fig. 7 where also the material properties are given. For the properties of the ¨ interface it is referred to Kruger et al. [18–20], Xu et al. [37]. Plane-stress conditions are assumed. Due to symmetry only one half of the system is analyzed, the FE mesh is given in Fig. 7(c). The beam thickness is assumed to be only 1 mm, since no out-ofplane actions are considered. The 200 elements are used for concrete and 68 elements for the interface, respectively. The number and location of slip nodes may change during optimization depending on the actual fiber geometries. Thus, mesh adaptation for slip nodes is carried out after each structural analysis. The fiber geometry is approximated by a symmetric biquadratic Be´zier-spline, see Fig. 7(b). As shown in Fig. 7(a) the origin O of the parametric element is defined at the edge of the beam. Due to symmetry the number of design variables is reduced, i.e. the location of the control points p3 and p4 is coupled to p1 and p0 , respectively. Further simplification of the fiber geometry can reduce the number of design variables. Firstly, the y-coordinate of p1 is set equal to that of p2 . Secondly, the x-coordinate of p1 is placed always at the center between the x-coordinates of p0 and p2 . Thus the number of design variables for a single fiber is three, i.e. s1 , s2 , and s3 , see Fig. 7(b). The initial set of the design variables is (i) s1 ¼ 0:075 (i.e. the x-coordinate of p0 is 0:075  400 mm) for all fibers and (ii) s2 and s3 are assumed to be 0.15, 0.38, 0.62, and 0.85 for each fiber.

Fig 9. Hanging deep beam.

The total number of design variables is 12 (3  4 fibers). Taking into account that thick concrete covers for textile fibers are obsolete, we adopt the lower bound sL ¼ 0:01 and the upper one sU ¼ 0:99 for s2 and s3 of all fibers. For the design variable s1 , the lower and upper bounds are set to sL ¼ 0:01 and sU ¼ 0:4, respectively. The analysis is carried out with a displacement controlled method; the control point c is at the lower center of the beam. For comparison the structure is optimized based on either a linear elastic or the damage model. The prescribed nodal displacement u^ (y-direction) at the control point is either 0.005 or 0.4 mm. The fiber volume is kept constant (1.4%) during the optimization leading to a fiber thickness tf ¼ 0:4 mm. Firstly we optimize the ductility of the structure for a linear elastic response, which means maximizing the overall stiffness of the structure. Fig. 8(a) shows the optimized fiber layout. The figure on the right side of Fig. 8(a) introduces the stress distribution of fibers. After optimization the two fibers are shifted to the upper part in compression and the two others to the lower

ARTICLE IN PRESS 412

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

12

2

8 1 4 0

0 0

0.005

0

0.02

0.04

0.06

Fig 10. Results of optimization for (a) a linear elastic and (b) a materially nonlinear response.

Fig 11. Problem description of the third example and parametric element.

25 20 15 10 5 0 0

0.1

Fig 12. Results of optimization for (a) a linear elastic and (b) a materially nonlinear response.

0.2

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

one carrying the tension force, respectively. The two upper fibers wind up with almost the same location. As a result an increase of 14% of ductility could be obtained. Fig. 8(b) shows the optimized fiber layout in the materially nonlinear situation applying the damage model. Also the damage distribution of concrete is displayed; fibers are not yet damaged at this stage. After optimization one fiber is shifted to the upper part and the three others to the lower one. These three fibers prevent the concrete from a premature damage propagation. Compared to the elastic case the fibers are more curved which is structurally reasonable. As a result an increase of 44% of ductility could be obtained. One could expect that also the fourth fiber moves to the lower part. This suggests that the achieved solution does not represent the global minimum, a consequence of the underlying nonconvex optimization problem. 8.2. Optimization for a hanging deep beam As second numerical example a hanging deep beam is chosen as displayed in Fig. 9. The material properties of concrete, carbon fibers and interface as well as loading condition and mesh follow the previous example. Due to symmetry only one half of the system is analyzed under plane-stress conditions. Again a symmetric biquadratic Be´zier-spline is adopted to approximate each fiber geometry. 180 elements are used for concrete and 50 elements for the interface. The initial set of design variables for all three fibers in the parametric space is (i) s1 ¼ 0:025 and (ii) s2 and s3 equal 0.25, 0.5 or 0.75. The fiber volume is kept constant (1.1%) during the optimization. Fig. 10(a) and (b) show the optimized fiber layouts based on a linear elastic and a materially nonlinear response, respectively. The prescribed displacement u^ at the control point c under the load is either 0.005 or 0.06 mm. For the linear elastic case the upper straight fiber reduces the compressive deformation and the middle curved fiber reflects the cable effect between the fixed supports. For the damage case, the upper and middle fibers are utilized to resist the damage propagation of concrete in the vicinity of the supports. As a result an increase of 13% of ductility was obtained. 8.3. Optimization for a splitting plate The third example is a splitting plate shown in Fig. 11. Again the same material properties, loading condition, mesh and initial assumption of fiber geometry are used as in the previous two examples. The 124 elements are used for concrete and 24 elements for the interface. For the present example the parametric element is restricted to the area below the cutout section, see Fig. 11; this means that fibers cannot be located in the non-design space. The initial set of the design variables is: (i) s1 ¼ 0:025 and (ii) s2 and s3 are 0.25, 0.50 or 0.75. The fiber volume is kept constant (0.74%) during the optimization. Fig. 12(a) and (b) show the optimized fiber layouts based on a linear elastic and a damage response, respectively. The prescribed displacement at control point c is either 0.002 or 0.2 mm. For the linear elastic case, the upper straight fiber controls the tensile deformation around the reentrant corners. The middle and lower fibers reduce the compressive deformation. As a result an increase of 5% of ductility was obtained. In case of damage the middle fiber is also shifted to the upper part to resist the damage propagation of concrete together with the upper fiber. These two fibers are damaged at the

413

prescribed displacement. The location of the lower fiber stays in the lower part of the plate although we could expect that it also moves as well to the cutout area (see comment on local minimum for first example). Anyhow an increase of 99% of ductility could be obtained.

9. Conclusions Shape optimization is applied to fiber geometry for textile fiber reinforced composites. The two brittle materials, namely concrete matrix and fibers, get the necessary ductility from the interface behavior of the two constituents. It could be shown that the main purpose of the present study, namely to increase the structural ductility of FRC with respect to the geometrical layout of continuous fibers, was successful. For this objective, it is of course not sufficient to base the optimization process on a linear material model; thus it is mandatory to consider material nonlinearities in the optimization process. As an example, we applied a damage model to both constituents together with a nonlinear interface model. The formulation could be extended to other material models and other composites. It can also be enhanced by material optimization varying the fiber content in a design element as it was introduced in Kato et al. [16]. For the sensitivity analysis a variational semi-analytical direct method was applied. The strategy of sensitivity analysis using a ‘total’ strain or displacement was discussed. This approach does not accumulate errors as it is the case when an ‘incremental’ sensitivity approach is applied, which is for example necessary for plasticity (e.g. Lipka et al. [21], Maute et al. [23]). In semi-analytical methods errors of sensitivities tend to increase when distinct ‘rigid body rotations’ appear, see for example Olhoff and Rasmussen [26], Olhoff et al. [27], Cheng and Olhoff [9], Mlejnek [25], Boer and Keulen [4,5], Keulen and Boer [17], and Bletzinger et al. [3]. This tendency was also observed in the present study when the structural response reaches the postpeak range, where several elements are severely damaged and other elements may already be in the unloading phase, similar to a ‘plastic hinge deformation’. The unloaded elements eventually will encounter rigid body rotations which in turn may lead to inaccurate sensitivities. However in the present study the structural damage is not driven so far into complete failure.

Acknowledgments The present study is supported by Grants of the ‘‘Deutsche Forschungsgemeinschaft’’ DFG (German Research Foundation) within the Research Projects Ra 218/19 and Ra 218/21. This support is gratefully acknowledged.

Appendix A. Transformation matrices The transformation matrices introduced in the text are listed here. First the usual rotation matrix Td is given " # " # cosðxG ; xL Þ cosðyG ; xL Þ l1 m1 ¼ ; ð72Þ Td ¼ cosðxG ; yL Þ cosðyG ; yL Þ l2 m2 where lk and mk ðk ¼ 1; 2Þ are introduced as abbreviations. The strain transformation matrix Te transforms the ‘vector’ eG with the components of the global strain tensor to the local one eL

ARTICLE IN PRESS 414

J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415 c

so that 2

e11

3

7 e eL ¼ 6 4 e22 5 ¼ T eG : 2e12

ð73Þ

Incidentally, the stress transformation matrix Ts which transforms the global stress ‘vector’ rG to the local one rL follows as 2 3

s11 6s 7 rL ¼ 4 22 5 ¼ Ts rG ; s12 where 2

un/reloading. Thus E is a term which controls damage and the loading condition simultaneously. The same situation also holds for the damage formulation of fibers. In the above linearization procedure the following two relations are derived from Eqs. (80), (22) and from Eq. (81) @rc @rc @e~ c c c ~ ; ¼ c v ¼E N @e @e~ v @e

@ecv ðeÞ @ec @ec c ¼ F Bc ; ¼ vc @e @d @d

ð84Þ

ð74Þ Appendix C. Linearization of gradient enhanced damage model for fiber

3

2

m21 l1 m1 l 6 12 7 6 7 T ¼ 4 l2 m22 l2 m2 5 2l1 l2 2m1 m2 l1 m2 þl2 m1 3 2 2 m21 2l1 m1 l1 7 6 2 7 Ts ¼ 6 m22 2l2 m2 5 4 l2 l1 l2 m1 m2 l1 m2 þ l2 m1 e

and

The derivatives of Eqs. (27) and (28) with respect to nodal displacements d and nodal non-local strains e lead to four stiffness matrices for the gradient enhanced damage model for fiber ð75Þ nele Z f [ @f int;u T ¼ Bf ðTe1 ÞT Cfed;L Te1 Bf dOf f @d e¼1 O f

The stress transformation matrix Ts can be replaced by the relationship ðTs Þ1 ¼ ðTe ÞT .

Kfdd ¼

[ Z

nfele

¼

O

e¼1

Appendix B. Linearization of gradient enhanced damage model for concrete

T

f

Bf Cfed;G Bf dOf ;

ð85Þ

Kfde ¼

nele Z f [ @f int;u T f f ~ dOf ; ¼ Bf ðTe1 ÞT E L N f @e e¼1 O

ð86Þ

Kfed ¼

nele Z f [ @f int;e ~ f ÞT F f Te Bf dOf ; ¼ ðN L 1 f @d e¼1 O

ð87Þ

Kfee ¼

nele Z f [ @f int;e f f ~ f ÞT N ~ f dOf ¼ ½cðB~ ÞT ðTd1 ÞT Td1 B~ þ ðN f @e e¼1 O

f

The derivatives of Eqs. (27) and (28) with respect to nodal displacements d and nodal non-local strains e lead to the following four stiffness matrices for the gradient enhanced damage model for concrete nele Z c [ @f int;u T ¼ Bc Cced Bc dOc ; Kcdd ¼ c @d O e¼1 c

nele Z c [ @f int;u T c ~c ¼ ¼ Bc E N dOc ; c @e O e¼1

ð76Þ f

c

Kcde

nele Z c [ @f int;e ~ c ÞT F c Bc dOc ; ¼ ðN c @d e¼1 O

f

ð77Þ

c

Kced ¼

nele Z c [ @f int;e c c ~ c ÞT N ~ c dOc ; ¼ ¼ ½cðB~ ÞT B~ þ ðN c @e e¼1 O

ð78Þ

ð79Þ

with

f

@rc @rc @Dc @kc E  c ¼ ; @Dc @kc @e~ cv @e~ v c

F 

EL 

FL 

c

@ecv @ec @I1 @ec @J ¼ v c þ v 2c ; c @e @I1 @e @J2 @e

f

ð88Þ

@sfL

¼

@ ~ fv;L

@sfL @Df @kf ; @Df @kf @e~ fv;L

e

ð89Þ

@efv;L @efL

:

ð90Þ

@sfL @sf @e~ f ¼ fL L ¼ Cfed;L Te1 Bf ; @d @e~ L @d

@J2 1 @ecaa ecbb 1 @ecab ecab 1 ¼  ¼ ecaa dkl eckl ; c @e 6 @eckl 2 @eckl 3

where d is the Kronecker symbol. is the matrix of the ‘secant’ material tensor for isotropic elasto-damage with rc ¼ ð1Dc ÞCcel ec ;

f

ð81Þ

ð82Þ

@rc ¼ ð1Dc ÞCcel @ec

f

In the linearization procedure the following two relations are derived from Eq. (26) and from Eqs. (89) and (22),

Cced

Cced ¼

f

~ ÞT N ~ dOf ; ½cðB~ ÞT B~ þðN

ð80Þ

and @ec @I1 ¼ cii ¼ dik dil ¼ dkl ; c @e @ekl

f e¼1 O

with f

c

Kcee

[ Z

nfele

¼

@sfL @sf @e fv;L f ¼ fL ¼ ELNf ; @e @e v;L @e

ð91Þ

The two remaining relations are derived from Eqs. (90), (26) and from Eq. (31), respectively, @efv;L ðefL Þ @ef @efL f ¼ F L Te1 Bf ; ¼ v;L @d @efL @d

@ðcre~ fv;L Þ @tfL f ¼ cTd1 B~ : ¼ @e @e

ð92Þ

ð83Þ

where Ccel denotes the matrix of the elastic material tensor, and D is the damage parameter. The second term on the right hand side of Eq. (80) is zero unless damage is initiated in an element. The third term is equal to either unity for loading or zero for

Appendix D. Linearization of interface element The derivatives of Eq. (29) with respect to nodal displacements d and nodal slip parameters d lead to the stiffness matrices for

ARTICLE IN PRESS J. Kato, E. Ramm / Finite Elements in Analysis and Design 46 (2010) 401–415

the interface Kfdd ¼

f @f int;u

@d

[ Z

nfele

¼

e¼1

T

O

f

Bf ðTe1 ÞT Cfed;L BdOf ;

ð93Þ

nele Z i [ @f int;i T ¼ B Cfed;L Te1 Bf dOf ; f @d e¼1 O f

Kfdd ¼

@f int;i @d

[ Z

niele

i

Kidd ¼

¼

e¼1

T

O

f

B Cfed;L BdOf þ

ð94Þ

Z O

i

 T N kL NdOi ;

ð95Þ

where the following two relations are derived from Eq. (26) and from Eq. (23), respectively, @sfL @d

¼

@sfL @efL ¼ Cfed;L B; @efL @d

@siL @d

¼

@siL @uiL ¼ kL N: @uiL @d

ð96Þ

Stiffness matrix Kfdd has already been introduced in Eq. (85). References [1] S. Balakrishnan, D.W. Murray, Finite element prediction of reinforced concrete behavior, Structural Engineering Report no. 138, University of Alberta, Edmonton, Alberta, Canada, 1986. [2] F. Barzegar, S. Maddipudi, Generating reinforcement in FE modeling of concrete structures, J. Struct. Eng. 120 (5) (1994) 1656–1662. [3] K.-U. Bletzinger, M. Firl, F. Daoud, Approximation of derivatives in semianalytical structural optimization, Comput. Struct. 86 (2008) 1404–1416. [4] H. de Boer, F. van Keulen, Error analysis of refined semianalytical design sensitivities, Struct. Optim. 14 (1997) 242–247. [5] H. de Boer, F. van Keulen, Refined semi-analytical design sensitivities, Int. J. Solids Struct. 37 (2000) 6961–6980. [6] W. Brameshuber, T. Brockmann, J. Hegger, M. Molter, Untersuchungen zum textilbewehrten Beton, Beton 52 (2002) 424–429. ˜ate, Structural shape sensitivity analysis for nonlinear [7] G. Bugeda, L. Gil, E. On material models with strain softening, Struct. Optim. 17 (1999) 162–171. [8] T.Y. Chang, H. Taniguchi, W.F. Chen, Nonlinear finite element analysis of reinforced concrete panels, J. Struct. Eng. ASCE 113 (1) (1987) 122–140. [9] G. Cheng, N. Olhoff, Rigid body motion test against error in semi-analytical sensitivity analysis, Comput. Struct. 46 (3) (1993) 515–527. [10] A.E. Elwi, T.M. Hrudey, Finite element model for curved embedded reinforcement, J. Eng. Mech. 115 (4) (1989) 740–754. [11] G. Hofstetter, H.A. Mang, Work-equivalent node forces from prestress of concrete shells, in: T.J.R. Hughes, E. Hinton (Eds.), Finite Element Methods for Plate and Shell Structures, Formulations and Algorithms, vol. 2, Pineridge Press, Swansea, 1986, pp. 312–347. [12] F. Huber, Nichtlineare dreidimensionale Modellierung von Beton- und ¨ Baustatik, Universitat ¨ Stahlbetontragwerken, Ph.D. Thesis, Institut fur Stuttgart, Germany, 2006 (/http://www.ibb.uni-stuttgart.deS). [13] A. Hund, Mehrskalenmodellierung des Versagens von Werkstoffen mit ¨ Baustatik, Universitat ¨ Stuttgart, Mikrostruktur, Ph.D. Thesis, Institut fur Germany, 2007 (/http://www.ibb.uni-stuttgart.deS). [14] B.L. Karihaloo, D. Lange-Kornbak, Optimization techniques for the design of highperformance fiber-reinforced concrete, Struct. Multidisc. Optim. 21 (2001) 32–39. [15] J. Kato, A. Lipka, E. Ramm, Preliminary investigation for optimization of fiber reinforced cementitious composite structures. in: C.A. Mota Soares, et al. (Eds.), Proceedings of Third European Conference on Computational Mechanics, Solids, Structures and Coupled Problems in Engineering, ECCM, Lisbon, Portugal, 2006.

415

[16] J. Kato, A. Lipka, E. Ramm, Multiphase material optimization for fiber reinforced composites with strain softening, Struct. Multidisc. Optim. online published (2008) doi:10.1007/s00158-008-0315-7. [17] F. van Keulen, H. de Boer, Rigorous improvement of semi-analytical design sensitivities by exact differentiation of rigid body motions, Int. J. Numer. Methods Eng. 42 (1998) 71–91. ¨ [18] M. Kruger, J. Ozbolt, H.-W. Reinhardt, A discrete bond model for 3D analysis of textile reinforced and prestressed concrete elements, Otto-Graf-Journal 13 (2002) 111–128. ¨ [19] M. Kruger, J. Ozbolt, H.-W. Reinhardt, A new 3D discrete bond model to study the influence of bond on the structural performance of thin reinforced and prestressed concrete plates. in: H.-W. Reinhardt, et al. (Eds.), Proceedings of High Performance Fiber Reinforced Cement Composites (HPFRCC4), RILEM, Ann Arbor, USA, 2003, pp. 49–63. ¨ [20] M. Kruger, S. Xu, H.-W. Reinhardt, J. Ozbolt, Experimental and numerical studies on bond properties between high performance fine grain concrete ¨ and carbon textile using pull out tests, in: Beitrage aus der Befestigungstechnik und dem Stahlbetonbau, Festschrift Professor R. Eligehausen ¨ Stuttgart, Germany, 2002, pp. 151–164. Universitat [21] A. Lipka, S. Schwarz, E. Ramm, Topology optimization of three-dimensional structures with consideration of elastoplastic structural response, in: Z. Waszczyszyn, et al. (Eds.), Proceedings of Second European Conference on Computational Mechanics, Solids, Structures and Coupled problems in Engineering, ECCM, Cracow, Poland, 2001. ¨ [22] K. Maute, Topologie- und Formoptimierung von dunnwandigen Tragwerken ¨ Baustatik, Universitat ¨ Stuttgart, Germany, 1998. Ph.D. Thesis, Institut fur [23] K. Maute, S. Schwarz, E. Ramm, Adaptive topology optimization of elastoplastic structures, Struct. Optim. 15 (1998) 81–91. [24] J. Mazars, G. Pijaudier-Cabot, Continuum damage theory—application to concrete, J. Eng. Mech. 115 (1989) 345–365. [25] H.P. Mlejnek, Accuracy of semi-analytical sensitivities and its improvement by the natural method, Struct. Optim. 4 (1992) 128–131. [26] N. Olhoff, J. Rasmussen, Study of inaccuracy in semi-analytical sensitivity analysis—a model problem, Struct. Optim. 3 (1991) 203–213. [27] N. Olhoff, J. Rasmussen, E. Lund, A method of exact numerical differentiation for error elimination in finite-element-based semi-analytical shape sensitivity analysis, Mech. Struct. Mach. 21 (1) (1993) 1–66. [28] S.N. Patnaik, J.D. Guptill, L. Berke, Merits and limitations of optimality criteria method for structural optimization, Int. J. Numer. Methods Eng. 38 (1995) 3087–3120. [29] R.H.J. Peerlings, R. de Borst, W.A.M. Brekelmans, J.H.P. de Vree, Gradient enhanced damage for quasi-brittle materials, Int. J. Numer. Methods Eng. 39 (1996) 3391–3403. [30] R.H.J. Peerlings, R. de Borst, W.A.M. Brekelmans, M.G.D. Geers, Gradientenhanced damage modelling of concrete fracture, Mech. Cohes. Frict. Mater. 3 (1998) 323–342. [31] R.H.J. Peerlings, Enhanced damage modelling for fracture and fatigue Ph.D. Thesis, Technische Universiteit Eindhoven, The Netherlands 1999. [32] D.V. Phillips, O.C. Zienkiewicz, Finite element nonlinear analysis of concrete structures, Proc. Inst. Civ. Eng. Part 2 61 (3) (1976) 59–88. [33] S. Schwarz, K. Maute, E. Ramm, Topology and shape optimization for elastoplastic structural response, Comput. Methods Appl. Mech. Eng. 190 (2001) 2135–2155. [34] J. Stegmann, E. Lund, Discrete material optimization of general composite shell structures, Int. J. Numer. Methods Eng. 62 (2005) 2009–2027. [35] M. Stolpe, J. Stegmann, A Newton method for solving continuous multiple material minimum compliance problems, Struct. Multidisc. Optim. 35 (2008) 93–106. [36] J.H.P. de Vree, W.A.M. Brekelmans, M.A.J. Gils, Comparison of nonlocal approaches in continuum damage mechanics, Comput. Struct. 55 (1995) 581–588. ¨ [37] S. Xu, M. Kruger, H.-W. Reinhardt, J. Ozbolt, Bond characteristics of carbon, alkali resistance glass, and aramid textiles in mortar, J. Mater. Civil Eng. ASCE 16 (4) (2004) 356–364.