An isogeometric numerical study of partially and fully implicit schemes for transient adjoint shape sensitivity analysis

Abstract

In structural design optimization involving transient responses, time integration scheme plays a crucial role in sensitivity analysis because it affects the accuracy and stability of transient analysis. In this work, the influence of time integration scheme is studied numerically for the adjoint shape sensitivity analysis of two benchmark transient heat conduction problems within the framework of isogeometric analysis. It is found that (i) the explicit approach (\beta = 0) and semi-implicit approach with \beta<0.5 impose a strict stability condition of the transient analysis; (ii) the implicit approach (\beta=1) and semi-implicit approach with \beta > 0.5 are generally preferred for their unconditional stability; and (iii) Crank–Nicolson type approach (\beta=0.5) may induce a large error for large time-step sizes due to the oscillatory solutions. The numerical results also show that the time-step size does not have to be chosen to satisfy the critical conditions for all of the eigen-frequencies. It is recommended to use \beta \approx \mathrm{0.75} for unconditional stability, such that the oscillation condition is much less critical than the Crank–Nicolson scheme, and the accuracy is higher than a fully implicit approach.

Key words： isogeometric shape optimization; design-dependent boundary condition; transient heat conduction; implicit time integration; adjoint method

Introduction

Designing structures under thermal loadings is an important issue in practical engineering problems [¹,²]. Numerical design optimizations for such problems often involve design sensitivity analysis, which is demonstrated in numerous articles, mostly for steady state heat conduction problems, such as in Refs. [³–⁷] with topology optimizations and in Refs. [⁸–¹⁰] with shape optimizations.

Structural design sensitivity analyses involving transient thermal responses have been studied since the late 1970s with the work of Ref. [¹¹], where approximation concepts were adopted and only critical time points were used. Following this, some fundamental technical issues for thermal sensitivity analyses were discussed in Ref. [¹²], which include the computational efficiency of explicit and implicit approaches. It was demonstrated that the implicit approaches are computationally more efficient than the explicit approach, because the critical stability condition required in the latter results in relatively small time-steps. There are a few works in literature, e.g., Refs. [¹³,¹⁴], adopting an explicit time integration scheme, though mostly utilized for simple heat conduction problems. Most design sensitivity analyses involving transient responses adopt either the semi- or fully implicit approach, e.g., in Refs. [¹⁵,¹⁶] with a Crank–Nicolson type semi-implicit approach, in Refs. [¹⁷,¹⁸] with general implicit approaches, and in Refs. [¹⁹,²⁰] with a fully implicit approach. In Refs. [⁹,¹⁰,²¹–²⁶], variational continuum adjoint methods for shape sensitivity analysis of transient heat conduction problems are presented without elaborating on the time integration schemes. There are also cases where the problems associated with typical time integration approaches may not appear, e.g., in Ref. [²⁷] where fully analytical solution can be obtained, or in Refs. [²⁸,²⁹] where precise time integration scheme is used.

In general, analytical structural shape sensitivity analysis can be complicated (see Refs. [⁸,³⁰,³¹]). To this end, isogeometric analysis (IGA) can significantly reduce the difficulties of performing shape sensitivity analysis, due to its ability to preserve exact geometrical features and high-order continuities [³²,³³]. These features are attractive for the development of an integrated design frame work for curved beams (e.g., in Refs. [³⁴–³⁷]), shells (e.g., in Refs. [³⁸,³⁹]) and general curved structures (e.g., in Refs. [⁴⁰–⁴⁵]). The ease of achieving multiple resolutions, and the high order shape functions of IGA, also promote the development of topology optimization, e.g., in Refs. [⁴⁶–⁵²]. A generalized shape optimization method combining level set method and finite cell method [⁵³] for structural designs using IGA can be found in Refs. [⁵⁴,⁵⁵]. More information of IGA-based design optimization can be found in a recent review in Ref. [⁵⁶].

The IGA reduces the numerical error induced by spatial discretizations in a shape sensitivity analysis. However, for time-dependent problems, the time integration scheme has a critical influence on the accuracy of sensitivity. Such time-dependent problems include the design of lattice structures incorporating heat conduction considerations [⁵⁷,⁵⁸], to give interesting thermal behaviors [⁵⁹]. In this work, the partially and fully implicit time integration schemes for shape sensitivity analysis of transient heat conduction problems are studied numerically within the framework of IGA. The findings provide a reference on the appropriate time integration schemes for problems with transient responses, which also include design optimization problems in the nonlinear deformation, e.g., Ref. [⁶⁰].

Problem statement and adjoint shape sensitivity analysis

The problem considered is the design of a given structure made from an isotropic material with linear thermal conduction properties, with domain {\mathrm{\Omega }}^s and boundary {\mathrm{\Gamma }}^s as shown in Fig. 1. The variable s denotes the design stages that are sequentially modified from the referential/initial design {\mathrm{\Omega }}^{\mathrm{0}}. The location function of a material point \mathit{p} in {\mathrm{\Omega }}^s is denoted as \mathit{x}[\mathit{p},\mathrm{\hspace{0.17em}}s]. The temperature at location \mathit{x}, time t and design stage s is denoted by \theta [\mathit{x},\mathrm{\hspace{0.17em}}t\mathrm{;}\mathrm{\hspace{0.17em}}s]. A general design objective functional defined over a time interval \mathcal{T}=[\mathrm{0},\mathrm{\hspace{0.17em}}T] can be characterized as

(1)\mathcal{J}:={\displaystyle {\int }_{\mathrm{0}}^T}\varsigma \left[t\right]({\displaystyle {\int }_{{\mathrm{\Omega }}^s}{\psi }_{\omega }\left[\theta \right[\mathit{x},t\mathrm{;}s\left]\right]}\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}^s}{\psi }_{\mathit{\gamma }}\left[\theta \right[\mathit{x},t\mathrm{;}s],q[\mathit{x},t\mathrm{;}s\left]\right]\mathrm{d}\mathrm{\Gamma }})\mathrm{d}t,

in which the time characteristic function \varsigma \left[t\right] is defined as

