Isogeometric Shape Optimization for Quasi-static and Transient Problems

Isogeometric Shape Optimization for

Quasi-static and Transient Problems

865229

789461

9

Zhenpei Wang

Zhenpei W

ang

Isogeometric Shape Optimization for Quasi-static and T

I

SOGEOMETRIC

S

HAPE

O

PTIMIZATION FOR

Q

UASI

-

STATIC AND

T

RANSIENT

P

ROBLEMS

I

SOGEOMETRIC

S

HAPE

O

PTIMIZATION FOR

Q

UASI

-

STATIC AND

T

RANSIENT

P

ROBLEMS

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K. C. A. M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op woensdag 20 januari 2015 om 15.00 uur

door

Zhenpei WANG

Master of Science, Northwestern Polytechnical University, Xi’an, China geboren te Baofeng, Henan Province, China

promotor: Prof. dr. C. Bisagni copromotor: Dr. S. R. Turteltaub Composition of the doctoral committee: Rector Magnificus

Prof. dr. C. Bisagni Dr. S. R. Turteltaub

voorzitter

Technische Universiteit Delft, promotor Technische Universiteit Delft, copromotor

Independent members: Prof. dr. -ing. K. -U. Bletzinger Prof. dr. P. Duysinx

Prof. dr. ir. F. van Keulen Dr. M. M. Abdalla Dr. ir. H. J. M. Geijselaers Prof. dr. ir. R. Benedictus

Technische Universität München, Germany Université de Liège, Belgium

Technische Universiteit Delft Technische Universiteit Delft Universiteit Twente

Technische Universiteit Delft, Reserve member

Keywords: Isogeometric analysis, Shape optimization, Continuous adjoint

method, Quasi-static, Transient heat conduction, Normalization approaches, Discretization/mesh-dependency

ISBN 978-94-6186-522-9

Email: wangzhenpei@gmail.com; z.p.wang@foxmail.com An electronic version of this dissertation is available at

http://repository.tudelft.nl/.

C

ONTENTS

Nomenclature xi

1 Introduction 1

1.1 Background. . . 1

1.1.1 Isogeometric shape optimization . . . 1

1.1.2 Structural optimization for time-dependent problems. . . 2

1.1.3 Adjoint shape sensitivity analysis . . . 3

1.1.4 Shape updating scheme . . . 3

1.2 Scope and aims. . . 4

1.3 Layout and notation scheme . . . 5

2 Isogeometric analysis 7 2.1 NURBS . . . 7

2.1.1 B-spline . . . 7

2.1.2 NURBS. . . 8

2.1.3 Basic properties of the basis functions. . . 11

2.1.4 Mean value property of a B-spline basis function over its local support . . . 11

2.1.5 Geometry operations . . . 13

2.2 NURBS-based Isogeometric analysis . . . 14

2.3 Isogeometric analysis for quasi-static mechanical problems . . . 15

2.3.1 Quasi-static problems . . . 15

2.3.2 Weak formulation . . . 16

2.3.3 Isogeometric discretization . . . 16

2.4 isogeometric analysis for transient heat conduction problems . . . 18

2.4.1 Transient heat conduction problems. . . 18

2.4.2 Weak formulation . . . 19

2.4.3 Isogeometric discretization . . . 19

3 Transport relations considering discontinuities 21 3.1 Continuous design parametrization . . . 21

3.2 Material and spatial design derivatives . . . 22

3.3 Parametrization of functions with discontinuities. . . 24

3.4 Volume transport relations with discontinuities. . . 25

3.4.1 Volume divergence theorem. . . 25

3.4.2 Derivative of a volume integral. . . 26

3.4.3 Volume transport relation considering discontinuities. . . 26 vii

3.5 Surface transport relations with discontinuities. . . 28

3.5.1 Surface divergence theorem . . . 28

3.5.2 Derivative of a surface integral. . . 28

3.5.3 Surface Transport relation considering discontinuities. . . 29

4 Normalization approaches for the descent search direction in isogeomet-ric shape optimization 31 4.1 Continuous and discrete gradients and search directions. . . 31

4.1.1 Continuous shape gradient and search direction. . . 31

4.1.2 Discrete shape gradient and search direction . . . 32

4.2 Discretization-dependency. . . 33

4.3 Normalization of the search directions . . . 37

4.3.1 Consistent discretization. . . 37

4.3.2 Standard normalization . . . 37

4.3.3 A simple example on the quadratic norm induced by the dis-cretization. . . 38

4.3.4 Diagonally-lumped mapping matrix (DLMM) normalization . . 39

4.3.5 Simplified DLMM normalization. . . 40

4.4 Verification. . . 42

4.4.1 Volume minimization of an initially rectangular domain. . . 42

4.4.2 Volume minimization of an initially elliptic plate . . . 42

4.5 Performance of normalization approach . . . 45

4.5.1 3D fillet shape design optimization . . . 45

4.5.2 Heat conduction problem . . . 48

4.6 Discussions and remarks . . . 51

4.6.1 Discussions on the cases that normalization is unnecessary. . . 51

4.6.2 Remarks. . . 52

5 Shape design optimization for quasi-static problems 55 5.1 Design problem statement . . . 55

5.2 Continuous adjoint shape sensitivity analysis. . . 58

5.2.1 Augmented functional. . . 59

5.2.2 Shape derivatives of Ψ. . . 60

5.2.3 Shape derivatives of the equilibrium constraints. . . 61

5.2.4 Adjoint model . . . 62

5.2.5 Local shape derivative. . . 65

5.2.6 Gradient of constrained, quasi-static objective functional. . . . 66

5.3 Discrete isogeometric shape sensitivity. . . 67

5.3.1 Sensitivity analysis based on the global shape derivative. . . 67

5.3.2 Sensitivity analysis based on the local gradient . . . 70

5.4 Iterative descent method . . . 70

5.5 Benchmark examples. . . 70

5.5.1 Benchmark static problems . . . 72

5.5.2 Verification of the time-dependent problems . . . 77

5.6 Application examples. . . 79

CONTENTS ix

6 Shape design optimization for transient problems 89

6.1 Design problem statement . . . 89

6.2 Continuous adjoint shape sensitivity analysis. . . 91

6.2.1 Augmented functional. . . 91

6.2.2 Shape derivatives of ˜J. . . 92

6.2.3 Transient adjoint system. . . 94

6.2.4 Continuous adjoint sensitivity . . . 95

6.3 Discrete isogeometric shape sensitivity. . . 97

6.3.1 Isogeometric analysis and design discretization space. . . 97

6.3.2 Discrete shape gradients with respect to design control points. . 97

6.3.3 Normalization of the search direction (sensitivity weighting) . . 98

6.3.4 Iterative descent method. . . 98

6.4 Numerical verifications and applications. . . 98

6.4.1 Verification example. . . 98

6.4.2 Shape optimization of a plunger. . . 102

6.4.3 Shape optimization of a thermal protection panel. . . 106

6.5 Remarks. . . 108

7 Conclusions and Recommendations 111 7.1 Conclusions and highlights. . . 111

7.1.1 Mean value property of B-spline basis. . . 111

7.1.2 Continuous shape sensitivity analysis considering discontinu-ities . . . 111

7.1.3 Normalization approaches for descent search direction in Iso-geometric shape optimization . . . 112

7.1.4 Isogeometric shape optimization for quasi-static problems . . . 112

7.1.5 Isogeometric shape optimization for transient problems . . . . 113

7.2 Future work. . . 114

References 115

Summary 123

Curriculum Vitae 127

Acknowledgements 129

