Delft University of Technology
Shape optimization and optimal control for transient heat conduction problems using an
isogeometric approach
Wang, Z.; Turteltaub, Sergio; Abdalla, Mostafa DOI
10.1016/j.compstruc.2017.02.004
Publication date 2017
Document Version
Accepted author manuscript Published in
Computers & Structures
Citation (APA)
Wang, Z., Turteltaub, S., & Abdalla, M. (2017). Shape optimization and optimal control for transient heat conduction problems using an isogeometric approach. Computers & Structures, 185, 59-74.
https://doi.org/10.1016/j.compstruc.2017.02.004 Important note
To cite this publication, please use the final published version (if applicable). Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
Shape optimization and optimal control for transient heat
conduction problems using an isogeometric approach
Zhen-Pei Wang, Sergio Turteltaub, Mostafa Abdalla
Faculty of Aerospace Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The Netherlands Z.P.Wang@tudelft.nl , S.R.Turteltaub@tudelft.nl, M.M.Abdalla@tudelft.nl
Abstract
This work is concerned with the development of a framework to solve shape optimization prob-lems for transient heat conduction probprob-lems within the context of isogeometric analysis (IGA). A general objective functional is used to accommodate both shape optimization and passive control problems under transient conditions. An adjoint sensitivity analysis, which accounts for possible discontinuities in the objective functional, is performed analytically and subsequently discretized within the context of IGA. The gradient of the objective functional is used in a descent algorithm to solve optimization problems. Numerical examples are presented to validate and demonstrate the capacity to manage thermal fields under transient conditions.
Keywords: Shape optimization, Isogeometric analysis, Transient heat conduction, Passive optimal control, Adjoint method
1. introduction
Engineering devices used for thermal management in applications where the temperature and/or the heat exchanged play an important role, such as heat exchangers, cooling components, heat sinks or thermal protection layers, are typically designed to control the maximum or min-imum temperature that a system should be exposed to or to guarantee a certain heat exchange rate. Given some performance requirements, a judicious choice of materials and shapes is made and the original design is subsequently optimized to improve its performance. This optimization process is typically based on steady-state conditions. Nevertheless, in many technologically-relevant applications, thermal conditions fluctuate during operation. In these cases, an active control system is sometimes used for thermal management. The drawback of an active control system is its additional cost, which may hinder its usefulness. An attractive alternative is passive control, which relies on a suitably-designed system that takes a priori into account fluctuations in the thermal fields. In particular, passive control may be achieved through a shape optimization procedure that finds an optimal shape under transient heat conduction conditions.
Shape optimization has received renewed attention due to development of Isogeometric Anal-ysis (IGA) as an alternative to the (classical) finite-element method (FEM) [1]. The main advan-tage of IGA is that the numerical analysis is carried out with the same shape functions used in most commercially-available Computer-Aided Design (CAD) programs. Consequently, the bot-tleneck of converting a CAD-generated design into an FEM-ready mesh is avoided [2]. CAD
software typically uses non-uniform rational B-splines (NURBS) [3] to represent the design ge-ometry. In fact, prior to the development of IGA, shape optimization has also relied on NURBS for design purposes (see, e.g., Braibant and Fleury [4], Espath et al. [5] and Wang and Zhang [6]). More recently, the integration of IGA and NURBS-based shape optimization has been ex-ploited to improve the overall efficiency and accuracy of the method, particularly in terms of an unified design and analysis parametrization that provides an enhanced sensitivity analysis (see, e.g., Wall et al. [7], Qian [8]).
Isogeometric shape optimization has been used for curved beam structures in Nagy et al. [9] and Nagy et al. [10], vibrating membranes in Manh et al. [11], fluid mechanics in Nørtoft and Gravesen [12], pulsatile ventricular assist devices in Long et al. [13], shells in Nagy et al. [14], Kiendl et al. [15] and Ha [16], photonic crystals in Qian and Sigmund [17], composite fiber orientation in [18, 19], heat conduction problems in Yoon et al. [20], Stokes flow problems in Park et al. [21]. Isogeometric shape optimization using boundary element method is presented in Li and Qian [22] and Lian et al. [23]. The aforementioned contributions illustrate the wide range of applications of the isogeometric shape optimization approach. However, the method has been limited to static or steady-state conditions. Recently, the isogeometric shape optimization method was extended in Wang and Turteltaub [24] to analyze quasi-static processes with slowly-varying external loads.
In the present work, one goal is to extend the isogeometric approach to carry out shape design and passive control of problems governed by a transient behavior, in particular transient heat conduction. Previous contributions related to inverse shape design work for transient heat conduction, which have relied mostly on a classical FEM analysis, can be found in Dulikravich [25], Jarny et al. [26], Huang and Chaing [27] and Korycki [28]. Shape sensitivity analysis
for the linear or nonlinear transient heat conduction problems can be found in Haftka [29],
Tortorelli and Haber [30], Tortorelli et al. [31], Słuz˙alec and Kleiber [32], Kleiber and Sluz˙alec [33], Dems and Rousselet [34, 35], Gu et al. [36] and Korycki [37]. Shape optimization work of thermoelastic problems with transient thermal field can be found in Kane et al. [38], Gao and Grandhi [39] and Song et al. [40]. Other related design sensitivity analysis may be found in Haftka and Grandhi [41], Van Keulen et al. [42], Choi and Kim [43] and Nanthakumar et al.[44, 45].
A second goal in the present work is to extend the design sensitivity for situations where the objective functional is measured only in selected areas of the design region and during selected intervals of the design time period. This modeling feature is useful in cases where only portions of the design region and/or specific time intervals of the analysis period are relevant for the overall shape design. In that context, the sensitivity analysis requires an extended version of the transport theorems to account for discontinuities in the objective functional.
The structure of the present work is as follows: the general formulation of the shape optimal design and passive control considering the jump conditions in the objective functional is pre-sented in a continuous setting in Sec. 2 (i.e., continuous 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 anal-ysis is developed in Sec. 3. Necessary details are included to treat the discontinuities of the characteristic functions that represent the regions (in space and/or time) where the objective functional is defined. The isogeometric analysis and design discretization is discussed briefly in Sec. 4. The design problem and its sensitivity are subsequently discretized in Sec. 5 within the context of IGA. This section also includes the algorithm used to numerically solve shape opti-mization problems. To validate the methodology, two benchmark problems, namely a minimum
p x = x[p; s]
^
Ω0
Ωs
Initial design Design s
Γts Γus Γ0 t Γ0 u
Figure 1: Family of design domains Ωsgenerated through design functions ˆx[·; s].
surface problem and a passive temperature control problem at pre-determined time are presented in Sec. 6.1. To further illustrate the range of possible applications, two examples are presented in Sec. 6.2 and Sec. 6.3, namely 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. Finally, some concluding remarks are given in Sec. 7.
2. Problem statement
2.1. Design function and heat conduction problem
Consider an isotropic, homogeneous and linearly thermally-conducting material that occu-pies a region Ωs with boundary Γs as shown in Fig. 1. The superscript s is a continuous
s-calar parameter that represents different configurations (or states) of the region. Each configu-ration corresponds to a design of the structure. The state of the region with s = 0, i.e. Ω0, is called the referential or initial design. A material design point p ∈ Ω0is mapped to a position
x = ˆx[ p, s] ∈ Ωsby a design function ˆx as indicated in Fig. 1. For simplicity, the mapping ˆx is henceforth also denoted as x and the meaning of the symbol may be inferred from the context (i.e., position or design mapping).
In the domain Ωsand during a time interval T = [0, T ], with T a given final analysis time,
the governing equation for a transient heat conduction problem can be expressed as l [θ[x, t]] := ρc∂θ[x, t]
∂t −k∇
2
θ[x, t] − Q[x, t] = 0, (x, t) ∈ Ωs× T, (1)
where θ = θ[x, t] is the temperature at point x and time t, c > 0 is the heat capacity, ρ > 0 is the mass density, k > 0 is the thermal conductivity, Q = Q[x, t] is the inner heat generation rate per unit volume (volumetric heat supply) and ∇2is the Laplacian operator. For subsequent use, the governing equation is written in terms of an operator l as defined in (1). It is assumed that there are only three types of boundary conditions on the boundary Γs, which is divided as follows:
Γs= Γθs∪ Γqs∪ Γes where Γs
θrepresents portions where the temperature is specified, Γ
s
qcorresponds to portions where the contact heat supply is given and Γs
eis the part of the boundary where heat is exchanged with 3
the environment through convection. The boundary and initial conditions are as follows: θ = ˆθ on Γsθ× T q · n = − ˆq on Γsq× T q · n = −qe= h(θ − θe) on Γes× T θ[x, 0] = θ0[x] in Ωs (2)
with the heat flux vector q given by Fourier’s model, i.e.,
q = −k∇θ .
In (2), ˆθ = ˆθ[x, t] and ˆq = ˆq[x, t] are the specified temperature and contact heat supply on Γs
θ× T
and Γs
q× T, respectively, n = n[x] is the unit outward normal vector to the boundary, h is the convection coefficient, θe=θe[t] is the convective exchange temperature (ambient temperature) and θ0[x] is the initial temperature field in the domain. For notational convenience, the contact heat supply on Γseis denoted as qe. The contact heat supplies ˆq and qeare given as the negative of the normal heat flux, hence they take positive values if heat flows from the external environment into the system and negative values otherwise.
2.2. Transient shape optimization problem
In order to develop a versatile framework for shape optimization problems, which can accom-modate various situations within a unified formulation, a general objective functional is defined such that it can measure the performance on pre-selected parts where the physical phenomenon occurs, both in space and time. This is achieved through the combination of two local per-formance functions, denoted as ψω and ψγ and defined, respectively, in Ωs and Γs, and three
characteristic functions, denoted as ω, γ and ς, that specify, respectively, the (sub)-parts of the region Ωs, boundary Γsand time interval T where the performance is monitored. In particular, the overall objective functional J is defined as
J[s] :=
Z T
0
ς[t]Ψ[t; s]dt with Ψ[t; s] := Ψω[t; s] + Ψγ[t; s] (3)
where Ψωis a functional defined over the region Ωsas
Ψω[t; s] :=
Z
Ωs
ω[ p, t]ψω[θ[x, t; s]] dΩ, (4)
and Ψγis a functional defined on the boundary Γsas
Ψγ[t; s] :=
Z
Γs
γ[ p, t]ψγθ[x, t; s], q[x, t; s] dΓ . (5)
The functions ψω and ψγ, appearing in (4) and (5), are the local performance functions. The
function ψωmeasures the performance in Ωsbased on the temperature field θ while the function
ψγ measures the performance on Γs based on the temperature θ and the contact heat supply
q = −q · n, where n is the outward unit vector on the boundary. The parameter s in the argument of the functions and functionals shown in the equations above indicate that these are evaluated for a given design function ˆx[·, s]. The material design point p used in the argument of the
p x[ p; s+δs] ^ x[p; s] ^ Initial design Γ0 γ Γγs+δs Γ Γ Γγs Ω0 ω Ωs+δsω Ωsω Γγs Ωsω ∂Ωs = ω Γωs ∂Γγs=ϒs T Tς 0 Time intervals T Current design
Figure 2: Illustration of subdomain Ωsω, boundary portion Γγs and time subinterval Tςwhere the
objective functional J is measured.
characteristic functions ω and γ should be interpreted as p = ˆp[x, s], where ˆp is the inverse of the design function ˆx. Following the approach proposed in [43] and [24], the functions ω and γ are associated with, respectively a (sub)domain Ωs
ω⊆ Ωsand a portion of the boundary Γγs ⊆ Γs
in which the performance function is measured.
The functions ω and γ correspond to fixed design points in the reference design domain. Correspondingly, the subdomain Ωωs ⊆ Ωsand the portion Γγs ⊆ Γsare assumed to be obtained
through a mapping of the referential subdomain Ω0ω⊆ Ω0and the referential portion Γ0γ ⊆ Γ0. For
increased design flexibility, the subdomain Ω0
ωand the portion Γ0γmay themselves be prescribed
as functions of time. Hence, the characteristic functions are defined as ω[ p, t] = 1, if p ∈ Ω0 ω(t) 0, otherwise , γ[ p, t] = 1, if p ∈ Γ0 γ(t) 0, otherwise . (6)
Similarly, the characteristic function ς that appears in (3) is associated with a time (sub)-interval Tς⊆ T such that ς[t] = 1, if t ∈ Tς 0, otherwise . (7)
In a typical optimization problem, the design may be subjected to a resource constraint that limits the total volume of the design region, i.e.,
Σ[s] := Z
Ωs
dΩ ≤ Σ0, (8)
where Σ0is a given upper bound. In addition, it is often required to consider constraints of the following form:
ˆ
dc[x, s] ≤ 0 c = 1, . . . , Nc, (9)
where the scalar functions ˆdcrepresent pointwise constraints (such as upper or lower bounds) and
Ncis the total number of constraints. The optimization problem consists on finding an admissible design function ˆx[·, s∗], which satisfies the constraints (8) and (9), such that J[s∗] = min
s J[s].
3. Continuous adjoint sensitivity analysis 3.1. Material and spatial design derivatives
Continuous shape modifications in a shape optimization problem can be understood as a s-mooth sequence of deformations from an initial (referential) design as shown in Fig. 1. The dependence of functions and functionals on the design function ˆx[·, s] is indicated with a depen-dence on the design parameter s. Similar to the notion of material and spatial time derivatives used in continuum mechanics, it is possible to define a material design derivative and a spatial design derivative. The material design derivative (at constant p) is denoted here with either a superimposed open ring, (˚·), or with the operator D(·)/Ds (also known as the total design deriva-tive), while the spatial design derivative (at constant x) is indicated using an apostrophe, (·)′, or
using the operator ∂(·)/∂s (partial design derivative). For functions that only depend on s, the derivative with respect to s is indicated as d(·)/ds. The material and spatial design derivatives of a function h are related as follows:
˚h = Dh Ds = h
′+
(∇h) ν, (10)
where ∇ refers to the gradient with respect to x and ν is the so-called design velocity. The design velocity of a material design point p is the material design derivative of the design function
x = ˆx[ p, s], i.e.
ν := ˚x = Dx
Ds. (11)
The spatial design derivative and the spatial gradient commute, i.e.,
(∇h)′= ∇(h′) = ∇h′. (12)
3.2. Transport relations
For a volume integral of a generic smooth function f defined in Ωs and a characteristic function ω as given in (6), the volume transport theorem is extended as
d ds Z Ωs ω f dΩ = Z Ωs ω f′dΩ + Z Γs ω f νndΓ + Z Ss ω [[ω]] f νndΓ, (13)
where νn = ν · n is the normal design velocity with n the outward normal vector on the boundary,
and Ss
ω = ∂Ωsω− Γs is the interior boundary of the sub-region (see [46]). The notation [[·]] =
(·)+−(·)−indicates the jump of a quantity across a surface/line/point of discontinuity, with (·)±
representing the values on the ± sides according to the normal vector (pointing towards the +side).
For a boundary integral of a generic smooth function f defined on Γs and a characteristic
function γ as given in (6), the boundary transport theorem is extended as d ds Z Γs γ f dΓ = Z Γs γ f′+γ (∇ f · n) νn−κ f νn dΓ + Z Υs [[γ]] f νmdΥ, (14)
where κ is the total curvature of Γs (twice the mean curvature of a two-dimensional surface embedded in three space dimensions), Υs=∂Γs
γis the edge of the boundary Γγs, νm= ν · m, with
m being the normal vector pointing outward of ∂Γs
γand divΓrefers to the surface divergence (see
[47], [48] and [49]).
3.3. Augmented functionals
The optimal design process is to find a design such that the objective functional is minimized (or maximized) while simultaneously satisfying all the design constraints. The governing equa-tion of the transient heat conducequa-tion problem shown in (1) and the associated load and boundary conditions in (2) can be treated as an equality constraint in variational form.
In order to include the constraints in the sensitivity analysis, it is convenient to introduce an augmented objective functional (Lagrangian) defined as
L[s, Λc, ΛΣ] := ˜J[s] + Nc X c=1 hΛc, ˆdc[s]iΩs+ ΛΣ(Σ[s] − Σ0) , (15) where ˜ J[s] = J[s] + Z T 0 hl[θ], ϑiΩsdt. (16)
In the preceding equations, the inner product h·, ·iΩsrefers to an integral over the region Ωsof the
product of the arguments. In (15), the inequality constraints given in (8) and (9) are multiplied by the corresponding Lagrange multipliers Λc = Λc[x] and ΛΣ, respectively. Moreover, using
integration by parts, the transient heat conduction equation (1) is included in (16) as a constraint in the following form:
hl[θ], ϑiΩs = Z Ωs ρc∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ ! dΩ − Z Γs qϑ dΓ = 0 .
The Lagrange multiplier field ϑ = ϑ[x, t] corresponds to the so-called adjoint temperature field. Observe that since the constraint must be satisfied, then the values of J and ˜J, as well as the values of their shape derivatives, essentially coincide. This feature is used in the adjoint method to compute the shape derivative of the functional.
3.4. Shape derivatives of ˜J
In view of (16) and since the variables t and s are independent of each other, the design derivative of ˜Jcan be written as
d ˜J ds = dJ ds + Z T 0 d dshl[θ], ϑiΩsdt, (17)
which consists of two parts: one is the design derivative of the original objective functional and the other one is the design derivative of the equality constraint hl[θ], ϑiΩs = 0. However, in order
to directly compute dJ/ds, it is required to determine the shape derivatives of the temperature, θ′, and the contact heat supply, q′, which typically are computationally expensive to obtain
(us-ing, e.g., a direct perturbation method). Instead, following the adjoint method, the strategy is to manipulate the right hand side of (17) in order to express it in a convenient format (i.e., elimi-nating θ′and q′), thus allowing an efficient computation of d ˜J/ds. In turn, d ˜J/ds delivers the desired expression of dJ/ds without the terms θ′and q′. The steps to convert the right hand side of (17) into a convenient format are indicated in the subsequent sections.
3.4.1. Shape derivatives of the objective functional
Bearing in mind that the contact heat supply on Γesis qe= −h(θ − θe), it follows that
q′e= −hθ′ and ∇qe= −h∇θ on Γse. (18)
Using the transport relations shown in (13) and (14) and in view of (12) and (18), the shape derivative of the objective functional J can be obtained as
dJ ds = Z T 0 ς Z Ωs ωψω,θθ′dΩdt + Z T 0 ς Z Γs ωψωνndΓ + Z Ss ω [[ω]]ψωνndΓ ! dt + Z T 0 ς Z Γs q γψγ,θθ′dΓ + Z Γs θ γψγ,θˆθ′dΓ + Z Γs e γψγ,θθ′dΓ dt + Z T 0 ς Z Γs q γψγ,qqˆ′dΓ + Z Γs θ γψγ,qq′dΓ − Z Γs e γhψγ,qθ′dΓ dt + Z T 0 ς Z Γs γ∇ψγ·n νn−ψγκνn dΓdt + Z T 0 ς Z Υs [[γ]]ψγνmdΥdt . (19)
3.4.2. Shape derivative of the governing equation.
The second term on the right hand side of (16) involves the heat conduction problem, treated as an equality constraint integrated in space and time. In view of obtaining a useful expression for this term, consider first the shape derivative of hl[θ], ϑiΩs, which, using the transport relations
(13) and (14) and in view of (12), can be expressed as d dshl[θ], ϑiΩs = Z Ωs ρc∂θ ′ ∂tϑ + ρc ∂θ ∂tϑ ′−Q′ ϑ − Qϑ′ ! dΩ + Z Ωs k∇θ′· ∇ϑ + k∇θ · ∇ϑ′ dΩ + Z Γs ρc∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ ! νndΓ − Z Γs q ˆ q′ϑ + ˆqϑ′ dΓ − Z Γs θ q′ϑ + qϑ′ dΓ − Z Γs e q′ϑ + qϑ′ dΓ − Z Γs ((∇ (qϑ) · n) νn−(qϑ) κνn) dΓ . (20) Noting that hl[θ], ϑ′i
Ωs = 0 and using (18), it follows from (20) that
d dshl[θ], ϑiΩs= Z Ωs ρc∂θ ′ ∂t ϑ − Q ′ ϑ + k∇θ′· ∇ϑ ! dΩ + Z Γs ρc∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ ! νndΓ − Z Γs q ˆq′ϑ dΓ − Z Γs θ q′ϑ dΓ + Z Γs e hθ′ϑ dΓ − Z Γs ((∇ (qϑ) · n) νn−(qϑ) κνn) dΓ . (21) Since the transient heat conduction equation (1), expressed in strong form as l[θ] = 0, is satisfied for every x ∈ Ωs and every t ∈ T , then the relation (21) holds for every t ∈ T . Integrating the first term in the first integral on the right side of (21) over the time interval T = [0, T ], exchanging the order of integration and subsequently integrating by parts in time provides the following relation:
Z T 0 Z Ωs ρc∂θ ′ ∂tϑdΩdt = Z Ωs ρcθ′ϑt=T t=0dΩ − Z T 0 Z Ωs ρcθ′∂ϑ ∂tdΩdt. (22) 8
It is convenient to introduce the adjoint time variable τ := T − t and the convention that an arbitrary function f may refer to either a function of t or τ, i.e., f = f (t) = f (τ) with τ related to t as indicated (i.e., the proper argument may be inferred from the context). In that case, it holds that Z T 0 f dt = Z T 0 f dτ and Z T 0 ∂ f ∂tdt = − Z T 0 ∂ f ∂τdτ . Following this, (22) becomes
Z T 0 Z Ωs ρc∂θ ′ ∂tϑdΩdt = Z Ωs ρcθ′ϑt=T t=0dΩ + Z T 0 Z Ωs ρcθ′∂ϑ ∂τdΩdτ . (23)
Integrating by parts over Ω the third term in the first integral on the right side of (21) (in order to transfer the gradient from θ to ϑ), one has
Z Ωs k∇θ′· ∇ϑ dΩ = − Z Ωs k∇2ϑθ′dΩ + Z Γs q∗θ′dΓ (24)
where q∗ is the adjoint contact heat supply on the boundary Γs, which is related to the adjoint
heat flux q∗as follows:
q∗:= −q∗·n q∗:= −k∇ϑ . (25)
3.5. Transient adjoint system
Integrating (19) in time and in view of (17), (21), (23) and (24), the design derivative of the augmented functional defined in (16) can be expressed as
d ˜J ds = Φ1+ Φ2, (26) where Φ1:= Z T 0 Z Ωs ρc∂ϑ ∂τ −k∇ 2ϑ + ςωψ ω,θ ! θ′dΩdτ + Z Ωs ρcθ′ϑt=T t=0dΩ + Z T 0 Z Γs q ςγψγ,θ+ q∗ θ′dΓ + Z Γs θ ςγψγ,q−ϑ q′dΓ dτ + Z T 0 Z Γs e h ϑ +ςγψγ,θ h −ςγψγ,q ! + q∗ ! θ′dΓ ! dτ (27) and Φ2:= Z T 0 ς Z Γs ωψωνndΓ + Z Ss ω [[ω]]ψωνndΓ + Z Υs [[γ]]ψγνmdΥ ! dt + Z T 0 Z Γs ρc∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ ! νndΓ ! dt + Z T 0 Z Ωs −Q′ϑdΩ + Z Γs θ ςγψγ,θ+ q∗ ˆθ′dΓ + Z Γs q ςγψγ,q−ϑˆq′dΓ dt + Z T 0 Z Γs ςγ∇ψγ·n νn−ψγκνn −(∇ (qϑ) · n) νn+(qϑ) κνn dΓdt. (28) 9
All the terms that contain either θ′or q′have been grouped in Φ
1(except ˆθ′and ˆq′which cor-respond to derivatives of prescribed boundary conditions) while the remaining terms have been collected in Φ2. The initial temperature θ[x, 0] = θ0[x] is henceforth assumed to be the same for all designs s (in an Eulerian sense). Correspondingly, the spatial shape derivative of the initial condition is zero, i.e., θ′[x, 0] = 0 for all s. As a consequence of this assumption, the term ρcθ′ϑ
appearing in Φ1is zero at t = 0.
In order to eliminate the terms appearing in Φ1 that involve the implicitly-dependent shape derivatives θ′and q′, the following adjoint system is required:
ρc∂ϑ ∂τ −k∇ 2 ϑ − Q∗= 0 with Q∗= −ςωψω,θ (x, τ) ∈ Ωs× T ϑ = ˆϑ with ϑ = ςγψˆ γ,q (x, τ) ∈ Γθs× T q∗·n = − ˆq∗ with ˆq∗= −ςγψγ,θ (x, τ) ∈ Γqs× T q∗·n = −q∗e= h(ϑ − ϑe) with ϑe= − ςγψγ,θ h +ςγψγ,q (x, τ) ∈ Γ s e× T ϑ[x, 0] = 0 x ∈ Ωs. (29)
In (29), the adjoint volumetric heat supply Q∗in Ωs× T, the adjoint temperature ˆϑ in Γs
θ× T,
q es
the adjoint contact heat supply qˆ∗in Γs × T and the adjoint ambient temperature ϑein Γ × T are known once the primary field equation for θ has been solved. Observe that the characteris-tic function ς needs to be evaluated at τ = T − t. Moreover, the “initial” adjoint temperature corresponds to τ = 0 (i.e., for the final time t = T of the primary problem). The adjoint prob-lem (29) is therefore formally a transient heat conduction problem from τ = 0 to τ = T. In the numerical implementation, the boundary conditions shown in (29) are evaluated based on the solved primary system temperature field θ[x, t] at each time step. Imposing these time- and temperature-dependent boundary conditions, especially the essential boundary conditions, needs
to be careful in isogeometric analysis due to the non-interpolatory of the NURBS basis.
Numer-ical implementation aspects can be found in [50].
Choosing an adjoint temperature ϑ that satisfies (29), the term Φ1vanishes, i.e.,
Φ1= 0 , (30)
hence the shape derivative of the objective functional is provided by the term Φ2given in (28).
3.6. Continuous adjoint sensitivity
In order to obtain the final expression for the shape derivative, further simplifications for the term Φ2can be carried out. First, the spatial gradient of ψγcan be expressed as
∇ψγ=ψγ,θ∇θ + ψγ,q∇q . (31)
In view of (31) and (29)2, it follows that the terms containing gradients in the fourth integral on the right hand side of (28) can be written as
Z Γs ςγ∇ψγ·n νn−(∇ (qϑ) · n) νn dΓ = Z Γs ςγψγ,θ∇θ − q∇ϑ ·nνndΓ + Z Γs q ςγψγ,q−ϑ(∇ ˆq · n) νndΓ − Z Γs e hςγψγ,q−ϑ(∇θ · n) νndΓ. (32) 10
Substituting (32) in (28), in view of (30) and using the fact that d ˜J/ds = dJ/ds, it follows from (26) that the shape derivative of the objective functional is
dJ ds = Z T 0 Z Ωs −Q′ϑdΩ + Z Γs q ςγψγ,q−ϑ ˆ q′+(∇ ˆq · n) νn dΓ dt + Z T 0 Z Γs θ ςγψγ,θ+ q∗ ˆθ′dΓ − Z Γs e hςγψγ,q−ϑ (∇θ · n) νndΓ dt + Z T 0 ς Z Ss ω [[ω]]ψωνndΓ + Z Υs [[γ]]ψγνmdΥ ! dt + Z T 0 Z Γs ςωψω+ρc ∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ ! νndΓdt + Z T 0 Z Γs ςγψγ,θ(∇θ · n) − ψγκ −q∇ϑ · n + (qϑ) κνndΓdt. (33)
The above equation is a general expression of the continuous adjoint gradient for the opti-mization problem with an objective functional given by (3) and a temperature field that satisfies the transient heat conduction problem given by (1) and (2). This general expression can be applied for problem with design-dependent temperature, contact heat supply and/or volumetric heat supply. For simplicity, however, it is henceforth assumed that the externally-prescribed heat supplies and temperature are such that ˆθ′ = 0, ˆq′= 0, ∇ ˆq · n = 0 and Q′= 0. Under the above-mentioned assumptions, the shape derivative of the objective functional J for a transient heat conduction problem given in (33) may be expressed as
dJ ds ˆθ′, ˆq′, ∇ ˆq·n, Q′=0 = Z Γs Z T 0 gdt ! · νdΓ + Z Ss ω Z T 0 [[ω]]ψωndt ! · ν dΓ + Z Υs Z T 0 [[γ]]ψγmdt ! · ν dΥ , (34) where g = ςωψω+ρc ∂θ ∂tϑ + k∇θ · ∇ϑ − Qϑ − γeh ςγψγ,q−ϑ (∇θ · n) +ςγψγ,θ(∇θ · n) − ψγκ −q∇ϑ · n + (qϑ) κn (35)
and γein (35) is the characteristic function that is equal to 1 on the boundary Γeand 0 otherwise. For regular design points on the boundary (i.e., points where the fields are continuous), the un-constrained local shape gradient is given by the time integral of g. For singular points (including edges), the local gradient contains additional terms that provide information in the tangential direction, as indicated in (34).
4. Isogeometric analysis and design discretization
In general, the NURBS-based geometrical model can be expressed as
x =X
I
RIxI (36)
where RI is a NURBS function associated with the control point of coordinates xI and the index
I runs over all control points that characterize the geometry.
The basic idea behind isogeometric analysis is to employ the NURBS basis functions used in CAD to define the geometry also as shape functions for analysis. This approach provides a unified environment between computer added design and finite element analysis. Under this framework, the temperature field can be discretized using NURBS shape functions as follows:
θ =X
I
θIRI, I = 1, 2, · · · , a, (37)
where θI is the temperature coefficient associated to the Ith control point and a is the number of
the control points used for analysis.
The representation (37), together with the corresponding approximations for the test func-tions, the volumetric and contact heat supplies and the heat flux fields, are substituted in a weak formulation of problem (1) in order to obtain a system of equations for the unknowns θI, i.e.,
C∂θ
∂t + Kθ = f (38)
where θ = {θ1, θ2, · · · } is the vector of temperature coefficients, C is the global capacitance matrix, K is the global conductance matrix and f is the global heat supply vector.
Introducing a parameter β ∈ [0, 1] and a time step ∆t in the finite-difference approximation in time ∂(·)/∂t ≈ (1/∆t)β(·)t+∆t+ (1 − β)(·)t, the fully-discrete system of equations can be expressed from (38) as
Aθt+∆t= b (39)
with
A := C + β∆tK b := ∆t f + (C − (1 − β)∆tK) θt
and where θt+∆t corresponds to an approximation of θ at time t + ∆t, which is obtained from the approximation θt by solving (39). The cases β = 0, 0.5, 1 correspond, respectively, to the
forward Euler, Crank-Nicolson and backward Euler methods. The isogeometric approach can be used to solve the primary and adjoint problems in order to determine the fields required to compute shape derivatives.
The discretization of the analysis and design can be carried out in two different levels (see Fig. 3). Typically, in the analysis space, it is required to first carry out a refinement in order to have sufficient accuracy. This can be achieved using knot refinement. In the design space, using a coarse mesh can reduce the design parameters and simplify the process of updating the locations of the interior control points. The links between the design and analysis discretization spaces can be summarized as follows (also see [9, 24, 23]):
1. The design model is updated in the design discretization space
2. The analysis discretization space is a refined space from the design discretization space
through NURBS knot refinement
3. The sensitivity analysis is carried out in the design discretization space but using the field
variables from the analysis discretization space
x x
(a) Design discretization space (b) Analysis discretization space
Figure 3: Schematic of geometrical representation using NURBS: (a) Coarse representation to define the design geometry (design control points), (b) Refinement used for analysis (analysis
control points) that preserve the design geometry.
5. Isogeometric shape design optimization
5.1. Discrete shape gradients with respect to design control points
At the design level, the discretization typically uses a coarser mesh compared with the dis-cretization used for analysis. The design shape is expressed using (36). Following the definition in (11), the design velocity can be discretized as follows:
ν[ p; s] = ˚ˆx[ p; s] =X I RI[ p]dx I[s] ds . (40) Substituting (40) in (34) gives dJ ds ˆθ′, ˆq′, ∇ ˆq·n, Q′=0 =X I ∂J ∂xI · dxI ds . (41) where J,xI = ∂J ∂xI = Z Γs RI Z T 0 gdt ! dΓ + Z Ss ω RI Z T 0 [[ω]]ψωndt ! dΓ + Z Υs RI Z T 0 [[γ]]ψγmdt ! dΥ. (42) Similarly, for the volume constraint, it can be deduced that
Σ,xI = ∂Σ ∂xI = Z Γs RIndΓ . (43)
After solving the primary and adjoint system, all the unknown state variables involved in (42) and (43) can be evaluated. The normal vector and curvature terms can be calculated directly using relations from differential geometry.
5.2. Normalization of the search direction
In isogeometric shape optimization, the descent direction predicted from the discrete shape gradient is strongly dependent on the discretization. This discretization-dependency can slow down the convergence speed and may lead the process into a sub-optimal solution. The source of this discretization-dependency can be traced back to the lack of consistency with the local steep-est descent search direction in the continuous formulation. This inconsistency can be alleviated using a normalization approach as presented in the work of Kiendl et al. [15] and Wang et al.
0 0 0 0.005 0.01 0.015 0.02 0.025 0 0.005 0.01 0.015 0.02 0.025
Design control points Control points 1 C 2 C 3 C 4 C 5 C 6 C (a) (b) Γ1 Γ2 0.01 0.02 0.01 0.02
Figure 4: (a) Initial shape of a 2D plate with a circular orifice and (b) the NURBS parameterization for a quarter of the plate (values in m).
[51]. The normalization scheme employed in the present work is a diagonally-lumped matrix mapping-based approach [51] (also referred to as “sensitivity weighting” in [15]), which results in ¯ J,xI = J,xI WxI , (44)
where WxI is a factor that accounts for the support of NURBS function associated to the control
point xI, i.e.,
WxI =
Z
Γs
RIdΓ. (45)
5.3. Iterative descent method.
With the sensitivity analysis given in (44), the shape optimization problem presented above can be solved using different iterative algorithms (e.g., classical steepest descent or sequential quadratic programming). For simplicity, a descent algorithm that incorporates the global con-straint through a penalty-like formulation is summarized in Algorithm 1, where the design index s in the continuous formulation is formally replaced by an iteration number n and the increment
δs is formally replaced by a step size α.
6. Numerical examples 6.1. Verification examples 6.1.1. Minimum surface problem
Consider a 2D plate with a circular orifice in it as shown in Fig. 4(a). The plate has an initial temperature of θ0 = 100◦C and is placed in an environment with an ambient temperature of θe= 0◦C. The plate is made of a material with a mass density of ρ = 7800 kg/m3, a heat capacity of c = 420 J/(kg·◦C) and a thermal conductivity coefficient of k = 20 W/(m·◦C). The convection coefficient on the external boundary Γ2is h = 50 W/(m2·◦C) and on the internal boundary Γ1it is assumed that no heat is exchanged, i.e., the contact heat supply is zero ( ˆq = 0). In this case the steady state, in which ∂θ/∂t vanishes, is reached only for the uniform temperature limit case as t → ∞ that corresponds to the situation when no more heat is exchanged between the plate and
Algorithm 1 Descent algorithm with global constraint Initialize (n = 0)
Choose a initial design domain Ω(n=0)with design control points for its boundary, a step size α > 0, a penalty factor βp> 0 and a tolerance ǫ and set Λ(0)Σ = 0
Main loop (n ≥ 0)
while |J(n+1)− J(n)|/J(n)≥ǫ do
Solve the transient problem (39), compute adjoint boundary conditions and solve the adjoint problem (29)
for All control points I on the design boundary do
Compute the gradient of the objective and the volume constraint with respect to discrete variables from (42) and (43), respectively Compute the normalized search direction from (44)
end for
Constrained minimization Initialize (m = 0), Λ(n+1,0)Σ = Λ(n)Σ
while |Λ(n+1,m+1)Σ − Λ(n+1,m)Σ |/Λ(n+1,m)Σ > ǫ for Λ(n+1,m)Σ ,0 do
for All control points I on the design boundary do
Update the location of the design control points for sub-iteration m: xI(n+1,m)= xI(n)−α J(n) ,xI + Λ (n+1,m) Σ Σ (n) ,xI
and denote new location also asxI(n+1,m)
end for
Update the volume of the design region Σ(n+1,m)and Λ(n+1,m+1)Σ Λ(n+1,m+1)Σ = max0, Λ(n+1,m)Σ +βpΣ(n+1,m)− Σ0
Check convergence for case Λ(n+1,m+1)Σ = 0 separately
m ← m + 1 end while
Set Λ(n+1)Σ = Λ(n+1,m+1)Σ andxI(n+1)=
xI(n+1,m)
Update internal control points based on boundary control points n ← n + 1
end while
its surroundings. In practice, with the initial design and the given data, the plate cools down and reaches a state of near-uniform temperature of 0◦C, to within a tolerance of 10−3 ◦C, after 7200 seconds.
Suppose that one would like to find the optimal shape of the external boundary Γ2 such that the heat exchange rate with the external environment is minimized while using the same amount of material as the original design. This benchmark problem can be alternatively thought of as a minimum surface problem that admits a simple solution, namely a circular external shape (which corresponds to the minimum surface for Γ2). The verification consists of formulating the problem as a minimization of the heat exchange rate and finding the circular shape starting from the initial shape shown in Fig. 4(a). It should be noted that this problem, formulated in terms of the heat exchanged, cannot be solved based on the steady-state problem since the steady-state solution for any shape corresponds to zero temperature and zero heat exchange rate. Instead, this problem can be solved using a transient-state design approach. To achieve this, the problem can be formulated as a minimization of the total heat exchanged during a given time interval with an objective function defined over the external boundary Γ2as follows:
J = Z T 0 Ψγdt with Ψγ= Z Γ2 h(θ − θe)dΓ (46)
subjected to a volume constraint
Σ = Σ0,
where Σ0is the volume of the initial design. In this problem, Ψγ≥ 0 since the surface temperature
satisfies θ ≥ θe.
For the optimization problem, it is assumed that the model has horizontal and vertical sym-metries, so it is sufficient to consider only a quarter of the domain as shown in Fig. 4(b). Cor-respondingly, zero contact heat is assumed on the lateral boundaries x1 = 0 and x2 = 0. The locations of six control points, denoted as CI, I = 1, . . . , 6, are chosen as the discrete design
variables. The shape optimization is carried out with a design time T = 300 s. The transient solutions, for both the primary and adjoint systems, are obtained using the isogeometric analysis framework while the integration over time for the computation of the gradient is done using the trapezoidal rule. The adjoint equations corresponding to this problem can be derived from (29) as ρc∂ϑ ∂τ −k∇ 2 ϑ − Q∗= 0 with Q∗= 0 (x, τ) ∈ Ωs× T q∗·n = −q∗e= h(ϑ − ϑe) with ϑe= −1 (x, τ) ∈ Γse× T ϑ[x, 0] = 0 x ∈ Ωs. (47) 0 2
The characteristic function γ[p, t] is equal to 1 for all times if p ∈ Γ = Γ2 and zero other-wise. It is worth mentioning that the numerical sensitivity obtained from the method presented here agrees well with those calculated using finite differences, as is shown in Table 1. Due to symmetry, the control point C1is only allowed to move horizontally while C6is only allowed to
move vertically. The corresponding sensitivity analysis regarding to the restricted directions are
omitted.
Using the iterative optimization approach summarized in Algorithm 1, the optimization pro-cess converges within 5 iterations to within a small tolerance of 10−4as shown in Fig. 5(a). The corresponding optimal shape is shown in Fig. 5(b) in comparison to the initial and the exact
0 0.005 0.01 0.015 0.02 0.025 0 0.005 0.01 0.015 0.02 0.025 0 5 10 15 20 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5 5.1 5.2× 10 4 Objective funcional [J] Iteration step x1 [m] Optimal shape Initial shape Circular shape (a) (b) x2 [m]
Figure 5: (a) Iteration history and (b) the optimal shape.
Table 1: Comparison between adjoint gradient and finite difference gradient for the transient heat dissipation problem.
CI Component Adjoint gradient Finite difference gradient Relative difference
C1 1 1.3528×105 1.3665×105 1.01% 2 ∼ ∼ ∼ C2 1 2.4255×105 2.4427×105 0.71% 2 -0.3706×105 -0.3733×105 0.75% C3 1 8.9649×105 8.9665×105 0.018% 2 3.0640×105 3.0542×105 0.32% C4 1 3.0640×105 3.0542×105 0.32% 2 8.9649×105 8.9665×105 0.018% C5 1 -0.3706×105 -0.3733×105 0.75% 2 2.4255×105 2.4427×105 0.71% C6 1 ∼ ∼ ∼ 2 1.3528×105 1.3665×105 1.01%
circular shape. As can be seen from Fig. 5(b), the obtained optimal shape agrees closely with the circular shape, which provides a validation of the design methodology. The control points and
the corresponding weights of the initial and optimal designs are presented in Tabel 2. The knot
vectors used for the two directions in the index space are ξ = [0 0 0 1/3 1/2 2/3 1 1 1] and η =
[0 0 0 1 1 1], respectively.
6.1.2. Passive temperature control problem
In this example, a heat flux q = 10kW/m2 applied to the left side (Γ1) of a 2D plate, as shown in Fig. 6. The rest of the boundary is isolated such that the contact heat flux on Γ2, 3
can be neglected. The material properties of the plate are assumed to be the same as in the
previous example. The initial temperature of the plate is assumed to be 0◦C. With the heat flux
on the left side, the temperature on the right side of the plate (Γ2) reaches a (nearly) uniform temperature of about 36◦C at t = 400 s. By modifying the shape of the upper boundary (Γ3) of the plate, it is possible to control the temperature on the right side (Γ2) such that it reaches a
Table 2: Geometry information of the heat dissipation problem
Initial design Optimal design
I(i, j) Locations Weights I(i, j) Locations Weights
xI 1 x I 2 x I 1 x I 2 (1, 1) 0.0100 0.0000 1.00 (1, 1) 0.0100 0.0000 1.00 (2, 1) 0.0100 0.0026 0.90 (2, 1) 0.0100 0.0026 0.90 (3, 1) 0.0080 0.0061 0.85 (3, 1) 0.0080 0.0061 0.85 (4, 1) 0.0061 0.0080 0.85 (4, 1) 0.0061 0.0080 0.85 (5, 1) 0.0026 0.0100 0.90 (5, 1) 0.0026 0.0100 0.90 (6, 1) 0.0000 0.0100 1.00 (6, 1) 0.0000 0.0100 1.00 (1, 2) 0.0150 0.0000 1.00 (1, 2) 0.0150 0.0000 1.00 (2, 2) 0.0150 0.0039 0.90 (2, 2) 0.0150 0.0039 0.90 (3, 2) 0.0121 0.0091 0.85 (3, 2) 0.0121 0.0091 0.85 (4, 2) 0.0091 0.0121 0.85 (4, 2) 0.0091 0.0121 0.85 (5, 2) 0.0039 0.0150 0.90 (5, 2) 0.0039 0.0150 0.90 (6, 2) 0.0000 0.0150 1.00 (6, 2) 0.0000 0.0150 1.00 (1, 3) 0.0200 0.0000 1.00 (1, 3) 0.0235 0.0000 1.00 (2, 3) 0.0200 0.0100 1.00 (2, 3) 0.0230 0.0092 1.00 (3, 3) 0.0238 0.0213 1.00 (3, 3) 0.0172 0.0163 1.00 (4, 3) 0.0213 0.0238 1.00 (4, 3) 0.0163 0.0172 1.00 (5, 3) 0.0100 0.0200 1.00 (5, 3) 0.0092 0.0230 1.00 (6, 3) 0.0000 0.0200 1.00 (6, 3) 0.0000 0.0235 1.00
target distribution at a desired time. In particular, suppose that the target at t = T = 400 s is a uniformly distributed temperature of 40◦C on the right side. The objective of this problem can be formulated as follows: J = Z T 0 ςΨγdt with Ψγ= Z Γ2 (θ − ˇθ)2dΓ (48)
where ς = δ[t − t0] is the characteristic function for the time interval, δ[t − t0] is the Dirac delta function with t0 = T = 400 s and θˇ = 40◦C is the target temperature. The adjoint equations corresponding to this problem can be derived from (29) as
ρc∂ϑ ∂τ −k∇ 2 ϑ − Q∗= 0 with Q∗0 (x, τ) ∈ Ωs× T q∗·n = − ˆq∗ with ˆq∗= −2ς(θ − ˇθ) (x, τ) ∈ Γsq× T ϑ[x, 0] = 0 x ∈ Ωs. (49)
The numerical sensitivity obtained from the method presented here agrees well with those
calcu-lated using finite differences, as is shown in Table 3. The optimization process converges
within 6 iterations to within a small tolerance of 10−4as shown in Fig. 7(a). The corresponding optimal shape and the temperature contours at t = 400 s are plotted in Fig. 7(b), from which it can be seen that the temperature at the right side matches the target temperature very well. To
Table 3: Comparison between adjoint gradient and finite difference gradient for temperature control problem.
I(i, j) Component Adjoint gradient Finite difference gradient Relative difference
(3,2) 1 13.9820 14.3893 2.91% 2 55.9280 57.5575 2.91% (3,3) 1 13.9571 14.3640 2.91% 2 55.8286 57.4564 2.91% 0 0.01 0.02 0.03 0.04 0 0.01 0.02 Insulating material Design control points Control points (a) (b) Γ3 0 0.01 0.02 0.03 0.04 0 0.01 0.02 36 38 40 42 44 46 48 50 52 t = 400 s oC Γ2 Γ4 Γ1
Figure 6: (a)The plate immersed in the isolation material and (b) the temperature contour at t = 400s.
a nearly uniform temperature distribution, the shape optimization procedure effectively creates a uni-directional heat flow pattern on the right side of the domain while the target temperature can be reached using less material compared to the initial design (i.e., in this case the resource con-straint is not active). The control points and the corresponding weights of the initial and optimal designs are presented in Table 4. The knot vectors used for the two directions in the index space
are ξ = [0 0 0 1 1 1] and η = [0 0 0 1/2 1 1 1], respectively.
6.2. Shape optimization of a plunger
The shape optimization methodology is tested in this section with an example drawn from
the literature, namely a two-dimensional plunger that is designed to form a television bulb panel
(a) (b) 0 0.01 0.02 0.03 0.04 0 0.01 0.02 36 38 40 42 44 46 48 50 52 t = 400 s 1 2 3 4 5 6 7 8 9 10 0 0 0.04 0.08 0.12 0.16 Objective functional [ oC 2 m 2] oC
Figure 7: (a) Iteration history and (b) the optimal shape and the temperature contour at t = 400s. 19
Table 4: Geometry information of the temperature control problem
Initial design Optimal design
I(i, j) Locations Weights I(i, j) Locations Weights
xI 1 x I 2 x I 1 x I 2 (1, 1) 0.0000 0.0000 1.00 (1, 1) 0.0000 0.0000 1.00 (2, 1) 0.0000 0.0100 1.00 (2, 1) 0.0000 0.0100 1.00 (3, 1) 0.0000 0.0200 1.00 (3, 1) 0.0000 0.0200 1.00 (1, 2) 0.0100 0.0000 1.00 (1, 2) 0.0100 0.0000 1.00 (2, 2) 0.0100 0.0088 1.00 (2, 2) 0.0100 0.0088 1.00 (3, 2) 0.0100 0.0175 1.00 (3, 2) 0.0100 0.0174 1.00 (1, 3) 0.0300 0.0000 1.00 (1, 3) 0.0300 0.0000 1.00 (2, 3) 0.0300 0.0063 1.00 (2, 3) 0.0300 0.0063 1.00 (3, 3) 0.0300 0.0125 1.00 (3, 3) 0.0289 0.0081 1.00 (1, 4) 0.0400 0.0000 1.00 (1, 4) 0.0400 0.0000 1.00 (2, 4) 0.0400 0.0050 1.00 (2, 4) 0.0400 0.0050 1.00 (3, 4) 0.0400 0.0100 1.00 (3, 4) 0.0400 0.0100 1.00 0 40 80 120 160 200 0 20 40 60 80 100 Control points Design control points
Γ1 Γ2 Γ3 Γ4 1 C 2 C 3 C C4 C5 θe1 θe3 Symmetric
Figure 8: A plunger model and its NURBS control points (distances in mm).
from molten glass [52], as is shown in Fig. 8. Although this example pertains to a technology that is not currently widely used, it is an existing reference that serves to illustrate the main features of the method and may be applicable to similar situations. As shown in the figure, the workpiece is brought into contact with a surface at a temperature of θe3 = 1000◦C on the boundary Γ3. The boundary Γ1 is a cooling surface that exchanges heat with a cooling fluid at a temperature of θe1 = 0◦C. The quality of the product surface is affected by temperature fluctuations on the contact surface Γ3 [52]. An optimization is formulated where the shape of the surface Γ1 is modified in order to reduce the temperature variance on the contact surface Γ3. To this end, denote as ˇθ = ˇθ[t] the average temperature along Γ3at time t, i.e.,
ˇθ = 1 |Γ3| Z Γ3 θdΓ with |Γ3| = Z Γ3 dΓ 20
120 140 160 180 200 0 20 40 60 80 100 1 C 100 80 60 40 20 0 3 C 2 C 4 C 5 C 120 140 160 180 200 0 20 40 60 80 100 100 80 60 40 20 0
Figure 9: The irregular geometry and mesh during the iteration without restricting the movements of the control points.
and define a performance functional as J = Z T 0 Ψγdt with Ψγ= Z Γ3 (θ − ˇθ)2dΓ, s.t. x1I ≥20 mm, xI2≥20 mm I = 1, . . . , 5 (50)
where xiI, i = 1, 2 are the coordinates of the control points CI, I = 1, . . . , 5 as shown in the
figure. The locations of the control points CI that characterize the cooling boundary Γ1are used
as design variables. For simplicity, the design control points C1 and C2 are only allowed to move horizontally, C4and C5are only allowed to move vertically, and C3is allowed to move both horizontally and vertically. The restrictions of these movements is to avoid the mesh and geometry irregularity. Without these restrictions, the geometry and mesh may become irregular,
as is shown in Fig. 9. The convection coefficients are h = 3.15 × 10−4W/(mm2·◦C) on Γ1and h
= 2.88 × 10−4W/(mm2·◦C) on Γ3, respectively. The thermal conductivity coefficient is k =
27.52 × 10−3W/(mm ·◦C) and the effective heat capacity is (ρc) = 2.288 × 10−3J/(mm3·◦C). On the boundaries Γ2and Γ4it is assumed that there is no heat exchange and the initial temperature on the whole domain is taken as 0◦C).
In order to study the influence of the design time T, the optimization problem was solved
with two distinct values, namely T = 500 s and T = 1000 s. The sensitivity obtained using the continuous adjoint method was verified with the finite difference method, which are shown in Table 5 and Table 6, respectively. The optimal shapes obtained are plotted in Fig. 10 for the two design times. As may be seen in the figure, the design for T = 500 s has a more complex shape than the design for T = 1000 s, which may be ascribed to the influence of the details of the transient states. In contrast, the design for T = 1000 s is closer to the steady state limit case, which is more heavily influenced by stationary conditions. The corresponding convergence histories are plotted in Fig. 11(a) for T = 500 s and Fig. 11(b) for T = 1000 s. These convergence histories also reflect the complexity of the designs in the sense that generally more iterations are required to converge for designs that are more influenced by the transient states (in this case more iterations are required to converge for the design for T = 500 s than for T = 1000
s). The control points and the corresponding weights of the initial and optimal designs are
presented in Table 7. The knot vectors used for the two directions in the index space are ξ = [0 0
0 0.2 0.4 0.5 0.7 1 1 1] and η = [0 0 0 1 1 1], respectively.
The space-averaged temperature fluctuations on the boundary Γ3 of the initial and the two optimal design shapes during the first 1000 seconds are plotted in Fig. 12(a) in terms of the
0 20 40 60 80 100 120 140 160 180 200 0 10 20 30 40 50 60 70 80 90 100 T = 1000 s T = 500 s
Figure 10: Optimal shapes for design times T = 500s and T = 1000s (distances in mm)
O je ct ive func ti ona l [ oC 2. m m 2. s ] 0 5 10 15 20 25 4 4.5 5 5.5 6 6.5 7 7.5 × 108 0 5 10 15 20 25 30 2. 2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8× 108 (a) (b)
Iteration step Iteration step
O bj ec ti ve func ti ona l [ oC 2. m m 2. s ] T = 1000 s T = 500 s
Figure 11: Iteration history of the optimization for design times (a) T = 500 s and (b) T = 1000 s.
corresponding transient temperature variance, which is evaluated as ˜θ[t] = qΨγ[t]/|Γ3|,
where, as before, |Γ3|denotes the length of the boundary Γ3. The quantity ˜θ[t] measures, at each time t, the average deviation of the local temperature from the average temperature (a larger value indicating an undesirable fluctuation).
From the plots in Fig. 12(a), it can be seen that the optimal shapes for both T = 500 s and T = 1000 s have a lower variance than the initial shape throughout the time interval 0 ≤ t ≤ 1000 s except during the first 150 s for the case T = 1000 s, in which the variance is slightly larger. However, the time-averaged variance for both designs is lower than the original design. The optimal shape for the case T = 500 s provides the lowest temperature variance during the first 410 seconds, but after that the temperature variance of the case T = 1000 s becomes the lowest. This result is consistent with the expected influence of the design parameter T on the design, namely that a design performs better during its own design interval [0, T ] but not necessarily during other times intervals (either larger or shorter than the design interval).
In terms of local temperature fluctuations, as shown in Fig. 12(b), the optimization proce-dure reduces the ratio of the maximum temperature θmaxon Γ3to the average temperature ˇθ and, simultaneously, increases the ratio of the minimum temperature θminon Γ3to the average temper-ature ˇθ, as compared with the corresponding ratios for the original design. The evolution of the average temperature ˇθ on Γ3 is shown in Fig. 12(c), which indicates that both optimal designs generally decrease the average temperature compared with the average temperature of the initial
0 100 200 300 400 500 600 700 800 900 1000 0 10 20 30 40 50 60 0 100 200 300 400 500 600 700 800 900 1000 0.9 1 1.1 1.2 1.3 1.4 0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 (a) (b) (c) Initial shape Optimal shape (T=500 s) Optimal shape (T=1000 s) Initial shape Optimal shape (T=500 s) Optimal shape (T=1000 s) Initial shape Optimal shape (T=500 s) Optimal shape (T=1000 s) θ (t) θ(t) max θ (t) θ(t) min Time t [s] Time t [s] Time t [s] T em pe ra ture va ri anc e [ C] T em pe ra ture ra ti o A ve ra ge t em pe ra ture [ C ]
Figure 12: (a) The temperature variance ˜θ, (b) the ratios of the maximum and minimum temperatures to the average temperature and (c) the transient average temperature on the boundary Γ3for 0 ≤ t ≤ 1000 s for the initial and optimal (T = 500 s and T = 1000 s) designs.
design.
The example presented in this section illustrates the significant effect that the transient re-sponse has on the optimal design and highlights the importance of making a judicious choice of the design parameters depending on the desired performance of the design.
6.3. Shape optimization of a thermal protection panel
In this example, the structural shape and the heating rate of a ballistic re-entry vehicle are considered, as shown in Fig. 13(a) and Fig. 14, which is adapted from [53]. The vehicle has two parts: the recovery tip (forward part) and the destructible frustum (aft part, see Fig. 13(a)). The outer boundary Γ1has a heavy heating rate during re-entry, as shown in Fig. 14 for 0 ≤ t ≤ 50s, which is used as (time-dependent) boundary condition. For times greater than 50 s, it is assumed that the heating rate is negligible in comparison to the rate for 0 ≤ t ≤ 50s. In order to protect the internal equipments inside the tip, a layer of low thermal conductivity material is used to insulate the heat from outside the tip. The cross section of the recovery tip and its NURBS parameterization are shown in Fig. 13(b). The thermal conductivity coefficient of the thermal insulation material is k = 8.0 × 10−5W/(mm·◦C), the heat capacity is c = 514J/(kg·◦C) and the material density is ρ = 0.22 × 10−6kg/mm3. The initial temperature is assumed to be 0◦C, the convection coefficient on Γ2is taken h = 6.0 × 10−3W/(mm2· ◦C), the ambient temperature on Γ2is assumed to be 0◦C and on the remaining boundaries the heat exchanged is taken as zero.
With the goal of minimizing the heat transferred through the boundary Γ2, an objective func-tional is formulated as J = Z T 0 Z Γ2 h(θ − θe)dΓdt (51) 23
R=158.75
262.5 378.9
14
θ
Recovery tip Destructible frustum
L −3000 −250 −200 −150 −100 −50 0 20 40 60 80 100 120 140 160 180 200 Control points Design control points
1 C 2 C 3 C 4 C 5 C Γ2 Γ1 (a) (b)
Figure 13: (a) The sketch of a re-entry vehicle (adapted from [53]) and (b) the cross section of the recovery tip skin and its NURBS parameterization (distances in mm).
0 10 20 30 40 50 0 60120 180240 300360 0 0.5 1 1.5 2 Bounda ry lengt h [mm] Time [s] H ea ti ng ra te [J /(m m 2. s )]
Figure 14: The time-dependent heating rate along the boundary Γ1(data from [53]).
Table 5: Comparison between adjoint gradient and finite difference gradient for the plunger design with T = 500s.
I Component Adjoint gradient Finite difference gradient Relative difference C1 1 -6.0855×105 -6.1078×105 0.36% 2 0 2.6388 ∼ C2 1 -3.1106×105 -3.0308×105 2.57% 2 3.3747×105 3.4454×105 2.10% C3 1 2.0409×106 2.0891×106 2.36% 2 1.1000×106 1.1451×106 4.10% C4 1 0.8443×105 0.8658×105 2.54% 2 -6.2055×105 -6.2875×105 1.32% C5 1 0 6.3077 ∼ 2 -1.5463×106 -1.5391×106 0.47%
subjected to a volume constraint
Σ ≤ Σ0,
where Σ0 is the volume of the initial design. The design time is chosen as T = 180 s, which includes an initial time interval with a large heating rate followed by a time interval with a lower rate in order to take both loading scenarios into account. The locations of five design control points, CI, I = 1, . . . , 5, which characterize the internal boundary Γ2, are chosen as design
variables. Control point C1is only allowed to move horizontally while control point C5 is only allowed to move vertically. The external boundary Γ1is kept fixed since its optimization depends mostly on aerodynamic requirements. Further, it is assumed that the heating rate on the external boundary, as shown in Fig. 14, is the same regardless of the shape of the internal boundary. As in previous examples, the sensitivity obtained using the continuous adjoint method was verified with the finite difference method, which is shown in Table 8. It should be noted that for the sensitivity of the design control point C5, extra terms need to be included into the sensitivity formulation due to the geometry discontinuity involved. The optimization is performed using the Algorithm 1 and the convergence history is shown Fig. 15(a), which indicates that after about six iterations the process has converged. The optimal shape obtained is plotted in Fig. 15(b). From the figure, it can be seen that the optimal shape requires a thicker layer close to the tip, which comes at the expense of the thickness in the rear part in order to preserve the total volume. With this change in shape, the heat conducted inside the front section can be reduced about 11.5% compared with the original design, as shown in Fig. 15(a). The control points and the corresponding weights of the initial and optimal designs are presented in 9. The knot vectors
used for the two directions in the index space are ξ = [0 0 0 1/2 1/2 1 1 1] and η = [0 0 0 1 1 1],
respectively.
7. Conclusions
In this work, the continuous adjoint shape sensitivity analysis for transient heat conduction problems is reformulated taking into consideration the discontinuities involved in the objective
Table 6: Comparison between adjoint gradient and finite difference gradient for the plunger design with T = 1000s.
I Component Adjoint gradient Finite difference gradient Relative difference C1 1 -6.6243×105 -6.6298×105 0.08% 2 0 2.6388 ∼ C2 1 -0.6969×106 -7.1838×106 3.09% 2 1.0752×106 1.0903×106 1.40% C3 1 7.0695×106 7.1734×106 1.47% 2 6.3920×106 6.5046×106 1.76% C4 1 0.3165×106 0.3213×106 1.52% 2 -1.2643×106 -1.2761×106 0.93% C5 1 0 3.7868 ∼ 2 -5.4290×105 -5.4290×105 2.94% (a) (b) Iteration step O bj ec ti ve func ti ona l [J ] −3000 −250 −200 −150 −100 −50 0 20 40 60 80 100 120 140 160 180 200 Initial shape Optimal shape 0 2 4 6 8 10 12 540 550 560 570 580 590 600 610
Figure 15: (a) Iteration history of the optimization and (b) the optimal shape.
functionals. The continuous sensitivity analysis, which is applicable to general shape optimiza-tion problems, is discretized in the context of isogeometric analysis. The control points of a NURBS description of the shape are used as design variables, which allows a seamless integra-tion between optimizaintegra-tion and analysis. The optimizaintegra-tion and the analysis are performed at two levels of discretization (coarse and fine, respectively), but no loss of geometrical information occurs in this process. The methodology was tested with benchmark problems, which also illus-trate the flexibility provided by the characteristic functions to measure the design performance in selected places and times. The transient isogeometric shape optimization was subsequently ap-plied to cases where thermal conditions fluctuate during operation such as the thermal protection system design for a reentry ballistic vehicle. In these examples, it is shown that shape optimiza-tion with accountability of transient states is an attractive approach in applicaoptimiza-tions where active control is not economically feasible. Furthermore, it is possible to combine shape optimization (as passive control) with an active control approach in view of increasing the efficiency of the thermal management. Following a similar approach, it is also possible to further extend the methodology and framework to include other transient situations, such as mechanical problems