(2)\varsigma \left[t\right]=\{\begin{array}{ll}\mathrm{1},\hfill & \mathrm{if}\mathrm{\hspace{0.17em}}t\in {\mathcal{T}}_{\varsigma }\subseteq \mathcal{T},\hfill \\ \mathrm{0},\hfill & \mathrm{otherwise,}\hfill \end{array}

and q[\mathit{x},t\mathrm{;}s] is the heat flux on the boundary.

Fig.1 Schematics of initial design at s=\mathrm{0} (left) and updated design at s (right).

Full size|PPT slide

The transient heat conduction problem within a time interval \mathcal{T}=[\mathrm{0},T] is governed by

(3)\{\begin{array}{ll}l\left[\theta [\mathit{x},t]\right]:=\rho c{\displaystyle \frac{\partial \theta [\mathit{x},t]}{\partial t}}-k{\nabla }^{\mathrm{2}}\theta [\mathit{x},t]-Q[\mathit{x},t]=\mathrm{0}\mathrm{\hspace{1em}}\hfill & (\mathit{x},t)\in {\mathrm{\Omega }}^s\times \mathcal{T},\hfill \\ \theta [\mathit{x},t]=\widehat{\theta }[\mathit{x},t]\mathrm{\hspace{1em}}\hfill & (\mathit{x},t)\in {\mathrm{\Gamma }}_{\theta }^s\times \mathcal{T},\hfill \\ \mathit{q}[\mathit{x},t]\cdot \mathit{n}\left[\mathit{x}\right]=-k\nabla \theta [\mathit{x},t]\cdot \mathit{n}\left[\mathit{x}\right]=-\widehat{q}[\mathit{x},t]\mathrm{\hspace{1em}}\hfill & (\mathit{x},t)\in {\mathrm{\Gamma }}_q^s\times \mathcal{T},\hfill \\ \mathit{q}[\mathit{x},t]\cdot \mathit{n}\left[\mathit{x}\right]=-k\nabla \theta [\mathit{x},t]\cdot \mathit{n}\left[\mathit{x}\right]=-q_e[\mathit{x},t]=h\left(\theta \right[\mathit{x},t]-{\theta }_e[t\left]\right)\mathrm{\hspace{1em}}\hfill & (\mathit{x},t)\in {\mathrm{\Gamma }}_e^s\times \mathcal{T},\hfill \\ \theta [\mathit{x},\mathrm{0}]={\theta }_{\mathrm{0}}\left[\mathit{x}\right]\mathrm{\hspace{1em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{1em}}\hfill & (\mathit{x},\mathrm{0})\in {\mathrm{\Omega }}^s,\hfill \end{array}

where c>\mathrm{0}, \rho >\mathrm{0}, and k>\mathrm{0} are the heat capacity, mass density, and thermal conductivity, respectively, Q denotes the body heat generation rate of each volume unit, \mathit{q}=-k\nabla \theta and q=-\mathit{q}\cdot \mathit{n} are the heat flow inside {\mathrm{\Omega }}^s and heat flux through {\mathrm{\Gamma }}^s, respectively, \mathit{n} represents the unit outward normal vector on the boundary, h denotes the heat convection coefficient at ambient environment, {\theta }_e is the ambient temperature, {\theta }_{\mathrm{0}}\left[\mathit{x}\right] represents the initial temperature field in the domain, and \nabla and {\nabla }^{\mathrm{2}} are the gradient and the Laplacian operators, respectively. The problems considered in this work have an essential boundary condition \widehat{\theta } defined over {\mathrm{\Gamma }}_{\theta }^s, a natural boundary condition \widehat{q} defined over {\mathrm{\Gamma }}_q^s and a Robin boundary condition defined over {\mathrm{\Gamma }}_e^s.

The governing equation in Eq. (3) can be treated as equality constraint of the optimization problem with design objective formulated in Eq. (1). To satisfy these equality constraints imposed over the entire domain and time, we introduce an augmented functional

(4)\overline{\mathfrak{J}}\left[s\right]=\mathfrak{J}\left[s\right]+{\displaystyle {\int }_{\mathrm{0}}^T{\langle l\left[\theta \right],\vartheta \rangle }_{{\mathrm{\Omega }}^s}\mathrm{d}t},

in which \vartheta =\vartheta [\mathit{x},t] is the adjoint temperature field with the equality constraint nested in

(5){\langle l\left[\theta \right],\vartheta \rangle }_{{\mathrm{\Omega }}^s}={\displaystyle {\int }_{{\mathrm{\Omega }}^s}\left(\rho c{\displaystyle \frac{\partial \theta }{\partial t}}\vartheta +k\nabla \theta \cdot \nabla \vartheta -Q\vartheta \right)}\mathrm{d}\mathrm{\Omega }-{\displaystyle {\int }_{{\mathrm{\Gamma }}^s}q\vartheta \mathrm{d}\mathrm{\Gamma }=\mathrm{0.}}

With the constraint satisfied, we have \mathfrak{J}=\overline{\mathfrak{J}}. The transient adjoint temperature field can be solved with the following adjoint problem

(6)\{\begin{array}{ll}\rho c{\displaystyle \frac{\partial \vartheta }{\partial \tau }}-k{\nabla }^{\mathrm{2}}\vartheta -Q^{*}=\mathrm{0},Q^{*}=-\varsigma {\psi }_{\omega ,\theta }\hfill & (\mathit{x},\tau )\in {\mathrm{\Omega }}^s\times \mathcal{T},\hfill \\ \vartheta =\widehat{\vartheta },\widehat{\vartheta }=\varsigma {\psi }_{\mathit{\gamma },q}\hfill & (\mathit{x},\tau )\in {\mathrm{\Gamma }}_{\theta }^s\times \mathcal{T},\hfill \\ {\mathit{q}}^{*}\cdot \mathit{n}=-\hat{q}{}^{*},\hat{q}{}^{*}=-\varsigma {\psi }_{\mathit{\gamma },\theta }\hfill & (\mathit{x},\tau )\in {\mathrm{\Gamma }}_q^s\times \mathcal{T},\hfill \\ {\mathit{q}}^{*}\cdot \mathit{n}=-q_e^{*}=h\left(\vartheta -{\vartheta }_e\right),\begin{array}{l}\hfill \end{array}{\vartheta }_e=-{\displaystyle \frac{\varsigma {\psi }_{\mathit{\gamma },\theta }}{h}}+\varsigma {\psi }_{\mathit{\gamma },q}\hfill & (\mathit{x},\tau )\in {\mathrm{\Gamma }}_e^s\times \mathcal{T},\hfill \\ \vartheta [\mathit{x},\mathrm{0}]=\mathrm{0}\mathrm{\hspace{1em}\hspace{1em}}\hfill & (\mathit{x},\mathrm{0}){|}_{t=T}\in {\mathrm{\Omega }}^s\mathrm{\hspace{0.17em}},\hfill \end{array}

in which \tau =T-t is adjoint time, Q^{*} is adjoint volumetric heat supply, \widehat{\vartheta } and {\widehat{q}}^{*} are adjoint Dirichlet and Neumann boundary conditions, respectively, and {\vartheta }_e is the ambient temperature. Eventually, the shape sensitivity with respect to parameter s can be derived as

(7)$$\begin{gathered} {\displaystyle \frac{\mathrm{d}\mathcal{J}}{\mathrm{d}s}}={\displaystyle {\int }_{\mathrm{0}}^T\left({\displaystyle {\int }_{{\mathrm{\Omega }}^s}-Q^{\prime }\vartheta \mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_q^s}\left(\varsigma {\psi }_{\mathit{\gamma },q}-\vartheta \right)\left({\widehat{q}}^{\prime }+\left(\nabla \widehat{q}\cdot \mathit{n}\right){\mathit{\nu }}_n\right)\mathrm{d}\mathrm{\Gamma }}}\right)}\mathrm{d}t \\ +{\displaystyle {\int }_{\mathrm{0}}^T}\left({\displaystyle {\int }_{{\mathrm{\Gamma }}_{\theta }^s}}\left(\varsigma {\psi }_{\mathit{\gamma },\theta }-q*\right){\widehat{\theta }}^{\prime }\mathrm{d}\mathrm{\Gamma }-{\displaystyle {\int }_{{\mathrm{\Gamma }}_e^s}}h\left(\varsigma {\psi }_{\mathit{\gamma },q}-\vartheta \right)\left(\nabla \theta \cdot \mathit{n}\right){\mathit{\nu }}_n\mathrm{d}\mathrm{\Gamma }\right)\mathrm{d}t\mathrm{d}t \\ +{\displaystyle {\int }_{\mathrm{0}}^T}{\displaystyle {\int }_{{\mathrm{\Gamma }}^s}}\left(\varsigma {\psi }_{\omega }-\rho c{\displaystyle \frac{\partial \theta }{\partial t}}\vartheta +k\nabla \theta \cdot \nabla \vartheta -Q\vartheta \right){\mathit{\nu }}_n\mathrm{d}\mathrm{\Gamma }\mathrm{d}t \\ +{\displaystyle {\int }_{\mathrm{0}}^T}{\displaystyle {\int }_{{\mathrm{\Gamma }}^s}\left(\varsigma \left({\psi }_{\mathit{\gamma },\theta }\left(\nabla \theta \cdot \mathit{n}\right)-{\psi }_{\mathit{\gamma }}\kappa \right)-q\nabla \vartheta \cdot \mathit{n}+\left(q\vartheta \right)\kappa \right)}{\mathit{\nu }}_n\mathrm{d}\mathrm{\Gamma }\mathrm{d}t. \end{gathered}$$

For problems with no design-dependent boundary conditions and unit volume heating source, the sensitivity presented in Eq. (7) may be expressed simply as [¹⁶,⁶¹]

(8){{\displaystyle \frac{\mathrm{d}\mathcal{J}}{\mathrm{d}s}}|}_{{\widehat{\theta }}^{\prime },\mathrm{\hspace{0.17em}}{\widehat{q}}^{\prime },\mathrm{\hspace{0.17em}}\nabla \widehat{q}\cdot \mathit{n},\mathrm{\hspace{0.17em}}{Q}^{\prime }=\mathrm{0}}={\displaystyle {\int }_{{\mathrm{\Gamma }}^s}\left({\displaystyle {\int }_{\mathrm{0}}^T\mathit{g}}\mathrm{d}t\right)}\cdot \mathbf{\nu }\mathrm{d}\mathrm{\Gamma },

where

(9)\mathit{g}=(\varsigma {\psi }_{\omega }+\rho c{\displaystyle \frac{\partial \theta }{\partial t}}\vartheta +k\nabla \theta \cdot \nabla \vartheta -Q\vartheta -{\mathit{\gamma }}_{e}h\left(\varsigma {\psi }_{\mathit{\gamma },q}-\vartheta \right)\left(\nabla \theta \cdot \mathit{n}\right)+\varsigma \left({\psi }_{\mathit{\gamma },\theta }\left(\nabla \theta \cdot \mathit{n}\right)-{\psi }_{\mathit{\gamma }}\kappa \right){\displaystyle \frac{\partial \theta }{\partial t}}-q\nabla \vartheta \cdot \mathit{n}+\left(q\vartheta \right)\kappa )\mathit{n},

with

(10){\mathit{\gamma }}_e=\{\begin{array}{ll}\mathrm{1},\hfill & \mathrm{\hspace{1em}}\forall \mathit{x}\in {\mathrm{\Gamma }}_e,\hfill \\ \mathrm{0},\hfill & \mathrm{\hspace{1em}}\mathrm{otherwise}\mathrm{.}\hfill \end{array}

It should be noted that with slightly proper modifications, the presented shape sensitivity framework can be applicable to the popular level-set-based topology optimizations [⁶²,⁶³].

IGA for transient heat conduction problems

The basic idea of IGA is to discretize the temperature field using non-uniform rational B-splines (NURBS) shape functions R^I, i.e.,

(11)\theta ={\displaystyle \sum _I{R}^I{\theta }^I}=\mathit{R}\cdot \mathit{\theta },

with {\theta }^I as the temperature control variables, such that a system of discrete equations can be derived from the weak formulation of Eq. (3) or (6) as

(12)\mathit{C}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}+\mathit{K}\mathit{\theta }=\mathit{f},

in which \mathit{C} and \mathit{K} are capacitance and conductance matrices, respectively, \mathit{\theta } is the vector of unknown nodal temperature, and \mathit{f} is heat flux vector. This discretization approach is the basic concept of IGA, which integrates CAD modeling and finite element analysis (FEM).

Introducing a parameter \beta \in [\mathrm{0},\mathrm{1}], the temperature discretization of Eq. (12) in the time domain can be written as

(13)\mathit{C}{\displaystyle \frac{{\mathit{\theta }}^{t+\mathrm{\Delta }t}-{\mathit{\theta }}^t}{\mathrm{\Delta }t}}+\mathit{K}\left(\beta {\mathit{\theta }}^{t+\mathrm{\Delta }t}+\left(\mathrm{1}-\beta \right){\mathit{\theta }}^t\right)=\mathit{f}.

The special cases of β =0, 0.5, and 1, correspond to fully explicit (forward Euler), semi-implicit (Crank–Nicolson-type) and fully implicit (backward Euler) schemes, respectively. The transient primary and adjoint problems are solved with this isogeometric framework in order to determine the fields required for shape sensitivity analysis.

Numerical error induced by time integration schemes

The numerical error of adjoint shape sensitivity analysis from structural analysis is induced mainly by the discretization in space to approximate the temperature field using polynomial basis functions, and in time due to the time integration schemes. The spatial numerical error can be reduced by decreasing the mesh size at a cost of increasing computational time. To investigate the extent of numerical error arising from time discretization in a transient sensitivity analysis, it is necessary to discuss on the accuracy, stability and oscillations induced by the time integration scheme.

Accuracy

The time discretization scheme used in Eq. (13) is a generalized trapezoidal rule. The temperature variation in a time-step is approximated as

(14){\displaystyle \frac{{\theta }^{t+\mathrm{\Delta }t}-{\theta }^t}{\mathrm{\Delta }t}}=\left(\mathrm{1}-\beta \right){\dot{\theta }}^t+\beta {\dot{\theta }}^{t+\mathrm{\Delta }t}.

Using Taylor series expansion, we have

(15){\theta }^{t+\mathrm{\Delta }t}={\theta }^t+\mathrm{\Delta }t{\dot{\theta }}^t+{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathrm{2}}}{\mathrm{2}}}{\ddot{\theta }}^t+\mathcal{O}\left[\mathrm{\Delta }{t}^{\mathrm{3}}\right],