Papers and conference presentations 131

List of Figures 133

N

OMENCL ATURE

G

REEK

S

YMBOLS

α = step-size

β = numerical approach parameter for transient heat conduction problems

γ = boundary characteristic function

Γ = boundary

Γu, Γt = displacement and traction boundary

Γq, Γθ = types of boundary with specified heat flux and temperature, respectively

Γe = type of boundary with heat convection conditions

δ = variational operator

ǫ = strain field

ζ = NURBS parameter

η = NURBS parameter

θ, θ = temperature field

ϑ = adjoint temperature field

κ = mean curvature λ = Lagrangian multiplier Λ = Lagrangian multiplier ν = design velocity ξ = NURBS parameter ρ = material density σ = stress field Σ = volume function

ς = time characteristic function

τ = adjoint time-like parameter

Υ = boundary/border of the sub-boundaries

φ = objective functional

χ = NURBS parameter

ψ = objective functional

ω = domain/region characteristic function

Ω = domain/region

∇ = gradient operator

L

ATIN

S

YMBOLS

A = family of feasible design functions

B = strain-displacement matrix

C = global capacitance matrix

c = heat conductivity coefficient C _{= stiffness tensor}

dc = design constraint

D, ¯D = integration domain

E = material property matrix

F = load vector

F = objective density function

f = body force

f = a generic function defined in a domain field

g = local shape gradient (continuous)

gI, G = discrete shape gradient

h = a generic function defined on the boundary

h = a generic spatial vector field defined on the boundary ¯

h = a generic tangential vector field defined on the boundary J = objective functional

K = stiffness matrix

l = a functional that represents the governing equation

L = augmented Lagrangian functional

Mp+1 = fundamental spline function with a degree of p M = mass-matrix-liked matrix

m = a unit vector that is tangential to the surface and normal to its boundary

m = number of control points

N = B-spline basis function

n = outward normal vector

n = number of control points

p = degree of the NURBS

p = material point

q = contact heat flux

q = heat flux vector

r = degree of the NURBS

R, R = NURBS basis function

s = time-like design parameter for the continuous formulation

Sω = internal boundary of the sub-domain Ωω

Sa = total number of control points used for analysis

Sb = total number of control points for the design boundary

St = total number of control points

t , T = time

T = time interval [0, T ]

ACRONYMS xiii

u = displacement field

U = discrete displacement vector

w = NURBS weigh functions

x = location vector

S

UBSCRIPTS

a = analysis

b = boundary

c = capacitance; curve; continuous; constraint; control points

d = design

e = exchange

f = body force

i , j , k = dimensions; control points index

l , m = numbers of control points

n = normalized; number of control points

p, r = degrees of the spline basis

q = heat flux; degree of spline basis

s = surface

t = traction

v = volume

ω = characteristic function for sub-domain/volume

γ = characteristic function for sub-boundary

ς = characteristic function for time sub-interval

Σ = volume

S

UPERSCRIPTS

I = mapping functions of the control points index

i , j , k = dimensions; control points index

l , m, n = numbers of control points

s = time-like design parameter

T = transpose

ˆ

(∗) = prescribed fields ˇ

(∗) = artificially chosen target fields ˚

(∗) = material design derivatives (∗)′ = spatial design derivatives

˙

(∗) = time derivatives (∗)∗ = adjoint fields

A

CRONYMS

B-spline = basis spline

NURBS = Non-uniform rational B-spline IGA = isogeometric analysis

CAD = computer aided design CAE = computer aided engineering FEM = finite element method FE = finite element

DLMM = diagonally lumped mapping matrix SD = search direction

1

I

NTRODUCTION

1.1.

B

ACKGROUND

1.1.1.

I

Structures and structural components are typically designed in a computer-aided design (CAD) environment using compact and efficient representations of the ge-ometry. In currently-used CAD software, the geometrical representation relies on non-uniform rational B-splines (NURBS) [PT97]. The concept of B-splines have been studied for decades since it was introduced by Schoenberg [Sch46, CS47]. Early fundamental work [CS66,Cox72,DB72,Rie73] related to the properties of B-spline basis functions continued through the early 1970’s. Following that, Non-Uniform Rational B-Splines (NURBS) were introduced to provide better flexibility in geometric modeling by Versprille in his PhD thesis [Ver75]. Detailed explanations of the definitions and properties of B-splines and NURBS are widely available; e.g., [BFK84,PT97,Sch07].

Currently-used finite-element software packages are often based on Lagrange polynomials for analysis. Consequently, to carry out the analysis of a CAD-generated structure for shape optimization, it is required to convert its NURBS-based represen-tation into one compatible with Lagrange polynomials (i.e., generation of a mesh, see, e.g., [BF84], [HG86] and [WZ12]). The conversion from NURBS to Lagrange polynomials may be computationally costly and may cause a loss of valuable infor-mation. This problem is aggravated in the context of shape optimization since, for the shape sensitivity analysis, it is important to measure accurately the effect that a small change in geometry has on the performance of a structure. Although it is possible to carry out the sensitivity analysis directly on a finite-element mesh (see, e.g., [LBT11]), the resulting design is expressed in terms of a mesh that is not compatible with a CAD representation.

The properties of the NURBS basis functions make them possible to be used as the shape functions for analysis. Aimed at integrating computer aided design (CAD) and analysis, a NURBS-based finite element method, namely isogeometric analysis

1

(IGA) was proposed in the work of [HCB05] and [CHB09]. In IGA, the description of a structure’s geometry and its analysis under loading are carried out using a common set of shape functions, namely NURBS, which eliminates the obstacles and inaccuracies associated with the translation from geometry to analysis and vice-versa. This integration of geometry and analysis not only makes the analysis be capable of utilizing the excellent properties of NURBS, but also naturally powers the development of structural shape optimization. With the ability to preserve exact CAD geometrical descriptions and a enhanced sensitivity analysis as explained in [WFC08] and [CH09], isogeometric shape optimization has attracted significant attention. Some recent contributions include the design optimization of curved beams in [NAG10a,NA11], membranes in [MEGG11] and shells in [NIA13,KSWB14]. Shape optimization using a T-spline-based IGA method is presented in [HCC10] and further extended to topological optimization in [SKY10]. Some other contributions under the IGA framework include a variational formulation for stress constraints in [NAG10b], a formulation for the boundary element method in [LQ11], a sensitivity analysis with respect to both the positions and the weights of NURBS control points in [Qia10], a transformed-basis function to impose essential boundary conditions in [KYC13] and a regularization method for design sensitivities in [AFA13]. Non-structural applications have been analyzed in [QS11] for photonic crystals, in [NG13] for fluid mechanics and in [YHC13] within the context of heat conduction. A method to avoid mesh irregularity for the interior control points updating in isogeometric shape optimization is presented in [CC15].

1.1.2.

S

-

The above-mentioned contributions reflect the growing interest in shape optimiza-tion within the IGA framework, however, this research effort has been limited to static or steady-state loading conditions. One goal of the research presented here is to adapt and expand the formulation of shape optimization for time-dependent processes within the context of IGA. Previous work on structural design optimization under time-dependent load can be found in the context of sizing and topology

optimization for beams, trusses, dynamic systems and elastic solid structures. Using