and

(16){\dot{\theta }}^{t+\mathrm{\Delta }t}={\dot{\theta }}^t+\mathrm{\Delta }t{\ddot{\theta }}^t+\mathcal{O}\left[\mathrm{\Delta }{t}^{\mathrm{2}}\right].

From Eq. (15), it can be obtained that

(17){\displaystyle \frac{{\theta }^{t+\mathrm{\Delta }t}-{\theta }^t}{\mathrm{\Delta }t}}={\dot{\theta }}^t+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\ddot{\theta }}^t+\mathcal{O}\left[\mathrm{\Delta }{t}^{\mathrm{2}}\right].

The truncation error can thus be obtained from the difference between Eqs. (16) and (17) as

(18){\displaystyle \frac{{\theta }^{t+\mathrm{\Delta }t}-{\theta }^t}{\mathrm{\Delta }t}}-\left(\left(\mathrm{1}-\beta \right){\dot{\theta }}^t+\beta {\dot{\theta }}^{t+\mathrm{\Delta }t}\right)=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-\beta \right)\mathrm{\Delta }t{\ddot{\theta }}^t+\mathcal{O}\left[\mathrm{\Delta }{t}^{\mathrm{2}}\right].

This indicates that (i) the numerical error of the temperature rate \dot{\theta } is in the order of \mathrm{\Delta }t when \beta \ne \mathrm{0.5} and \mathrm{\Delta }{t}^{\mathrm{2}} when \beta =\mathrm{0.5}; and (ii) as \beta approaches 0.5, the magnitude of truncation error generally becomes smaller.

In general, if oscillations do not occur in the numerical solution, the Crank–Nicolson scheme (\beta =\mathrm{0.5}) induces the lowest truncation error. However, when the time-step is large, oscillating solutions may develop, which may induce significant errors.

Stability

The truncation error at each time-step accumulates in time. This leads to the issue of stability, of which the modal analysis is able to give an indication. In modal analysis, a separation of variable is imposed on the time-dependent temperature field \mathit{\theta }[\mathit{x},t] in terms of an orthonormal basis defined in space, {\mathit{\tau }}_i\left[\mathit{x}\right] and the corresponding temporal basis functions {\phi }_i\left[t\right], such that,

(19)\mathit{\theta }[\mathit{x},\mathrm{\hspace{0.17em}}t]={\displaystyle \sum _{i=\mathrm{1}}^n{\phi }_i}\left[t\right]{\mathit{\tau }}_i\left[\mathit{x}\right],

in which n corresponds to the dimension of matrices \mathit{K} and \mathit{C}.

Bearing in mind that the orthonormal basis {\mathit{\tau }}_i\left[\mathit{x}\right] satisfy

(20)\forall i\ne j,\begin{array}{l}\hfill \end{array}{\mathit{\tau }}_i^{\mathrm{T}}\mathit{K}{\mathit{\tau }}_j={\mathit{\tau }}_i^{\mathrm{T}}\mathit{C}{\mathit{\tau }}_j=\mathrm{0},

and substituting Eq. (19) to Eq. (12), we obtain

(21){\mathit{\tau }}_j^{\mathrm{T}}\mathit{C}{\mathit{\tau }}_j\dot{\phi }\left[t\right]+{\mathit{\tau }}_j^{\mathrm{T}}\mathit{K}{\mathit{\tau }}_j\phi \left[t\right]={\mathit{\tau }}_j^{\mathrm{T}}\mathit{f}.

Denoting

(22)c_j={\mathit{\tau }}_j^{\mathrm{T}}\mathit{C}{\mathit{\tau }}_j,\begin{array}{l}\hfill \end{array}{k}_j={\mathit{\tau }}_j^{\mathrm{T}}\mathit{K}{\mathit{\tau }}_j,\begin{array}{l}\hfill \end{array}\mathrm{and}\begin{array}{l}\hfill \end{array}{f}_j={\mathit{\tau }}_j^{\mathrm{T}}\mathit{f},

we have

(23)c_j{\dot{\phi }}_j\left[t\right]+k_j{\phi }_j\left[t\right]=f_j\left[t\right],

in which c_j, k_j, and f_j\left[t\right] are the generalized heat capacitance, generalized heat conductance and generalized heat flux for the jth natural mode, respectively. Compared to the coupled equations in Eq. (12), Eq. (23) is decoupled and thus facilitates stability analysis.

Assuming the temperature field \mathit{\theta }[\mathit{x},t] in the form of \mathit{\theta }[\mathit{x},t]=\mathit{\tau }\left[\mathit{x}\right]{\mathrm{e}}^{-\alpha t} with e as the Napier’s constant (\mathrm{e}\approx2.71828), and ignoring the load vector (external excitation) \mathit{f}, the system of equations in Eq. (12) becomes

(24)\mathit{C}\dot{\mathit{\theta }}+\mathit{K}\mathit{\theta }=\left(\mathit{K}-\alpha \mathit{C}\right)\mathit{\tau }{\mathrm{e}}^{-\alpha t}=\mathbf{0}.

For the global capacitance matrix \mathit{C}, there exists a lower triangular matrix \mathit{L} such that \mathit{C}=\mathit{L}{\mathit{L}}^{\mathrm{T}}. From Eq. (24), we obtain the following eigen problem

(25)\left(\overline{\mathit{K}}-\alpha \mathit{I}\right){\mathit{L}}^{\mathrm{T}}\mathit{\tau }=\mathrm{0},\mathrm{with}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\overline{\mathit{K}}={\mathit{L}}^{-\mathrm{1}}\mathit{K}{\mathit{L}}^{-\mathrm{T}}.

The n\times n matrix \overline{\mathit{K}} has n real non-negative eigenvalues ({\alpha }_{\mathrm{1}},\mathrm{\hspace{0.17em}}{\alpha }_{\mathrm{2}},\mathrm{\hspace{0.17em}}\mathrm{...},\mathrm{\hspace{0.17em}}{\alpha }_n) that correspond to n eigenvectors {\overline{\mathit{\tau }}}_i={\mathit{L}}^{\mathrm{T}}{\mathit{\tau }}_k. One can easily derive that \mathit{K}{\mathit{\tau }}_i={\alpha }_i\mathit{C}{\mathit{\tau }}_i. Referring to Eq. (22), it can be obtained that

(26)a_j={\displaystyle \frac{k_j}{c_j}}.

Eventually, Eq. (23) can be re-written as

(27){\dot{\phi }}_j\left[t\right]+{\alpha }_j{\phi }_j\left[t\right]={\displaystyle \frac{f_j\left[t\right]}{c_j}}.

Substituting the time differentiating scheme in Eq. (13) into Eq. (27),

(28){\phi }_j^{t+\mathrm{\Delta }t}-{\phi }_j^t+{\alpha }_j\mathrm{\Delta }t\left(\beta {\phi }_j^{t+\mathrm{\Delta }t}+\left(\mathrm{1}-\beta \right){\phi }_j^t\right)={\displaystyle \frac{\mathrm{\Delta }t{f}_j}{c_j}},

which can be rearranged to give

(29){\phi }_j^{t+\mathrm{\Delta }t}={\displaystyle \frac{\mathrm{1}-\left(\mathrm{1}-\beta \right){\alpha }_j\mathrm{\Delta }t}{\mathrm{1}+\beta {\alpha }_j\mathrm{\Delta }t}}{\phi }_j^t+\frac{\mathrm{\Delta }t}{\mathrm{1}+\beta {\alpha }_j\mathrm{\Delta }t}{\displaystyle \frac{f_j}{c_j}}.

To ensure that {\phi }_j^{t+\mathrm{\Delta }t} remains bounded in time, the recurrence factor should satisfy

(30)\left|{\displaystyle \frac{\mathrm{1}-\left(\mathrm{1}-\beta \right){\alpha }_j\mathrm{\Delta }t}{\mathrm{1}+\beta {\alpha }_j\mathrm{\Delta }t}}\right|\le \mathrm{1},\mathrm{\hspace{0.17em}}\forall j.

Note that the eigenvalues of a positive matrix are positive. The stability condition can be obtained as

(31)\left(\mathrm{1}-\mathrm{2}\beta \right){\alpha }_j\mathrm{\Delta }t\le \mathrm{2},\mathrm{\hspace{0.17em}}\forall j,

which indicates that the recurrence in Eq. (29) is unconditionally stable if \mathrm{0.5}\le \beta \le \mathrm{1}, and conditionally stable if . For the conditionally stable cases, the critical time-step for stability is given by

(32)\mathrm{\Delta }{t}_{\mathrm{sta}}={\displaystyle \frac{\mathrm{2}}{\left(\mathrm{1}-\mathrm{2}\beta \right)\overline{\alpha }}},

in which \overline{\alpha }=\min \left({\alpha }_j\right), \forall j.

Oscillations

Non-physical oscillatory results may occur depending on the discretization in time and space. For the temporal discretization, when the recurrence factor in Eq. (29) is negative, i.e., , the solution oscillates. The oscillations can reduce the accuracy for the thermal analysis and thus lead to a large error for the sensitivity analysis, especially when time-step size is large. The critical time-step size of a given eigenvalue {\alpha }_j is \mathrm{1}-(\mathrm{1}-\beta ){\alpha }_j\mathrm{\Delta }t>\mathrm{0}, i.e.,

(33)\mathrm{\Delta }t\le \mathrm{\Delta }{t}_{\mathrm{osc}}={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-\beta \right){\alpha }_j}}.

In literature such as Refs. [⁶⁴,⁶⁵], it is suggested to use a time-step size smaller than {\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-\beta \right)\overline{\alpha }}} with \overline{\alpha }\ge {\alpha }_j, \forall \mathrm{\hspace{0.17em}}j, to avoid the oscillatory solutions. However, this is not always true. Oscillatory results may still develop when smaller time-step sizes are used, which will be demonstrated later. In a typical modal analysis, the overall thermal responses are dominated by the low frequency modes corresponding to low eigenvalues {\alpha }_j. This is also indicated in the recurrence factor in Eq. (29) where a larger recurrence factor is associated with the low frequency modes (smaller {\alpha }_j). A time-step size that guarantees a certain number of non-oscillatory low frequency modes can produce a reasonably good solution since the contributions of the high frequency modes may not be too significant.

Oscillatory solutions cannot be avoided simply by using smaller time-step sizes, because the spatial discretization have an important role as well. For a continuum media, the solution comprises of an infinite number of frequencies. However, in a numerical approach, spatial discretization limits the extent where the high frequency modes can be captured. If the time-step is in the order of the (high) frequencies that cannot be captured by the spatial discretization, oscillations will occur. This problem is also explained in terms of non-dimensionalization of the transient heat conduction equation in Eq. (3) (see Ref. [⁶⁶]) and it is recommended to choose a time-step size using

(34)\mathrm{\Delta }t\ge {\displaystyle \frac{\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}\rho c}{k}},

with \mathrm{\Delta }\mathit{x} as characteristic spatial mesh size. These oscillations caused by the high frequencies are demonstrated in the numerical example later. Nevertheless, in most cases, the contributions of the high frequency modes are negligible, particularly so with sufficiently small spatial mesh size.

Numerical studies

Minimum boundary problem

Problem description

The plate shown in Fig. 2(a) is first heated up to a given temperature and then left in a low-temperature environment to cool down. It is desired to minimize the heat dissipation speed by varying boundary {\mathrm{\Gamma }}_{\mathrm{2}} for the same amount material used. The optimal solution for the minimum boundary problem is to have a circular outer boundary {\mathrm{\Gamma }}_{\mathrm{2}}. The problem can be expressed as

(35)\mathcal{J}={\displaystyle {\int }_{\mathrm{0}}^T{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{2}}}h\left(\theta -{\theta }_e\right)\mathrm{d}\mathrm{\Gamma }}}\mathrm{d}t,

with the volume fixed.

Fig.2 The initial plate design and the NURBS parameterization (values in m) [¹⁶,²⁰]. Problem parameters are: {\theta }_{\mathrm{0}}\left[\mathit{x}\right]=\mathrm{100}°C, \forall \mathit{x}\in {\mathrm{\Omega }}^s, {\theta }_e=\mathrm{0}°C, \rho =\mathrm{7800} kg/m³, c=\mathrm{420} J/(kg·°C), k=\mathrm{20} W/(m·°C) and h=\mathrm{50} W/(m²·°C).

Full size|PPT slide

The design time is chosen to be T=\mathrm{300} s. For simplicity, only a quarter of the plate is parametrized, as depicted in Fig. 2(b), using NURBS with knot vectors \mathit{\xi }=\left[\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{0}\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\mathrm{\hspace{0.17em}}\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{1}\right] and \mathit{\eta }=\left[\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{1}\mathrm{\hspace{0.17em}}\mathrm{1}\right], and control points shown in Table 1. The locations of six control points, denoted as C_I, I = 1, 2, …, 6, are chosen as design variables. Due to symmetry, control point‚C_{\mathrm{1}} is restricted to move only horizontally while C_{\mathrm{6}} only vertically.

Stability and oscillations in the transient analysis

The analysis model is refined with standard k-refinement approach using knot vectors \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{20}}},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{2}}{\mathrm{20}}},\mathrm{\hspace{0.17em}}\mathrm{...},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{19}}{\mathrm{20}}}\right] and \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{10}}},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{2}}{\mathrm{10}}},\mathrm{\hspace{0.17em}}\mathrm{...},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{9}}{\mathrm{10}}}\right] in the two orthonormal index directions, respectively. This eventually produce an isogeometric model with 288 control points, leading to matrices \mathit{C} and \mathit{K} with a dimension of \mathrm{288}\times \mathrm{288}. Following Eq. (25), the largest eigenvalue for matrix \overline{\mathit{K}} is 422.22. The critical time-step size of the stability condition in Eq. (32) versus the time integration scheme coefficient \beta is plotted in Fig. 3, which shows that the time-step size needs to be very small to ensure the analysis stability for . When the time-step size \mathrm{\Delta }t>\mathrm{0.00948} s with \beta =\mathrm{0.25}, the transient analysis of the heat conduction will become unbounded, which leads to failure of the sensitivity analysis.