a modal analysis, a methodology for solving sizing problems for trusses is developed in [FAH77] and extended in [HA84] to account for point-wise constraints. In [LA87], spring and damping forces are used as variables for design purposes. Using the equivalent static loads, the effect of dynamics loads have been analyzed in [CP02], [PLP05] and [KCP01]. Previous contributions related to inverse shape design work for transient heat conduction, which have relied mostly on a classical FEM anal-ysis, can be found in the work of Dulikravich [Dul88], Jarny et al. [JOB91], Huang and Chaing [HC08] and Korycki [Kor07]. Shape sensitivity analysis for the linear or nonlinear transient heat conduction problems can be found in the works of Haftka [Haf81], Tortorelli et al. [TH89,THL89], Słu˙zalec et al. [SK92,KS96], Dems et

al.[DR99a,DR99b], Gu et al. [GCZG02] and Korycki [Kor06]. Shape optimization work of thermoelastic problems with transient thermal field can be found in the work of Kane et al. [KKS91], Gao et al. [GG99] and Song et al. [SSM04]. Other related design sensitivity analysis may be found in[HG86], [VKHK05] and [CK05]. Within the context

1.1.BACKGROUND

1

of topology optimization, a fully-explicit dynamic problem in an elastic medium is considered in [Tur02]. While a comprehensive review of the research carried out on this topic is outside the scope of the present work, it is remarkable that shape optimization for time-dependent processes has received less attention despite its clear relevance in practical applications where variable loads are the norm rather than the exception.

1.1.3.

A

An important issue in structural shape optimization is design sensitivity analysis. Some commonly-used methods and approaches of design sensitivity analysis can be found in [CK05] and [VKHK05]. Recent contributions on design sensitivity analysis for static problems based on the concept of isogeometric analysis have been presented in [Qia10], [KYC13] and [AFA13] (see also the above-mentioned the references on isogeometric shape optimization). It is well-accepted that the adjoint method is highly efficient and accurate for problems involving a relatively large number of design variables. The general theory of design sensitivity analysis for structural problems using the adjoint method may be found in, among others, [CH83], [DM84] and [CC94].

Compared with shape optimization based on the traditional finite element method, isogeometric shape optimization can provide a simpler and more accurate sensitivity, especially when the sensitivity expressions depend on geometric properties such as curvature. This advantage is discussed in terms of enhanced sensitivity in the work of [CH09]. In the framework of NURBS discretization for design, full analytical sensitivities with respect to both the positions and the weights of NURBS control points are achievable ([Qia10]). Design sensitivity analysis is further studied using transformed basis functions for Kronecker delta property in the work of [KYC13].

1.1.4.

S

In isogeometric shape optimization, the control points, or some value associated with the control points, are commonly chosen as the design variables. The shape design sensitivity leads to the calculation of the discrete shape gradient with respect to control point variables. Depending on how a search direction is constructed from the discrete shape gradient, discretization - dependency of the optimization history might arise ([KSWB14]). This discretization - dependency can affect the convergence speed and may lead the optimization process into a sub-optimal solution. It should be noted here that this discretization - dependency is not limited to shape optimiza-tion based on isogeometric analysis; it also occurs in shape optimizaoptimiza-tion based on other analysis methods, such as the classical finite element formulation.

The discretization - dependency of the sensitivity analysis requires some tech-nique to find a reasonable search direction. In [NAG10a], a Sobolev semi-norm, referred to as "shape change norm", is introduced to balance the shape variation with the cost of constructing the Sobolev semi-norm and solving a system of equations. In [AFA13], the H1gradient method is used such that the discretization-dependency is avoided at the expense of solving a reshaping problem constructed by the H1 gradient method. In [IA04], [AT06], [STA07] and [RFS+14], a similar method called

1

"traction method" is used to secure the shape regularity at about the same price as the H1gradient method. Other approaches and mesh regularization strategies, such as the filtering method, are presented for discretization-free shape optimization in [BFLW10], [LBT11] and [FWB13] to avoid the shape irregularity, referred to as a “parameterization-free" optimization, however, all of these approaches require a system of linear equations to be constructed and solved. In [KSWB14], a "sensitiv-ity weighting" scheme is employed to obtain a discretization-independent search direction. The "sensitivity weighting" is computationally much simpler since it only requires a locally integrated ’weighting’ factor over the local support of the corresponding design control point. While the "sensitivity weighting" approach is shown to work in [KSWB14], the underlying reasons of its success and where the problem originates remain unclear.

The simplest and most intuitive search direction is the steepest descent direction. The use of the steepest descent direction, computed from the discrete shape gradient, makes the optimization history strongly dependent on the discretization. In this thesis, we analyzed the cause of the discretization - dependency of the steepest de-scent search direction from the underlying mathematical aspects and propose three normalization approaches to obtain a discretization - independent search direction, namely (i) a “standard” normalization, (ii) a diagonally-lumped mapping matrix (DLMM) normalization and (iii) a simplified DLMM normalization. The approaches proposed in the present contribution will also work for the finite element method based shape optimization, especially those methods using NURBS parameterization to describe the geometry (see, e.g., [BF84], [ZWY10], [ELA11], [WZ12] and [CZZG14]).

1.2.

S

COPE AND AIMS

The aim of the research discussed in this thesis is to investigate the shape opti-mization for some time-dependent problems under the isogeometric framework using continuous adjoint method. Isogeometric shape optimization for quasi-static mechanical problems and transient heat conduction problems is studied, and during the development of the optimization framework, some problems such as functional discontinuities, search directions for the discrete models and the mean value property of B-spline basis function, etc. are also studied. The key contribution of this work are listed below:

• Isogeometric shape optimization framework for quasi-static mechanical prob-lems using continuous adjoint shape sensitivity analysis is presented.

• Isogeometric shape optimization framework for transient heat conduction problems using continuous adjoint shape sensitivity analysis is presented. • Continuous shape sensitivity analysis considering functional discontinuities is

studied.

• The link between the search directions for continuous and discrete formula-tions in isogeometric shape optimization and the normalization approaches is analyzed.

1.3.LAYOUT AND NOTATION SCHEME

1

• The mean value property of B-spline basis function is derived using mathemat-ical induction method.

It is worth mentioning that, the shape optimization framework presented in this work is also applicable to the finite element (FE) based shape optimization and the continuous adjoint shape sensitivity analysis is also applicable to the steady or static problems with some modifications on the time-dependent terms.

1.3.

L

AYOUT AND NOTATION SCHEME

The structure of this work is as follows.

1. The general formulations of NURBS is presented in Chapter2, in which the mean value property of B-spline basis function is proposed and subsequently proved using mathematical induction method. Then, the general formulations of isogeometric analysis for quasi-static mechanical and transient heat con-duction problems are derived.

2. The continuous design parameterization of the shape design problem is pre-sented in Chapter3. Based on this, the transport relations considering the func-tional discontinuities for continuous shape design derivatives are proposed and derived.

3. The continuous and discrete search directions are presented in Sec.4.1of Chap-ter4. The discretization - dependency of the un-normalized search direction is illustrated in Sec.4.2using a volume minimization problem with different NURBS discretizations. A quadratic norm induced by the NURBS discretization is derived in Sec.4.3. Subsequently, a standard normalization approach for an optimization problem with a quadratic norm is proposed to obtain a discretiza-tion - independent search direcdiscretiza-tion. Using the concept of the lumped mass matrix and the partition of unity property of NURBS, two simpler approaches are also proposed. The suitability of the normalization schemes is verified in Sec.4.4using the same volume minimization problem presented in Sec.4.2and, additionally, using the sensitivity analysis in a volume minimization problem of a two-dimensional domain. Design problems of a three-dimensional fillet and a two-dimensional thermal isolating panel are presented in Sec.4.5 to demonstrate the superiority of the proposed normalization approaches. Some discussions and remarks are presented in Sec.4.6.