Fig.3 Critical time-step size of the stability condition versus the time integration scheme coefficient \beta.

Full size|PPT slide

The 288 eigenvalues of matrix \overline{\mathit{K}} are depicted in Fig. 4(a) with the corresponding critical time-step sizes for the oscillatory conditions in Eq. (33) of \beta =0.5 and \beta =0.75. A zoom in on the eigenvalues smaller than 20 are shown in Fig. 4(b). It is obvious that the critical time-step size with \beta =0.75 is twice of that with \beta =0.5. When time-step size is \mathrm{\Delta }t=\mathrm{1} s, the non-oscillatory eigenmodes correspond to for \beta =0.5 and for \beta =0.75. When time-step size is \mathrm{\Delta }t=\mathrm{0.1} s, the non-oscillatory eigenmodes corresponds to for \beta =0.5 and for \beta =0.75. The temperature histories of initial time-steps at points A and C_{\mathrm{1}} with different time-step sizes are plotted in Fig. 5 for \beta =\mathrm{0.5} and Fig. 6 for \beta =\mathrm{0.75}, respectively. It can be observed that (i) the oscillations of temperature for \beta =\mathrm{0.5} are much more noticeable than those for \beta =\mathrm{0.75}; and (ii) even with a time-step size that is smaller than the critical values of all eigenvalues, the oscillations still persist, as depicted in Figs. 5(n) and 6(n) and explained in Section 4.3.

Fig.4 Critical time-step size of the oscillatory conditions for (a) all 288 eigenvalues and (b) the first 87 eigenvalues smaller than 20, with \beta =0.5 and 0.75, respectively.

Full size|PPT slide

Fig.5 Temperature oscillations at point A and C_{\mathrm{1}} with different time-step sizes and \beta =0.5 for the first few iterative steps.

Full size|PPT slide

Fig.6 Temperature oscillations at point A and C_{\mathrm{1}} with different time-step sizes and \beta =0.75 for the first few iterative steps.

Full size|PPT slide

Tab.1 Initial locations of the design control points for the minimum boundary problem [¹⁶,²⁰]

I\left(i,j\right)	Location		Weight	‚I\left(i,j\right)	Location		Weight
	{\mathit{x}}_{\mathrm{1}}^I	{\mathit{x}}_{\mathrm{2}}^I			{\mathit{x}}_{\mathrm{1}}^I	{\mathit{x}}_{\mathrm{2}}^I
(1, 1)	0.0100	0.0000	1.00	‚(4, 2)	0.0091	0.0121	0.85
(2, 1)	0.0100	0.0026	0.90	‚(5, 2)	0.0039	0.0150	0.90
(3, 1)	0.0080	0.0061	0.85	‚(6, 2)	0.0000	0.0150	1.00
(4, 1)	0.0061	0.0080	0.85	‚(1, 3)	0.0200	0.0000	1.00
(5, 1)	0.0026	0.0100	0.90	‚(2, 3)	0.0200	0.0100	1.00
(6, 1)	0.0000	0.0100	1.00	‚(3, 3)	0.0238	0.0213	1.00
(1, 2)	0.0150	0.0000	1.00	‚(4, 3)	0.0213	0.0238	1.00
(2, 2)	0.0150	0.0039	0.90	‚(5, 3)	0.0100	0.0200	1.00
(3, 2)	0.0121	0.0091	0.85	‚(6, 3)	0.0000	0.0200	1.00

Referential sensitivity calculation

Using a relatively small time-step \mathrm{\Delta }t=0.01 s, the sensitivity analysis, termed {\mathit{G}}_f^I for design control point {\mathit{x}}^I, is computed using finite difference (FD) method for β=0.5, 0.75, and 1, respectively:

(36){\mathit{G}}_f^I:={\displaystyle \frac{\mathcal{J}\left[{\mathit{x}}^I+\delta \mathit{x}\right]-\mathcal{J}\left[{\mathit{x}}^I\right]}{\delta \mathit{x}}},

where \delta \mathit{x} is the perturbation of the locations of design control points. The FD computation for is omitted due to the unrealistic solutions caused by the instability.

The results presented in Table 2, show a close match between the three cases. The referential sensitivity analysis is calculated using the average value:

(37){\overline{\mathit{G}}}_f^I={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left({\mathit{G}}_{f\mathrm{1}}^I+{\mathit{G}}_{f\mathrm{2}}^I+{\mathit{G}}_{f\mathrm{3}}^I\right),

where {\mathit{G}}_{f\mathrm{1}}^I, {\mathit{G}}_{f\mathrm{2}}^I, and {\mathit{G}}_{f\mathrm{3}}^I are the FD gradient for β=0.5, 0.75, and 1, respectively. The relative difference of {\mathit{G}}_{fi} comparing to the referential sensitivity analysis is calculated using

(38){\mathit{D}}_{fi}={\displaystyle \frac{{\mathit{G}}_{fi}-{\overline{\mathit{G}}}_f}{\max {\overline{\mathit{G}}}_f}}\mathrm{.}

The L^{\mathrm{2}} norm of {\mathit{D}}_{fi} shows the magnitude of the difference between calculated sensitivity and the referential sensitivity. The L^{\mathrm{2}} norms of the relative difference for {\mathit{G}}_{f\mathrm{1}}, {\mathit{G}}_{f\mathrm{2}}, and {\mathit{G}}_{f\mathrm{3}} are \mathrm{1.2924}\times {\mathrm{10}}^{-\mathrm{5}}, \mathrm{1.6588}\times {\mathrm{10}}^{-\mathrm{5}}, and \mathrm{2.5184}\times {\mathrm{10}}^{-\mathrm{5}}, respectively, which provides a quantification of the close match between the three cases.

Tab.2 Sensitivity analysis using FD with \mathrm{\Delta }t=\mathrm{0.01} s for different \beta and the referential sensitivity of the minimum boundary problem

C_I	Component	FD			Referential FD
		\beta =\mathrm{0.5}	\beta =\mathrm{0.75}	\beta =\mathrm{1}
C_{\mathrm{1}}	1	135874.5067	135876.4566	135887.5961	135879.5198
C_{\mathrm{2}}	1	243335.5767	243333.5103	243350.4778	243339.8549
	2	–37170.7865	–37162.2082	–37166.9958	–37166.6635
C_{\mathrm{3}}	1	897323.9019	897326.7759	897308.0330	897319.5703
	2	306189.0384	306194.4444	306182.1844	306188.5557
C_{\mathrm{4}}	1	306189.1548	306194.2916	306176.7857	306186.7440
	2	897323.3635	897312.5441	897306.3377	897314.0818
C_{\mathrm{5}}	1	–37165.4896	–37167.1631	–37163.8816	–37165.5115
	2	243336.8791	243328.4244	243347.0800	243337.4611
C_{\mathrm{6}}	2	135881.9864	135871.1306	135886.7594	135879.9588

Numerical adjoint sensitivity analysis convergence with respect to the number of time-steps for different \beta

For β=0.5, 0.75, and 1, the adjoint shape sensitivity is computed for time-step sizes \mathrm{\Delta }t=60, 30, 15, 10, 5, 3, 2, 1, 0.5, 0.3, 0.1, and 0.01 s, corresponding to 5, 10, 20, 30, 60, 100, 150, 300, 600, 1000, 3000, and 30000 time-steps, respectively. The L^{\mathrm{2}} norms of the relative difference for the different numbers of time-steps are depicted in Fig. 7, showing all the three cases converge to a small value. It is also clear that the case of \beta =\mathrm{0.75} has an overall better performance than other two cases. The larger error associated with \beta =\mathrm{0.5} is mainly due to the oscillatory solutions induced by the time integration scheme, as demonstrated earlier in Section 5.1.2.

Fig.7 The L^{\mathrm{2}} norm of the relative difference of the adjoint sensitivity analysis versus number of time-steps for the minimum boundary problem.

Full size|PPT slide

The computational time of the adjoint sensitivity analysis for each time-step size on a Dell laptop with a four core CPU of i7-6600U is listed in Table 3. It can be found that as the step-size decreases, the computational cost increases dramatically, mainly caused by the increase of analysis time. This indicates that using time integration schemes with \beta<0.5 is not preferable as it requires a very smaller time-step size to keep the analysis stability.

Tab.3 Computational time of different time-step sizes for the minimum boundary problem

Time-step size/s	Number of time-steps	Computational time/s
60.00	5	0.4633
30.00	10	0.5676
15.00	20	0.6906
10.00	30	0.7631
5.00	60	1.1828
3.00	100	1.7065
2.00	150	2.0199
1.00	300	3.7158
0.50	600	6.8297
0.30	1000	11.2075
0.10	3000	36.3381
0.01	30000	429.7075

Numerical error of adjoint sensitivity analysis with respect to \beta at a given time interval discretization

To further investigate the numerical error of adjoint shape sensitivity with respect to \beta, the adjoint shape sensitivity is computed for \mathrm{\Delta }t=10, 3, and 0.3 s (corresponding to 30, 100, and 1000 time-steps), respectively, with different values of \beta ranging from 0.5 to 1. The numerical error versus coefficient \beta is plotted for these three cases with different time-step sizes in Fig. 8. For all three cases, it is observed that the numerical error increases when the value of \beta increases from 0.52 to 1. When \beta =\mathrm{0.5}, the numerical errors for the two coarser time discretizations are relatively big due to the oscillatory solutions. A larger value of \beta can help to reduce the error of the shape sensitivity analysis caused by the oscillatory solutions.

Fig.8 The L^{\mathrm{2}} norm of the relative difference of the adjoint sensitivity analysis versus \beta for the minimum boundary problem.

Full size|PPT slide

Plunger shape design problem

Problem description

Consider a plunger designed to form a television glass bulb panel as listed in Fig. 9 [¹⁶,²⁰,⁶⁷]. The model is parametrized with knot arrays \mathit{\xi }=[\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{0.2},\mathrm{\hspace{0.17em}}\mathrm{0.4},\mathrm{\hspace{0.17em}}\mathrm{0.5},\mathrm{0.7},\mathrm{\hspace{0.17em}}\mathrm{1},\mathrm{\hspace{0.17em}}\mathrm{1},\mathrm{\hspace{0.17em}}\mathrm{1}] and \mathit{\eta }=[\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{0},\mathrm{\hspace{0.17em}}\mathrm{1},\mathrm{\hspace{0.17em}}\mathrm{1},\mathrm{\hspace{0.17em}}\mathrm{1}], and control points listed in Table 4. The heat convection coefficient are h_{\mathrm{1}}=\mathrm{3.15}\times {\mathrm{10}}^{-\mathrm{4}} W/(mm²·°C) on {\mathrm{\Gamma }}_{\mathrm{1}} and h_{\mathrm{3}}=\mathrm{2.88}\times {\mathrm{10}}^{-\mathrm{4}} W/(mm²·°C) on {\mathrm{\Gamma }}_{\mathrm{3}}. Other related parameters are k=\mathrm{27.52}\times {\mathrm{10}}^{-\mathrm{3}} W/(mm²·°C), \rho c=\mathrm{2.288}\times {\mathrm{10}}^{-\mathrm{3}} J/(mm³·°C) and T=\mathrm{500} s.

Fig.9 NURBS parameterization of the initial plunger model (values in mm) [¹⁶,²⁰].

Full size|PPT slide

Tab.4 Initial locations of the design control points for the plunger design problem [¹⁶,²⁰]

I\left(i,j\right)	Location		Weight	‚I\left(i,j\right)	Location		Weight
	{\mathit{x}}_{\mathrm{1}}^I	{\mathit{x}}_{\mathrm{2}}^I			{\mathit{x}}_{\mathrm{1}}^I	{\mathit{x}}_{\mathrm{2}}^I
(1, 1)	0.00	100.00	1.00	‚(5, 2)	90.00	20.00	1.00
(2, 1)	0.00	80.00	1.00	‚(6, 2)	145.00	20.00	1.00
(3, 1)	0.00	30.00	1.00	‚(7, 2)	200.00	20.00	1.00
(4, 1)	0.00	0.00	0.71	‚(1, 3)	30.00	100.00	1.00
(5, 1)	30.00	0.00	1.00	‚(2, 3)	30.00	80.00	1.00
(6, 1)	140.00	0.00	1.00	‚(3, 3)	30.00	65.00	1.00
(7, 1)	200.00	0.00	1.00	‚(4, 3)	30.00	45.00	1.00
(1, 2)	15.00	100.00	1.00	‚(5, 3)	70.00	45.00	1.00
(2, 2)	15.00	80.00	1.00	‚(6, 3)	120.00	45.00	1.00
(3, 2)	15.00	65.00	1.00	‚(7, 3)	200.00	45.00	1.00
(4, 2)	15.00	20.00	1.00

The ambient temperature of boundary {\mathrm{\Gamma }}_{\mathrm{3}}, which contact the molten glass, is assumed to be {\theta }_{e\mathrm{3}}=\mathrm{1000}°C, while the ambient temperature of boundary {\mathrm{\Gamma }}_{\mathrm{1}}, which contact the cooling fluid, is assumed to be {\theta }_{e\mathrm{1}}=\mathrm{0}°C. Temperature difference along the fixed boundary {\mathrm{\Gamma }}_{\mathrm{3}} affects the quality of the television bulb, which makes it necessary to design boundary {\mathrm{\Gamma }}_{\mathrm{1}} such that the temperature difference along boundary {\mathrm{\Gamma }}_{\mathrm{3}} can be minimized. For this problem, an objective functional is introduced as

(39)\mathcal{J}={\displaystyle {\int }_{\mathrm{0}}^T{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{3}}}{(\theta -\tilde{\theta })}^{\mathrm{2}}\mathrm{d}\mathrm{\Gamma }}}\mathrm{d}t,

where the average temperature along {\mathrm{\Gamma }}_{\mathrm{3}}, \tilde{\theta }\left[t\right], is computed using

(40)\tilde{\theta }={\displaystyle \frac{\mathrm{1}}{\left|{\mathrm{\Gamma }}_{\mathrm{3}}\right|}}{\displaystyle {\int }_{{\mathrm{\Gamma }}_{\mathrm{3}}}\theta \mathrm{d}\mathrm{\Gamma }}.

Stability and oscillations in the transient analysis

The analysis model is refined with standard k-refinement approach using knot vectors \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{50}}},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{2}}{\mathrm{50}}},\mathrm{\hspace{0.17em}}\mathrm{...},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{49}}{\mathrm{50}}}\right] and \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{2}}{\mathrm{8}}},\mathrm{\hspace{0.17em}}\mathrm{...},\mathrm{\hspace{0.17em}}{\displaystyle \frac{\mathrm{7}}{\mathrm{8}}}\right] in the two orthonormal index directions, respectively. This eventually produces an isogeometric model with 520 control points and matrices \mathit{C} and \mathit{K} with a dimension of \mathrm{520}\times \mathrm{520}. Following Eq. (25), the maximum eigenvalue for matrix \overline{\mathit{K}} is 62.67. The critical time-step size of the stability condition in Eq. (32) versus the time integration scheme coefficient \beta is plotted in Fig. 10, where it can be observed that the time-step size needs to be very small to ensure the analysis stability for . When the time-step size \mathrm{\Delta }t>\mathrm{0.01778} s with \beta =\mathrm{0.25}, the transient analysis of the heat conduction will become unbounded, which leads to failure of the sensitivity analysis.

Fig.10 Critical time-step size of the stability conditions versus the time integration scheme coefficient \beta.

Full size|PPT slide

The 520 eigenvalues of matrix \overline{\mathit{K}} are depicted in Fig. 11(a) with the corresponding critical time-step sizes for the oscillatory conditions Eq. (32) of \beta =0.5 and 0.75. A zoom in on the eigenvalues smaller than 20 are shown in Fig. 11(b). It is obvious that the critical time-step sizes of the majority non-oscillatory eigenmodes are bigger than 0.1 s.

Fig.11 Critical time-step size of the oscillatory conditions for (a) all 520 eigenvalues and (b) the first 445 eigenvalues smaller than 20, with \beta =0.5 and 0.75, respectively.

Full size|PPT slide

Referential sensitivity analysis calculation

Similarly, the FD sensitivity is computed for β=0.5, 0.75, and 1 with a perturbation of \delta {\mathit{x}}_i={\mathrm{10}}^{-\mathrm{6}} and time-step size of \mathrm{\Delta }t=\mathrm{0.05} s. The FD computation for is omitted due to the unrealistic solutions caused by the instability.

The results are presented in Table 5, which shows that the sensitivities of these four cases are relatively close. The referential sensitivity analysis is calculated using the same approach as in Section 5.1.3. The L^{\mathrm{2}} norm of the relative difference for {\mathit{G}}_{f\mathrm{1}}, {\mathit{G}}_{f\mathrm{2}}, and {\mathit{G}}_{f\mathrm{3}} are \mathrm{0.9862}\times {\mathrm{10}}^{-\mathrm{4}}, \mathrm{0.2142}\times {\mathrm{10}}^{-\mathrm{4}}, and \mathrm{1.0319}\times {\mathrm{10}}^{-\mathrm{4}}, respectively, which conform the close match of the three cases and the average gradient {\overline{\mathit{G}}}_f can be used as the referential gradient.