4. The shape optimization for quasi-static problems in a continuous setting is presented in Sec.5.1 of Chapter 5. Subsequently, the continuous adjoint sensitivity analysis is developed in Sec.5.2. Necessary details are included to treat the discontinuity of characteristic functions. The structural problem and its sensitivity are subsequently discretized in Sec.5.3within the context of IGA. The algorithm used to solve shape optimization problems for quasi-static problems numerically is presented in Sec.5.4. Several benchmark problems are presented in Sec.5.5 to verify the proposed methodology. An optimal

1

control-like problem is presented in Sec.5.6to illustrate the range of possible applications offered by the method within the context of time-varying loading. 5. The general formulation of the shape optimal design and passive control considering the jump conditions in the objective functional for transient heat conduction problems is presented in a continuous setting in Sec.6.1(i.e., con-tinuous description of the design and analysis spaces). By assigning a Lagrange multiplier for the strong form of the governing equations at every material point during the transient time interval, the continuous adjoint sensitivity analysis is developed in Sec.6.2. Necessary details are included to treat the discontinuity of characteristic functions that represent the regions (in space and/or time) where the objective functional is defined. The isogeometric design discretization of the sensitivity is subsequently discretized in Sec.6.3 within the context of IGA. This section also includes the algorithm used to numerically solve shape optimization problems. Some numerical problem including the shape design of a plunger for a molten glass forming die and a thermal protection system (TPS) for a re-entry ballistic vehicle nose are presented in Sec.6.4to illustrate the range of possible applications.

6. Finally, the conclusions of the work presented in the thesis and recommenda-tions for future work are presented in Chapter7.

The notation scheme in this dissertation is as follows. Lightface italic Latin or Greek letters, e.g., a and α, are used to denote scalars or the elements of a vector. Boldface italic lowercase Latin or Greek letters, e.g., a and α, are used to denote vectors. Boldface italic capital Latin letters, e.g., A, in most cases, are used to denote to the second-order tensors. For the B-spline or NURBS basis functions, the boldface italic capital letters, N and R, are used for the vectors with the purpose of being consistent to the customary use. Blackboard bold capital letters, e.g., A, are used for the fourth-order tensor. The dot (inner) product of two vectors a and b and two second-order tensors A and B are denoted aTb or a · b and A · B , respectively. The

product of a second-order tensor A and a vector b is denoted Ab. The product of two second-order tensors A and B is denoted AB . The tensor product and cross product of two vectors a and b are denoted a ⊗ b and a × b, respectively.

2

I

SOGEOMETRIC ANALYSIS

Although NURBS-based geometrical objects have been used in the computer-aided design community for several decades, it is only relatively recently that they have gained wider acceptance as a tool for analysis. In this chapter, the conception, basic properties and geometry operations of NURBS are briefly introduced. Based on this, the isogeometric analysis for quasi-static and transient problems is presented. It is worth noting that in this chapter, the mean value property of the B-spline basis function is proposed and proved using mathematical induction method. This mean value property is useful for the normalization approaches of the steepest search direction presented in Chapter4.

2.1.

NURBS

2.1.1.

B-

A B-spline curve with a degree of p is defined as

x[ξ] =

n

X

i=1

Npi[ξ]xi 0 É ξ É 1 (2.1)

where x = [x1, x2, · · · , xd] is the location of a point embedded in d-dimensional

physical space corresponding to the parameter ξ defined in a parametric space, n is the number of the control points, Ni

p is the i th basis function of degree p, and xi =[x₁i, xi₂, · · · , xi_d] is the location of the i th control point. The i th B-spline basis function of degree p can be defined as