Tab.5 Sensitivity analysis using FD with \Delta t=\mathrm{0.05} s for different \beta and the referential sensitivity of the plunger design problem

C_I	Component	FD			Referential FD
		\beta =\mathrm{0.5}	\beta =\mathrm{0.75}	\beta =\mathrm{1}
C_{\mathrm{1}}	1	–6.2906×105	–6.2891×105	–6.2885×105	–6.2894×105
	2	€2.0000×10	€3.0000×10	1.0000×10	2.0000×10
C_{\mathrm{2}}	1	–3.2439×105	–3.2434×105	–3.2422×105	–3.2432×105
	2	€€€3.4358×105	€€€3.4351×105	€€€3.4350×105	€€€3.4353×105
C_{\mathrm{3}}	1	€€€2.0760×106	€€€2.0759×106	€€€2.0758×106	€€€2.0759×106
	2	€€€1.1210×106	€€€1.1210×106	€€€1.1210×106	€€€1.1210×106
C_{\mathrm{4}}	1	€€€8.6030×104	€€€8.6030×104	€€€8.6010×104	€€€8.6020×104
	2	–6.2714×105	–6.2713×105	–6.2714×105	–6.2714×105
C_{\mathrm{5}}	1	€2.0000×10	€1.0000×10	1.0000×10	€2.0000×10
	2	€€€1.5853×106	€€€1.5852×106	€€€1.5849×106	€€€1.5851×106

Numerical adjoint sensitivity analysis convergence with respect to the number of time-steps for different \beta

For \beta= 0.5, 0.75, and 1, the adjoint shape sensitivity is computed for time-step sizes \mathrm{\Delta }t=25, 10, 5, 2.5, 1, 0.5, 0.25, 0.1, and 0.05 s, corresponding to 20, 50, 100, 200, 500, 1000, 2000, 5000, and 10000 time-steps, respectively. The L^{\mathrm{2}} norms of the relative difference versus the number of time-steps are plotted in Fig. 12, where it can be observed that all three cases converge to a relatively small value. It is also clear that the case of \beta =\mathrm{0.5} has an overall better performance than other two cases.

Fig.12 The L^{\mathrm{2}}norm of the relative difference of the adjoint sensitivity analysis versus number of time-steps for the plunger design problem.

Full size|PPT slide

The computational time of the adjoint sensitivity analysis for each time-step size on a Dell laptop with a four core CPU of i7-6600U is listed in Table 6. It can be found that as the step-size decreases, the computational cost increases dramatically, mainly caused by the increase of analysis time. Similarly, this indicates that using time integration schemes with \beta<0.5 is not preferable as it requires a very smaller time-step size to keep the analysis stability.

Tab.6 Computational time of different time step-sizes for the plunger design case

Time-step size/s	Number of time-steps	Computational time/s
25.00	20	11.5991
10.00	50	27.7114
5.00	100	54.6364
2.50	200	108.7865
1.00	500	278.4205
0.50	1000	565.3433
0.25	2000	1110.6680
0.10	5000	2794.8650
0.05	10000	5704.8800

Numerical error of adjoint sensitivity analysis with respect to \beta at a given time interval discretization

To further investigate the numerical error of adjoint shape sensitivity with respect to \beta, the adjoint shape sensitivity is computed for \mathrm{\Delta }t=10, 1, and 0.1 s (corresponding to 50, 500, and 5000 time-steps), respectively, with different values of \beta ranging from 0.48 to 1. The numerical error versus coefficient \beta is plotted for these three cases with different time-step sizes in Fig. 13. From Fig. 13, it can be seen that for all three cases, the numerical error increases when the value of \beta increases from 0.5 to 1.

Fig.13 The L^{\mathrm{2}} norm of the relative difference of the adjoint sensitivity analysis versus \beta for the plunger design problem.

Full size|PPT slide

Conclusions

In this work, we investigate the numerical error of time integration scheme in adjoint shape sensitivity analysis for transient heat conduction problems. The accuracy, stability and oscillations in transient analysis, which are the main causes of numerical errors in time integration, are briefly discussed. The study is computed using IGA for adjoint shape sensitivity analysis of two benchmark transient heat conduction problems with design-dependent boundary conditions. In general, time integration approaches with coefficient are not recommended due to numerical stability concerns; Crank–Nicolson approach with \beta =\mathrm{0.5} may induce large error because of oscillatory solutions; semi-implicit approaches with \beta >\mathrm{0.5} are preferred; and fully implicit approach with \beta =\mathrm{1} has a lower accuracy than the semi-implicit approaches. Hence, a value around of \beta \approx \mathrm{0.75} is recommended.

Acknowledgements

The authors would like to thank Dr. Dan Wang from Institute of High Performance Computing (IHPC), A*STAR for the communications related to this work.ƒ

Open Access

This article is licensed under Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format, as long as appropriate credit is given to the original author(s) and the source, a link is provided to the Creative Commons license, and any changes made are indicated.

Images or other third-party materials in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Cssommons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

References

1 Li Q, Steven G P, Querin O M, Shape and topology design for heat conduction by evolutionary structural optimization. International Journal of Heat and Mass Transfer, 1999, 42(17): 3361–3371

2 Xie G, Liu Y, Sunden B, Computational study and optimization of laminar heat transfer and pressure loss of double-layer microchannels for chip liquid cooling. Journal of Thermal Science and Engineering Applications, 2013, 5(1): 011004

3 Sigmund O, Torquato S. Design of materials with extreme thermal expansion using a three-phase topology optimization method. Journal of the Mechanics and Physics of Solids, 1997, 45(6): 1037–1067

4 Gao T, Zhang W, Zhu J, Topology optimization of heat conduction problem involving design-dependent heat load effect. Finite Elements in Analysis and Design, 2008, 44(14): 805–813

5 Iga A, Nishiwaki S, Izui K, Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection. International Journal of Heat and Mass Transfer, 2009, 52(11–12): 2721–2732

6 Yaji K, Yamada T, Kubo S, A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions. International Journal of Heat and Mass Transfer, 2015, 81: 878–888

7 Xia Q, Xia L, Shi T. Topology optimization of thermal actuator and its support using the level set based multiple-type boundary method and sensitivity analysis based on constrained variational principle. Structural and Multidisciplinary Optimization, 2018, 57(3): 1317–1327

8 Choi K K, Kim N H. Structural Sensitivity Analysis and Optimization 1: Linear Systems. New York: Springer, 2005

9 Dems K, Rousselet B. Sensitivity analysis for transient heat conduction in a solid body-Part I: External boundary modification. Structural Optimization, 1999, 17(1): 36–45

10 Dems K, Rousselet B. Sensitivity analysis for transient heat conduction in a solid body-Part II: Interface modification. Structural Optimization, 1999, 17(1): 46–54

11 Haftka R T, Shore C P. Approximation Methods for Combined Thermal/Structural Design. NASA Technical Paper 1428. 1979

12 Haftka R T. Techniques for thermal sensitivity analysis. International Journal for Numerical Methods in Engineering, 1981, 17(1): 71–80

13 Greene W H, Haftka R T. Computational aspects of sensitivity calculations in transient structural analysis. Computers & Structures, 1989, 32(2): 433–443

14 Gao Z Y, Grandhi R V. Sensitivity analysis and shape optimization for preform design in thermo-mechanical coupled analysis. International Journal for Numerical Methods in Engineering, 1999, 45(10): 1349–1373

15 Haftka R T, Malkus D S. Calculation of sensitivity derivatives in thermal problems by finite differences. International Journal for Numerical Methods in Engineering, 1981, 17(12): 1811–1821

16 Wang Z P, Turteltaub S, Abdalla M M. Shape optimization and optimal control for transient heat conduction problems using an isogeometric approach. Computers & Structures, 2017, 185: 59–74

17 Michaleris P, Tortorelli D A, Vidal C A. Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity. International Journal for Numerical Methods in Engineering, 1994, 37(14): 2471–2499

18 Tortorelli D A, Haber R B, Lu S C Y. Design sensitivity analysis for nonlinear thermal systems. Computer Methods in Applied Mechanics and Engineering, 1989, 77(1–2): 61–77

19 Tortorelli D A, Haber R B. First-order design sensitivities for transient conduction problems by an adjoint method. International Journal for Numerical Methods in Engineering, 1989, 28(4): 733–752

20 Wang Z P, Kumar D. On the numerical implementation of continuous adjoint sensitivity for transient heat conduction problems using an isogeometric approach. Structural and Multidisciplinary Optimization, 2017, 56(2): 487–500

21 Kane J H, Kumar B L, Stabinsky M. Transient thermoelasticity and other body force effects in boundary element shape sensitivity analysis. International Journal for Numerical Methods in Engineering, 1991, 31(6): 1203–1230

22 Jarny Y, Ozisik M N, Bardon J P. A general optimization method using adjoint equation for solving multidimensional inverse heat conduction. International Journal of Heat and Mass Transfer, 1991, 34(11): 2911–2919

23 Kleiber M, Służalec A. Material derivative and control volume approaches to shape sensitivity analysis of nonlinear transient thermal problems. Structural Optimization, 1996, 11: 56–63

24 Dorai G A, Tortorelli D A. Transient inverse heat conduction problem solutions via Newton’s method. International Journal of Heat and Mass Transfer, 1997, 40(17): 4115–4127

25 Korycki R. Sensitivity analysis and shape optimization for transient heat conduction with radiation. International Journal of Heat and Mass Transfer, 2006, 49(13–14): 2033–2043

26 Huang C H, Chaing M T. A transient three-dimensional inverse geometry problem in estimating the space and time-dependent irregular boundary shapes. International Journal of Heat and Mass Transfer, 2008, 51(21–22): 5238–5246

27 Służalec A, Kleiber M. Shape optimization of thermo-diffusive systems. International Journal of Heat and Mass Transfer, 1992, 35(9): 2299–2304

28 Gu Y X, Chen B S, Zhang H W, A sensitivity analysis method for linear and nonlinear transient heat conduction with precise time integration. Structural and Multidisciplinary Optimization, 2002, 24(1): 23–37

29 Chen B, Tong L. Sensitivity analysis of heat conduction for functionally graded materials. Materials & Design, 2004, 25(8): 663–672

30 Haftka R T, Grandhi R V. Structural shape optimization—A survey. Computer Methods in Applied Mechanics and Engineering, 1986, 57(1): 91–106

31 van Keulen F, Haftka R T, Kim N H. Review of options for structural design sensitivity analysis, Part 1: Linear systems. Computer Methods in Applied Mechanics and Engineering, 2005, 194(30–33): 3213–3243

32 Cho S, Ha S H. Isogeometric shape design optimization: Exact geometry and enhanced sensitivity. Structural and Multidisciplinary Optimization, 2009, 38(1): 53–70

33 Qian X. Full analytical sensitivities in nurbs based isogeometric shape optimization. Computer Methods in Applied Mechanics and Engineering, 2010, 199(29–32): 2059–2071

34 Nagy A P, Abdalla M M, Gürdal Z. Isogeometric sizing and shape optimisation of beam structures. Computer Methods in Applied Mechanics and Engineering, 2010, 199(17–20): 1216–1230

35 Nagy A P, Abdalla M M, Gürdal Z. Isogeometric design of elastic arches for maximum fundamental frequency. Structural and Multidisciplinary Optimization, 2011, 43(1): 135–149

36 Liu H, Yang D, Wang X, Smooth size design for the natural frequencies of curved Timoshenko beams using isogeometric analysis. Structural and Multidisciplinary Optimization, 2019, 59(4): 1143–1162

37 Weeger O, Narayanan B, Dunn M L. Isogeometric shape optimization of nonlinear, curved 3D beams and beam structures. Computer Methods in Applied Mechanics and Engineering, 2019, 345: 26–51

38 Nagy A P, IJsselmuiden S T, Abdalla M M. Isogeometric design of anisotropic shells: Optimal form and material distribution. Computer Methods in Applied Mechanics and Engineering, 2013, 264: 145–162

39 Hirschler T, Bouclier R, Duval A, Isogeometric sizing and shape optimization of thin structures with a solid-shell approach. Structural and Multidisciplinary Optimization, 2019, 59(3): 767–785

40 Lian H, Kerfriden P, Bordas S. Implementation of regularized isogeometric boundary element methods for gradient-based shape optimization in two-dimensional linear elasticity. International Journal for Numerical Methods in Engineering, 2016, 106(12): 972–1017

41 Lian H, Kerfriden P, Bordas S. Shape optimization directly from CAD: An isogeometric boundary element approach using T-splines. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 1–41

42 Wang C, Xia S, Wang X, Isogeometric shape optimization on triangulations. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 585–622

43 Wang Z P, Poh L H, Dirrenberger J, Isogeometric shape optimization of smoothed petal auxetic structures via computational periodic homogenization. Computer Methods in Applied Mechanics and Engineering, 2017, 323: 250–271

44 Wang Z P, Poh L H. Optimal form and size characterization of planar isotropic petal-shaped auxetics with tunable effective properties using IGA. Composite Structures, 2018, 201: 486–502

45 Kumar D, Wang Z P, Poh L H, Isogeometric shape optimization of smoothed petal auxetics with prescribed nonlinear deformation. Computer Methods in Applied Mechanics and Engineering, 2019, 356: 16–43

46 Wang Y, Benson D J. Geometrically constrained isogeometric parameterized level-set based topology optimization via trimmed elements. Frontiers of Mechanical Engineering, 2016, 11(4): 328–343

47 Wang Y, Xu H, Pasini D. Multiscale isogeometric topology optimization for lattice materials. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 568–585

48 Xie X, Wang S, Xu M, A new isogeometric topology optimization using moving morphable components based on R-functions and collocation schemes. Computer Methods in Applied Mechanics and Engineering, 2018, 339: 61–90

49 Lieu Q X, Lee J. Multiresolution topology optimization using isogeometric analysis. International Journal for Numerical Methods in Engineering, 2017, 112(13): 2025–2047

50 Hou W, Gai Y, Zhu X, Explicit isogeometric topology optimization using moving morphable components. Computer Methods in Applied Mechanics and Engineering, 2017, 326: 694–712

51 Liu H, Yang D, Hao P, Isogeometric analysis based topology optimization design with global stress constraint. Computer Methods in Applied Mechanics and Engineering, 2018, 342: 625–652

52 Hao P, Yuan X, Liu C, An integrated framework of exact modeling, isogeometric analysis and optimization for variable-stiffness composite panels. Computer Methods in Applied Mechanics and Engineering, 2018, 339: 205–238

53 Guo Y, Ruess M. Nitsche’s method for a coupling of isogeometric thin shells and blended shell structures. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 881–905

54 Cai S Y, Zhang W H, Zhu J, Stress constrained shape and topology optimization with fixed mesh: A B-spline finite cell method combined with level set function. Computer Methods in Applied Mechanics and Engineering, 2014, 278: 361–387

55 Zhang W, Zhao L, Gao T, Topology optimization with closed B-splines and Boolean operations. Computer Methods in Applied Mechanics and Engineering, 2017, 315: 652–670

56 Wang Y, Wang Z P, Xia Z, Structural design optimization using isogeometric analysis: A comprehensive review. Computer Modeling in Engineering & Sciences, 2018, 117(3): 455–507

57 Xia L, Xia Q, Huang X, Bi-directional evolutionary structural optimization on advanced structures and materials: A comprehensive review. Archives of Computational Methods in Engineering, 2018, 25(2): 437–478

58 Meng L, Zhang W, Quan D, From topology optimization design to additive manufacturing: Today’s success and tomorrow’s roadmap. Archives of Computational Methods in Engineering, 2019 (in press)

59 Kaminski W. Hyperbolic heat conduction equation for materials with a nonhomogeneous inner structure. Journal of Heat Transfer, 1990, 112(3): 555–560

60 Xia L, Breitkopf P. Recent advances on topology optimization of multiscale nonlinear structures. Archives of Computational Methods in Engineering, 2017, 24(2): 227–249

61 Wang Z P, Turteltaub S. Isogeometric shape optimization for quasi-static processes. International Journal for Numerical Methods in Engineering, 2015, 104(5): 347–371

62 Xia Q, Shi T, Liu S, A level set solution to the stress-based structural shape and topology optimization. Computers & Structures, 2012, 90: 55–64

63 Xia Q, Shi T, Xia L. Stable hole nucleation in level set based topology optimization by using the material removal scheme of BESO. Computer Methods in Applied Mechanics and Engineering, 2019, 343: 438–452

64 Reddy J N, Gartling D K. The Finite Element Method in Heat Transfer and Fluid Dynamics. Boca Raton: CRC Press, 2001

65 Bergheau J M, Fortunier R. Finite Element Simulation of Heat Transfer. Hoboken: John Wiley & Sons, 2013

66 Carter W C. Lecture Notes on Mathematics for Materials Science and Engineers. MIT 3.016, 2012

67 Ho Lee D, Man Kwak B. Shape sensitivity and optimization for transient heat diffusion problems using the BEM. International Journal of Numerical Methods for Heat & Fluid Flow, 1995, 5(4): 313–326