N0i[ξ] := ( 1 if ξiÉξ < ξi+1 0 otherwise Npi[ξ] := ξ − ξi ξi+p−ξi N_p−1i [ξ] + ξi+p+1−ξ ξi+p+1−ξi+1 N_p−1i+1[ξ], (for p > 0) . (2.2) 7

2

where ξiis the i th element of a non-decreasing knot vector, i.e., ξ ={ξ1, ξ2, ξ3, ..., ξn+p+1}.

In general, the knot vector can be any non-decreasing real vectors, but it is convenient to normalize it such that 0 ≤ ξi ≤1 and ξ ∈ [0, 1]. The knot vector is said to be

uniform if the knots are equally spaced and open if the knots at each ends have p + 1

multiplicity, respectively. A uniform vector with p + 1 multiple equal knots at each end is referred to as open-uniform. In particular, for a given polynomial degree p, the so-called open knot vector has the following form:

ξ ={0, ..., 0 | {z } p+1 , ξp+2, ..., ξn, 1, ..., 1 | {z } p+1 }. (2.3)

The B-spline basis functions are plotted in Fig.2.1(a) for a open-uniform know vector {0, 0.25, 0.5, 0.75, 1} with a degree of 0, Fig.2.1(b) for a open-uniform know vector {0, 0, 0.25, 0.5, 0.75, 1, 1} with a degree of 1, and Fig.2.1(c) for a open-uniform know vector {0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1} with a degree of 2.

The multiplicity of the knots has important consequences on the smoothness of the curve being described. For a open knot vector, the internal knots, ξi, p +2 ≤ i ≤ n,

have a multiplicity of k, 1 ≤ k ≤ p, and the curve at ξ = ξi is Cp−k-continuous, e.g. see

Fig.2.2. The NURBS basis function Npi[ξ] is nonnegative and is nonzero only on the

subinterval [ξi, ξi+p+1]. This property is referred to as local support, which implies

that moving a control point xi only affects the part of the geometry that corresponds to the subinterval [ξi, ξp+i+1] in the parametric space.

2.1.2.

NURBS

In principle, it is possible to construct a curve based directly on B-splines. However, in order to increase the range of geometrical shapes that can be described with a given knot vector, it is customary in computer-aided design software to work with NURBS, which are a generalization of B-splines. For a given partition (knot vector), suppose that the B-splines Npi (of order p) have been constructed, with i = 1, . . . , n. Introduce a

positive scalar wifor each function i , termed the weighting factor of the i -th function

and i -th control point, and define a NURBS function Ripas

Rip[ξ] := Npi[ξ]wi w[ξ] , (2.4) where w[ξ] := n X j =1 Npj[ξ]wj.

A (one-dimensional) NURBS-based curve ˆxc, embedded in a d-dimensional

eu-clidean space, can be constructed as follows:

xc[ξ] :=

n

X

i=1

2.1.NURBS

2

Figure 2.1: B-spline basis functions of (a) a open-uniform know vector {0, 0.25, 0.5, 0.75, 1} with a degree of 0, (b) a open-uniform know vector {0, 0, 0.25, 0.5, 0.75, 1, 1} with a degree of 1, and (c) a open-uniform know vector {0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1} with a degree of 2.

0 0.25 0.5 1 0 1 _{R 1,2} R 2,2 R 3,2 R 4,2 R 5,2 R 6,2

Figure 2.2: B-spline basis functions of a open-uniform know vector {0, 0, 0, 0.25, 0.5, 0.5, 1, 1, 1} with a degree of 2 and multiple knots at 0.5.

2

where ˆxc[ξ] represents a d-dimensional position vector of a point in the curve for a

given parameter ξ ∈ [0, 1]. The vector-valued term xi in Eq.2.5is the position vector of the i -th control point, with the index i ranging throughout all the distinct control points. The position vectors of points on the curve as well as the control points are referred to an arbitrarily-chosen origin using a given orthonormal basis {e1, . . . , ed} in

the d-dimensional space.

NURBS-based descriptions of higher-dimensional geometrical objects, namely surfaces and volumes, may be obtained from a partition involving two scalar pa-rameters for surfaces and three for volumes. Denote as ξ and η the papa-rameters for surfaces and as ξ, η and ζ the parameters for volumes. In this context the partition of the parametric space can be achieved with vectors ξ = {ξ1, . . . , ξn+p+1} and η =

{η1, . . . , ηm+q+1} for surfaces and, additionally, with ζ = {ζ1, . . . , ζl +r +1} for volumes.

Making use of this partition, surfaces and volumes can be represented, respectively, as xs[ξ, η] := n X i=1 m X j =1 Ri,jp,q[ξ, η]xi,j , (2.6) and xv[ξ, η, ζ] := n X i=1 m X j =1 l X k=1 Rp,q,ri,j ,k[ξ, η, ζ]xi,j ,k , (2.7)

where xi,jin Eq.2.6and xi,j ,kin Eq.2.7represent the d-dimensional position vectors of the control points. The NURBS basis functions Ri,jp,q for surfaces and R

i,j ,k

p,q,r for

volumes are obtained as

Ri,jp,q[ξ, η] := Npi[ξ]M j q[η]wi,j Pn i′₌₁ Pm j′₌₁Ni ′ p[ξ]M j′ q[η]wi ′_,j′ . (2.8) and Ri,j ,kp,q,r[ξ, η, ζ] := Npi[ξ]M j q[η]Lkr[ζ]wi,j ,k Pn i′₌₁ Pm j′₌₁ Pl k′₌₁Ni ′ p[ξ]M j′ q[η]Lk ′ r [ζ]wi ′_,j′_,k′ , (2.9)

where the polynomial degrees of the B-splines in the parametric axes ξ, η and ζ are, respectively, p, q and r . To simplify the notation, introduce the following mapping function: I =     

I [i ], for 1D parametric space;

I [i , j ], for 2D parametric space;

I [i , j , k], for 3D parametric space;

(2.10)

This uniquely-determined index I is taken to range from 1 to St, where Stis the total

number of control points that characterize the geometry. Correspondingly, NURBS-based geometrical objects can be generally expressed as

ˆ x := St X I =1 RIxI, (2.11)

2.1.NURBS

2

with the scalar RI denoting the NURBS basis function associated with the I -th parametric control point and xI representing the corresponding vector-valued co-efficient. For simplicity, the arguments from the parametric space and the individual polynomial powers used have been suppressed.

2.1.3.

B

For the analysis and design purpose, some of the most important basic properties of NURBS basis functions are listed here.

P1 Non-negativity: Rip[ξ] ≥ 0 for all i , p, and ξ ∈ [0, 1].

P2 Partition of unity:Pn_i=1Rpi[ξ] = 1 for all ξ ∈ [0, 1].

P3 Local support and approximation: Rpi[ξ] = 0 for ξ ∉ [ξi, ξp+i+1]. Moving control

point xi or changing the weight wi affects the curve only on the portion corresponding to the interval [ξi, ξp+i+1].

P4 The B-spline basis functions, Npi[ξ], are special cases of the NURBS basis

functions, Rip[ξ], with wi=1 for all i .

P5 The derivative of a B-spline basis of degree p > 0 is a function of the basis functions of degree p − 1, i.e.,

Np′i= p ξi+p−ξi N_p−1i [ξ] − p ξi+p+1−ξi+1 N_p−1i+1[ξ], (p > 0) _(2.12)

Detailed derivation or proof for these properties can be found in [PT97]

2.1.4.

M

B-

LOCAL SUPPORT

In section Sec.4.3, the following property of B-spline basis is employed to obtain a simplified normalization approach for the search directions in isogeometric shape optimization: R Npidξ ξi+p+1−ξi = 1 p + 1. (2.13)

This formulation is termed mean value property of a B-spline basis function over its local support. In other words, given a B-spline with a degree of p, the mean value of every basis function over its local support is constant and equals to_p+11 .

We can prove this by induction on p. For p = 0, N₀i equals to 1 in its local support interval. It is obvious that the mean value of a basis over its local support satisfies Eq.2.13. For p = 1, N₁i is a linear function of ξ in its local support interval. It is also obvious that it satisfies Eq.2.13. Now assume that Eq.2.13is true for p − 1, p > 1. By

2

using Eq.2.2and Eq.2.12, we have Z Npidξ = Z µ ξ − ξi ξi+p−ξi Np−1i + ξi+p+1−ξ ξi+p+1−ξi+1 Np−1i+1 ¶ dξ = −ξi Z _Ni p−1 ξi+p−ξi dξ + ξi+p+1 Z _Ni+1 p−1 ξi+p+1−ξi+1 dξ + Z Ã _ξNi p−1 ξi+p−ξi + −ξN_p−1i+1 ξi+p+1−ξi+1 ! dξ =ξi+p+1−ξi p + 1 p Z Ã_pξNi p−1 ξi+p−ξi − pξN i+1 p−1 ξi+p+1−ξi+1 ! dξ =ξi+p+1−ξi p + 1 p Z ξN_p′idξ (2.14)

Using the product rule, (ξNpi)′=Npi+ξNp′i , it can be obtained that

Z ξNp′idξ = Z (ξNpi)′dξ − Z Npidξ (2.15)

If Npi[ξ] is smooth, it is true that ξNpi[ξ] is also smooth. Thus,

Z

(ξN_pi)′dξ = (ξN_pi)|ξ=ξ_ξ=ξn+p+1

1 =0. (2.16)

For the case in which Npi[ξ] is piece-wisely smooth, it still holds that

Zξn ξ1 (ξN_pi)′dξ = 0, since Zξn ξ1 (ξN_pi)′dξ = (ξN_pi)|ξ=ξ∗− ξ=ξ1+(ξN i p)| ξ=ξn ξ=ξ∗ + =0, (2.17)

in which ξ = ξ∗_{corresponds to the location where discontinuity occurs. This is also}

applicable to the basis functions corresponding to two end knots in an open knot

vector. Therefore, _Z

ξNp′idξ = −

Z

Npidξ. (2.18)

Substitute Eq.2.18into Eq.2.14, we can obtain that Z

Npidξ =

ξi+p+1−ξi

p + 1 . (2.19)

Recalling the fact that Npi[ξ] is nonzero only on the subinterval [ξi, ξi+p+1], the proof

can be completed. It should be noted that a similar property, termed unity of integral by Schoenberg in [CS66] for the fundamental spline functions, was presented as

Z∞ −∞

2.1.NURBS

2

in which the fundamental spline functions Mp+1[ξ; ξi] are shown to be frequency

functions with a degree of p (or an order of p + 1). Equation (2.20) is obtained by applying a theorem of Peano (see [Dav75]) to the divided differences used in the fundamental spline function definition. The B-spline basis functions defined using

Cox-de Boor recursion formula are the fundamental functions Mp+1normalized by a

coefficient ξi +p+1_p+1−ξi [Cox72,DB72,Sch07]. In view of this, and taking into account the local support, Eq.2.20is identical to Eq.2.13.

For the B-spline basis functions defined in 3D parametric space with parameters

ξ, η and ζ and the degrees of p, q and r , respectively, since the integrands Npi[ξ],

Mqj[η] and Lkr[ζ] are independent over the integral, based on Eq. 2.19, it can be

obtained easily that Zζk+r +1 ζk Zηj +q+1 ηj Zξi +p+1 ξi Npi[ξ]M j q[η]Lkr[ζ]dξdηdζ =ξi+p+1−ξi p + 1 Zζk+r +1 ζk Zηj +q+1 ηj Mqj[η]Lkr[ζ]dηdζ =µ ξi+p+1−ξi p + 1 ¶ µ ηj +q+1−ξj q + 1 ¶ Zζk+r +1 ζk Lkr[ζ]dζ =µ ξi+p+1−ξi p + 1 ¶ µ ηj +q+1−ηj q + 1 ¶ µ ζk+r +1−ζk k + 1 ¶ . (2.21)

Therefore, the mean value of the basis over its local support in 3D parametric space is Rζk+r +1 ζk Rηj +q+1 ηj Rξi +p+1 ξi N i p[ξ]M j q[η]Lkr[ζ]dξdηdζ ¡ ξi+p+1−ξi ¢ ¡ ηj +q+1−ηj ¢ (ζk+r +1−ζk) =_¡ 1 p + 1¢ ¡q + 1¢(r + 1). (2.22)

Similarly, the mean value in 2D parametric space can be obtained as Rηj +q+1 ηj Rξi +p+1 ξi N i p[ξ]M j q[η]dξdη ¡ ξi+p+1−ξi¢ ¡ηj +q+1−ηj¢ =_¡ 1 p + 1¢ ¡q + 1¢ . (2.23)

2.1.5.

G

In order to have more flexibilities for the geometry description or to obtain a better accuracy for analysis, it is always required to insert some extra control points, basis functions or knots into a given NURBS parameterization without losing any geometry information. This operation can be achieved by knot insertion or knot refinement.

KNOT INSERTION

Knot insertion is to add a single knot, or multiple knots sequentially, into a give knot

vector such that the geometry will be represented with the extended knot vector. Suppose that a NURBS curve xc[ξ] =Pn_i=1Rip[ξ]xi is defined over a knot vector ξ ={ξ1, ξ2, ξ3, · · · ξn+p+1}. Inserting a knot ˆξ ∈ [ξa, ξa+1) gives a new knot vector

2

¯

ξ ={ ¯ξ1=ξ1, · · · , ¯ξa=ξa, ¯ξa+1= ˆξ, ¯ξa+2=ξa+1, · · · , ¯ξn+1+p+1=ξn+p+1}. This new

knot vector ¯ξ, together with the corresponding basis functions ¯Ri_p[ξ], can be used to represent the curve as

xc[ξ] = n+1 X i=1 ¯ Rpi[ξ] ¯xi, (2.24)

where ¯xi correspond to a series of new control points and can be obtained through ¯ xi=αix¯i+(1 − αi) ¯xi−1, (2.25) in which αi=        1, i ≤ a − p ˆ ξ−ξi ξi +p−ξi, a − p + 1 ≤ i < a 0, i ≥ a (2.26)

It should be noted here that the above approach is essentially the same as the approach presented in [PT97] except that the index used here starts from 1 while the index used in [PT97] starts from 0.

KNOT REFINEMENT

Inserting multiple knots using knot insertion is not efficient. It is possible to insert multiple knots at once in a more efficient way. This approach is termed as knot

refinement.

Given a NURBS curve xc[ξ] =Pn_i=1Rip[ξ]xidefined over a knot vector ξ = {ξ1, ξ2, ξ3,

· · ·, ξn+p+1}. Inserting a vector of knots ˆξ = ˆξ1, ˆξ2, · · · , ˆξr, with ξa≤ ˆξi≤ ˆξi+1≤ξbfor

all i , gives a new knot vector ¯ξwith ξi ∈ ¯ξand ˆξi ∈ ¯ξ. The curve will be represented

using the new knot vector ¯ξ, the corresponding basis functions ¯Rip[ξ] and control

points ¯xi, i = 1, 2, · · · , n + r + p + 1. The knot refinement approach can be performed by the following steps:

a find the indices a and b for ξaand ξb

b evaluate the unchanged control points: ¯xi =xi, for i ≤ a − p and ¯xi+r =xi for

i ≥ b − 1

c compute the location of the r + p + b − a − 2 unknown control points using Eq.2.25and Eq.2.26with a forward approach starting at ¯xa−p+1or a backward approach starting at ¯xb+r −1

2.2.

NURBS-

BASED

I

SOGEOMETRIC ANALYSIS

The properties of the NURBS basis functions enable the possibilities of utilizing the NURBS basis as the shape functions of the finite element analysis, which is the underlying principle behind isogeometric analysis. In the isogeometric analysis, the unknown fields, such as the displacement in the mechanical problems or the temperature in the heat conduction problems, are approximated using the same type of representation as for the geometry. Beyond this, the isogeometric approach is

2.3.ISOGEOMETRIC ANALYSIS FOR QUASI-STATIC MECHANICAL PROBLEMS

2

similar to the isoparametric approach in the traditional finite element formulation. In general, relation Eq.2.11is used to define the geometry, which is often referred to as an “exact” geometry. The NURBS representation of the unknown fields, together with the corresponding approximations for gradient fields, are substituted in a weak formulation of partial differential problem in order to obtain a system of equations for the discrete unknowns. The procedure is analogous to the Bubnov-Galerkin method used in the traditional finite element method that uses Lagrange polynomials as shape functions. One notable difference, however, is that NURBS are generally not interpolatory, hence the imposition of essential boundary conditions requires special attention(see, e.g., [RSB+13],[GR15]). In particular, for prescribed non-zero essential boundary conditions, one cannot directly determine the so-called lift function, unless the prescribed boundary condition are defined a priori in terms of a NURBS-based representation.

The term “exact" in this context only implies that Eq.2.11is used as a definition and not, for example, as an approximation of an existing geometry. In contrast, the numerical analysis of a structural design defined by Eq.2.11in general provides only an approximation of the (exact) solution. In practice, the first step to construct the basis functions is to choose a partition of the parametric space, a process that is analogous to generating a mesh in the traditional finite-element method. The basis functions are subsequently constructed from a general formula applied to the specific partition of the parametric space. It is possible to insert knots (points, lines or surfaces) without affecting the geometry, which is a useful property to improve the accuracy of the numerical analysis if required.

2.3.

I

SOGEOMETRIC ANALYSIS FOR QUASI

-

STATIC MECHANI

-CAL PROBLEMS

2.3.1.

Q

-

In solid mechanics, a quasi-static process is a loading process that is instantaneously dependent on time but infinitely slowly performed such that the inertial effects are negligible. Consider a structural domain Ω with boundary Γ, as is shown in Fig.2.3. The loading functions ˆu and ˆt are assumed to vary slowly in time such that the

corresponding structural problem may be treated as a quasi-static process, i.e., the inertial term in the balance of linear momentum is neglected. Correspondingly, at each time t ∈ [0, T ], an equilibrium problem is solved. Furthermore, it is assumed that the structural component is made out of a homogeneous, linearly elastic material characterized by a stiffness tensor C. Consequently, the quasi-static structural problem is formulated as follows:

                 l [u] := div σ + f = 0 , σ = σT in Ω × [0,T ] ǫ =1/2¡∇u + (∇u)T¢ in Ω × [0,T ] σ = Cǫ in Ω × [0,T ] t = ˆt , t = σn on Γt×[0, T ] u = ˆu on Γu×[0, T ] (2.27)

2

p

Ω

x(p)

(t)

(t)

Γ

Γ

Γ

Figure 2.3: A structure occupied in domain Ω with boundary Γ under quasi-static loading.

where σ = σ[x, t ] is the stress tensor, ǫ = ǫ[x, t ] is the strain tensor, f = f [x, t ] is the body force, u = u[x, t ] is the displacement and n = n[x] is the outward normal unit vector on the boundary. The (fourth-order) tensor C is assumed to have major and minor symmetries. In particular, the major symmetry corresponds to CT = C, where the (major) transpose corresponds in Cartesian components to C_{i j kl}T =Ckl i j. Observe

that in Eq.2.27the divergence and the gradient operators, div and ∇, are taken with respect to position x.

2.3.2.

W

Multiplied by a smooth test function, ¯u, and integrated by parts, a weak formulation

of the problem Eq.2.271can be obtained by applying the divergence theorem and

other conditions in Eq.2.272∼5as follows:

〈l [u], ¯u〉Ω= Z Ω (−C∇2u · ¯u − f · ¯u) dΩ = Z Ω C∇u · ∇ ¯u dΩ − Z Ω f · ¯u dΩ − Z Γ t · ¯u dΓ = 0 , (2.28)

in which 〈·, ·〉Ωis a inner product integrated over Ω and the following relation holds

Z Γ t · ¯u dΓ = Z Γˆt ˆt · ¯u dΓ + Z Γuˆ t · ¯u dΓ. (2.29)

2.3.3.

I

The basic idea behind isogeometric analysis is to use the NURBS basis functions as the shape functions for interpolation of the displacement field. Under this framework, the displacement field can be discretized using NURBS parameterization as follows:

u =X

I

RIuI=RU , (2.30)

where uI is the associated displacement at the I th control point, U is the vector of nodal displacement

2.3.ISOGEOMETRIC ANALYSIS FOR QUASI-STATIC MECHANICAL PROBLEMS

2

Correspondingly, the displacement gradient with the above interpolation can be expressed as                                      ∂u1 ∂x1 ∂u1 ∂x2 ∂u1 ∂x3 ∂u2 ∂x1 ∂u2 ∂x2 ∂u2 ∂x3 ∂u3 ∂x1 ∂u3 ∂x2 ∂u3 ∂x3                                      =                     ∂R1 ∂x1 0 0 ∂R2 ∂x1 0 0 · · · ∂R1 ∂x2 0 0 ∂R2 ∂x2 0 0 · · · ∂R1 ∂x3 0 0 ∂R2 ∂x3 0 0 · · · 0 ∂R_∂x1 1 0 0 ∂R2 ∂x1 0 · · · 0 ∂R_∂x1 2 0 0 ∂R2 ∂x2 0 · · · 0 ∂R_∂x1 3 0 0 ∂R2 ∂x3 0 · · · 0 0 ∂R_∂x1 1 0 0 ∂R2 ∂x1 · · · 0 0 ∂R_∂x1 2 0 0 ∂R2 ∂x2 · · · 0 0 ∂R_∂x1 3 0 0 ∂R2 ∂x3 · · ·                     U = BU . (2.31)

The representation Eq.2.30, together with the corresponding approximations for the load and strain/stress fields, are substituted in Eq.2.28in order to obtain a system of equations for the unknowns uI:

K U = F (2.32)

where K is the global stiffness matrix and F is the load vector. The global stiffness matrix K can be calculated using

K =

Z

Ω

E BTB dΩ, (2.33)

where E is the coefficients of the material property. The numerical calculation of K follows the standard Gaussian integration rule based on the order of the NURBS used in the analysis (see [HCB05] for more details). The load vector F can be calculated using F = Ft+Ff, (2.34) where Ft= Z Γt RT{ˆt1, ˆt2, ˆt3}TdΩ. and Ff = Z Ω RT{ f1, f2, f3}TdΩ.

The above isogeometric analysis approach is used in Sec.5to solve both the primary and the adjoint quasi-static mechanical problems in the optimization process.

2

p

Ω

x(p)

(t)

(t)

Γ

Γ

Γ

(t)

Γ

Figure 2.4: A structure occupied in domain Ω with boundary Γ under transient heat conduction.

2.4.

ISOGEOMETRIC ANALYSIS FOR TRANSIENT HEAT CON

-DUCTION PROBLEMS

2.4.1.

T

Consider an isotropic, homogeneous and linearly thermally-conducting material that occupies a domain/region Ω with boundary Γ as shown in Fig.2.4. The governing partial differential equation for a transient heat conduction problem in domain Ω is

l£θ[x]¤:= ρc∂θ[x, t ]

∂t −k∇

2

θ[x, t ] − Q[x, t ] = 0, x ∈ Ω, (2.35)

where c is the material heat capacity, ρ is the material density, t is time, θ[x, t ] is a scalar function that represents temperature, Q[x, t ] is the inner heat-generating rate per unit volume, k is the thermal conductivity coefficient, and ∇ is the gradient operator. It is assumed that there are only three types of boundary conditions on the boundary Γ, i.e.

Γ = Γθ+ Γq+ Γe (2.36)

where Γθrepresents portions where the temperature is specified, Γqrepresents

por-tions where the heat flux is specified, and Γeis the convection boundary conditions.

The boundary and initial conditions are as follows:

θ = ˆθ on Γθ

−(k∇θ)Tn = qTn = − ˆq on Γq

−(k∇θ)Tn = qTn = −qe=h(θ − θe) on Γe

θ[x, 0] = θ0[x] x ∈ Ω

(2.37)

where vector q[x, t ] = −k∇θ is the heat flow at x ∈ Ω at time t, scaler q[x, t] is the heat flux at x on Γ, ˆθ[x, t] and ˆq[x, t ] are the specified temperature and heat flux on

Γθand Γq, respectively; n[x] is the unit outward normal vector to the boundary, h is

the convection coefficient, θe[t ] is the convective exchange temperature and θ0[x] is

2.4.ISOGEOMETRIC ANALYSIS FOR TRANSIENT HEAT CONDUCTION PROBLEMS

2

2.4.2.

W

Multiplied by a smooth test function, ¯θ, and integrated by parts, a weak formulation

of the transient heat conduction problem Eq.2.35can be obtained by applying the divergence theorem and the boundary conditions in Eq.2.37as follows:

〈l [θ], ¯θ〉Ω= Z Ω (ρc∂θ ∂tθ − k∇¯ 2 θ ¯θ − Q ¯θ) dΩ = Z Ω ρc∂θ ∂tθdΩ +¯ Z Ω k∇θ · ∇ ¯θ dΩ − Z Ω Q ¯θ dΩ − Z Γ q ¯θ dΓ = 0 , (2.38)

in which 〈·, ·〉Ωis a inner product integrated over Ω and the following relation holds

Z Γ q ¯θ dΓ = Z Γq ˆ q ¯θ dΓ + Z Γθ q ¯θ dΓ + Z Γe qeθ dΓ .¯ (2.39)

2.4.3.

I

Interpolate the temperature field using the NURBS basis functions to get

θ =X

I

RIθI=R · θ, (2.40)

where θI is the associated temperature at the I th control point, θ = {θ1, θ2, · · · }T and R = {R1, R2, · · · }. Correspondingly, the temperature gradient with the above interpolation can be expressed as

     ∂θ ∂x1 ∂θ ∂x2 ∂θ ∂x3      =     ∂R1 ∂x1 ∂R2 ∂x1 · · · ∂R1 ∂x2 ∂R2 ∂x2 · · · ∂R1 ∂x3 ∂R2 ∂x3 · · ·    θ = B θ. (2.41)

Following this, the heat flow fields inside the domain can be formulated as

q = −kB θ. (2.42)

The representation Eq.2.40, together with the corresponding approximations for the heat flux and heat flow fields, are substituted in Eq.2.38to obtain a system of equations for the unknowns θI:

C∂θ

∂t +K θ = f (2.43)

where , C is the global capacitance matrix, K is the global conductance matrix, θ is the vector of nodal temperature, and f is the heat flux vector. The global capacitance matrix can be calculated using

C =

Z

Ω

2

The global conductance matrix K can be calculated using

K = Kc+Ke, (2.45) where Kc= Z Ω kBTB dΩ, and Ke= Z Γe hRTdΩ.

The heat flux vector f can be calculated using

f = fq+fe+fQ, (2.46) where fq= Z Γq ˆ qRTdΩ. fe= Z Γe hθeRTdΩ. and fQ= Z Ω QRTdΩ.

The problem becomes a steady state problem, if the temperature field does not change with time, i.e. ∂θ_∂t =0. By differentiating the time, an approximated equation can be introduced with a parameter β ∈ [0, 1] as

C (θ

t +∆t_−θt

∆t ) + K (βθ

t +∆t_{+(1 − β)θ}t_{) = f . (2.47)}

In particular, β = 0 and 1 corresponds to a explicit method and a fully implicit method, respectively, while β = 0.5 corresponds to a semi-implicit (Crank-Nicolson-type) method. The isogeometric approach can be used to solve the primary and adjoint problems in order to determine the fields required to compute shape deriva-tives.

The above isogeometric analysis approach is used in Sec.6 to solve both the primary and the adjoint transient heat conduction problems in the optimization process.

3

T

RANSPORT REL ATIONS

CONSIDERING DISCONTINUITIES

The continuous shape sensitivity analysis requires the shape derivatives of the objective functionals and the constraint formulations. In this chapter, the mate-rial and spatial design derivatives are introduced based on the continuous design parametrization of the structure. Design transport relations with jump conditions are derived to deal with the functionals that contain discontinuities. The volume transport relation that considers the discontinuities is derived using the well known Reynolds’ transport relation and the divergence theorem. The surface transport relation that can deal with the discontinuities on the surface is derived following the definitions of the material and spatial design derivatives and the surface divergence theorem. These transport relations will be used later for the derivation of the continuous shape derivatives.

3.1.

C

ONTINUOUS DESIGN PARAMETRIZATION

A structure, or a structural component, can be characterized by a region Ωs _that

it occupies in space (see Fig.3.1). The scalar parameter s is used to identify different designs of the same structural component. The distinct regions Ωs_{, with}

corresponding boundaries Γs, are referred to simply as "design s". The region Ω0, corresponding to s = 0, is referred to as the initial (or referential) design. Points p in the region Ω0are referred to as referential design points. A design s is assumed to be generated from the initial design through a one-to-one smooth function ˆx[·; s],

referred to as the design function at design s. The design function maps a referential design point p in Ω0into a design point x in Ωs, i.e., x = ˆx[p; s], as shown in Fig.3.1. For each s, the mapping ˆx[·; s] represents a (re-)design of the referential design.

The design function may be interpreted as a deformation, however this is only a geometrical analogy and not a mechanical one since the function ˆx does not generate

stresses or changes in mass density.

3

p x = x[p; s]

^

Ω0

Ωs

Initial design Design s

Γ_ts

Γ_us Γ_t0

Γ0

u

Figure 3.1: Family of design regions Ωsgenerated through design functions ˆx[·; s].

3.2.

M

ATERIAL AND SPATIAL DESIGN DERIVATIVES

For design purposes, the shape modification, in analogy to the configuration change with time t in continuum mechanics, is assumed to change from a design s into a design s + δs smoothly, such that the function defined in the region or on the boundary is differentiable with respect to the design index s. In order to distinguish the material time derivative, which is with respect to time in continuum mechanics, the material design derivative, denoted in the sequel as ˚(·), is introduced to denote the total derivative of a field with respect to the design index s while holding the (referential) design point p fixed.

For a given vector (or tensor) field h defined on the referential design region Ω0, at a given time t ∈ T = [0,T ] and for designs indices s ∈ [0,∞), the material design derivative is expressed as ˚ h[p, t ; s] := ∂h ∂s[p, t ; s] ¯ ¯ ¯ ¯_p, (3.1)

where the notation emphasizes that the partial derivative is taken at constant p. In particular, the design velocity ν corresponds to the material design derivative of the design function ˆx, i.e.,

ν[p; s] := ˚ˆx[p; s] = ∂ ˆx ∂s[p; s] ¯ ¯ ¯ ¯ p . (3.2)

The design index s may be formally interpreted as a time-like parameter, which is convenient to derive useful formulas, however, it is important to distinguish s from the actual time variable t in a time-dependent process since these two quantities are

independent: observe that the design velocity ν defined in Eq.3.2is independent of the actual time t .

In analogy with the spatial time derivative in continuum mechanics, it is custom-ary to introduce the spatial design derivative, denoted (·)′, which corresponds to a derivative with respect to the design index s while holding the (current) design point

3.2.MATERIAL AND SPATIAL DESIGN DERIVATIVES

3

Ωs, during a time interval T and for designs indices s ∈ [0,∞), the spatial design derivative is expressed as h′[x, t ; s] := ∂h ∂s[x, t ; s] ¯ ¯ ¯ ¯ x , (3.3)

where, as before, the notation emphasizes that the partial derivative is taken at constant x. In this context, the notion of “current" refers to “current design". A smooth vector or tensor field h[x, t ; s], defined in the region Ωsduring a time interval T and for designs s ∈ [0, ∞), can also be expressed as a function defined in the

referential design region Ω0 through the connection x = ˆx[p; s]. By analogy, the

Lagrangian (or referential) design description of a field corresponds to a function of the referential design point p ∈ Ω0and the Eulerian (or current) description of the same field corresponds to a function of the current position x ∈ Ωs of the referential design point p, i.e., x = ˆx[p; s]. In order to simplify the notation, the same symbol h

will be used for the referential and current design descriptions of the function, i.e.,

h[x, t ; s] = h[p, t ; s] with x = ˆx[p; s] .

Correspondingly, again by analogy with time derivatives, the material design deriva-tive of h can be expressed as the sum of the spatial design derivaderiva-tive and a convecderiva-tive term, i.e.,

˚

h = h′+(∇h) ν , (3.4)

where ∇ represents the gradient with respect to x. For notational compactness, in Eq.3.4and in the sequel, the arguments of the functions are suppressed. Unless otherwise noted, the arguments used are [x, t ; s] ∈ Ωs×[0, T ] × [0, ∞).

It is worth mentioning that the spatial design derivative (·)′(which is taken with respect to s at constant x) and the spatial gradient ∇ (which is taken with respect to x at constant s) are operations that commute, i.e.,

(∇h)′= ∇(h′) = ∇h′. (3.5)

In contrast, the material design derivative ˚(·) (which is taken with respect to s at constant p) and the spatial gradient ∇ do not commute; instead, the following relation holds:

˚

(∇h) = (∇h)′+ ∇(∇h)ν = ∇ ˚h − (∇h) (∇ν) . (3.6) While developing the sensitivity analysis, one can in principle work with either the spatial or the material design derivatives (or a combination thereof ) as long as consistency is preserved. In this thesis, the spatial design derivatives are used to develop the sensitivity analysis. For material-design-derivatives-based derivations, one can find them in [CK05] and [WT15].

Henceforth, for a linearly elastic solid, the stress tensor σ is formally viewed as a function of the displacement u upon combination of the constitutive relation, the strain displacement relation and the major symmetry properties of the stiffness tensor C, i.e.,