RESEARCH ARTICLE

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

  • Zhen-Pei WANG 1,2 ,
  • Zhifeng XIE 3 ,
  • Leong Hien POH , 1
Expand
  • 1. Department of Civil and Environmental Engineering, National University of Singapore, Singapore 117576, Singapore
  • 2. Institute of High Performance Computing (IHPC), Agency for Science, Technology and Research (A*STAR), Singapore 138632, Singapore
  • 3. China Academy of Launch Vehicle Technology, Beijing Institute of Astronautical Systems Engineering, Beijing 100076, China

Received date: 17 Aug 2019

Accepted date: 26 Sep 2019

Published date: 15 Jun 2020

Copyright

2020 The Author(s) 2020. This article is published with open access at link.springer.com and journal.hep.com.cn

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 (β = 0) and semi-implicit approach with β<0.5 impose a strict stability condition of the transient analysis; (ii) the implicit approach (β=1) and semi-implicit approach with β > 0.5 are generally preferred for their unconditional stability; and (iii) Crank–Nicolson type approach (β=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 β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.

Cite this article

Zhen-Pei WANG , Zhifeng XIE , Leong Hien POH . An isogeometric numerical study of partially and fully implicit schemes for transient adjoint shape sensitivity analysis[J]. Frontiers of Mechanical Engineering, 2020 , 15(2) : 279 -293 . DOI: 10.1007/s11465-019-0575-5

Introduction

Designing structures under thermal loadings is an important issue in practical engineering problems [1,2]. 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. [37] with topology optimizations and in Refs. [810] with shape optimizations.
Structural design sensitivity analyses involving transient thermal responses have been studied since the late 1970s with the work of Ref. [11], 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. [12], 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. [13,14], 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. [15,16] with a Crank–Nicolson type semi-implicit approach, in Refs. [17,18] with general implicit approaches, and in Refs. [19,20] with a fully implicit approach. In Refs. [9,10,2126], 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. [27] where fully analytical solution can be obtained, or in Refs. [28,29] where precise time integration scheme is used.
In general, analytical structural shape sensitivity analysis can be complicated (see Refs. [8,30,31]). 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 [32,33]. These features are attractive for the development of an integrated design frame work for curved beams (e.g., in Refs. [3437]), shells (e.g., in Refs. [38,39]) and general curved structures (e.g., in Refs. [4045]). 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. [4652]. A generalized shape optimization method combining level set method and finite cell method [53] for structural designs using IGA can be found in Refs. [54,55]. More information of IGA-based design optimization can be found in a recent review in Ref. [56].
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 [57,58], to give interesting thermal behaviors [59]. 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. [60].

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 Ω s and boundary Γs as shown in Fig. 1. The variable s denotes the design stages that are sequentially modified from the referential/initial design Ω0. The location function of a material point p in Ω s is denoted as x[p, s]. The temperature at location x, time t and design stage s is denoted by θ[x, t;s]. A general design objective functional defined over a time interval T=[0,T] can be characterized as
J:= 0T ς [t]( Ωsψω [θ[x,t;s]]dΩ+ Γsψγ[θ[x,t;s], q[x ,t;s]]d Γ)dt,
in which the time characteristic function ς [t] is defined as
ς [t]={1,ift TςT, 0,otherwise,
and q[x, t;s] is the heat flux on the boundary.
Fig.1 Schematics of initial design at s =0 (left) and updated design at s (right).

Full size|PPT slide

The transient heat conduction problem within a time interval T=[0,T] is governed by
{l [θ[ x,t]]:=ρcθ[x, t]t k 2θ [x,t]Q[x, t]=0 (x, t)Ωs×T ,θ[ x ,t]=θ^[x, t]( x,t) Γ θs×T ,q[ x,t] n[x] =k θ[x, t]n[ x]=q^[x, t]( x,t) Γ qs×T ,q[ x,t] n[x] =k θ[x, t]n[ x]= qe [x,t]=h(θ[ x,t] θe[t]) (x, t)Γes× T, θ[x, 0]= θ0[x] (x, 0) Ω s,
where c>0, ρ>0, and k>0 are the heat capacity, mass density, and thermal conductivity, respectively, Q denotes the body heat generation rate of each volume unit, q=k θ and q =q n are the heat flow inside Ωs and heat flux through Γs, respectively, n represents the unit outward normal vector on the boundary, h denotes the heat convection coefficient at ambient environment, θe is the ambient temperature, θ0[x] represents the initial temperature field in the domain, and and 2 are the gradient and the Laplacian operators, respectively. The problems considered in this work have an essential boundary condition θ ^ defined over Γθs, a natural boundary condition q^ defined over Γqs and a Robin boundary condition defined over Γes.
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
J¯[s]=J[s]+ 0Tl [θ],ϑ Ωsd t,
in which ϑ =ϑ[ x,t] is the adjoint temperature field with the equality constraint nested in
l[θ ],ϑΩs = Ωs(ρcθtϑ+kθϑ Qϑ)d Ω Γs qϑdΓ=0.
With the constraint satisfied, we have J=J¯. The transient adjoint temperature field can be solved with the following adjoint problem
{ρ c ϑ τ k 2ϑ Q *=0, Q*=ςψω,θ (x,τ)Ωs× T, ϑ=ϑ^, ϑ^=ςψ γ,q (x,τ)Γθs×T, q * n=q ^ *, q^*=ς ψγ,θ( x,τ) Γ qs×T ,q *n= q e*=h(ϑ ϑ e) , ϑe=ς ψγ,θh+ς ψγ ,q(x, τ) Γes× T, ϑ[x, 0]=0  ( x,0)|t =T Ω s,
in which τ=Tt is adjoint time, Q* is adjoint volumetric heat supply, ϑ^ and q^* are adjoint Dirichlet and Neumann boundary conditions, respectively, and ϑe is the ambient temperature. Eventually, the shape sensitivity with respect to parameter s can be derived as
dJds= 0T( Ωs Q ϑ dΩ + Γqs (ς ψγ,qϑ) (q^+( q^n)ν n)dΓ)dt+0 T ( Γθs(ς ψγ,θ q*)θ^ d Γ Γes h(ς ψγ,qϑ)(θ n)νndΓ)dtdt +0 T Γs(ς ψωρc θtϑ+kθϑ Qϑ)νndΓ dt+ 0T Γs (ς( ψγ,θ( θn) ψγκ)q ϑn+ (qϑ)κ)νndΓdt.
For problems with no design-dependent boundary conditions and unit volume heating source, the sensitivity presented in Eq. (7) may be expressed simply as [16,61]
dJds| θ^ , q^, q^n, Q= 0= Γs (0 Tgdt ) νd Γ,
where
g =(ς ψω+ρc θtϑ+kθϑ Qϑ γ eh(ςψγ, qϑ)(θ n )+ς(ψ γ,θ( θn)ψγκ) θ t qϑ n+(q ϑ)κ )n,
with
γe={ 1, x Γe, 0,otherwise.
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 [62,63].

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 RI, i.e.,
θ = I RIθ I=Rθ,
with θ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
C θt+ Kθ= f,
in which C and K are capacitance and conductance matrices, respectively, θ is the vector of unknown nodal temperature, and 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 β[0,1], the temperature discretization of Eq. (12) in the time domain can be written as
Cθt +Δt θtΔt+K(β θt+Δt+ (1 β) θt)=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
θt+Δt θ tΔt=(1β) θ˙t +β θ˙t+Δt .
Using Taylor series expansion, we have
θt+Δt =θ t+Δt θ˙t + Δt22 θ¨t+ O[Δt3],
and
θ˙ t+Δ t=θ˙t+Δtθ¨ t+O[Δt2].
From Eq. (15), it can be obtained that
θt+Δt θ tΔt= θ˙t + Δt 2θ¨ t+O[Δt2].
The truncation error can thus be obtained from the difference between Eqs. (16) and (17) as
θt+Δt θ tΔt((1β)θ˙ t+β θ˙ t+Δ t)=( 12 β) Δt θ ¨t+O [Δ t2].
This indicates that (i) the numerical error of the temperature rate θ˙ is in the order of Δt when β0.5 and Δt2 when β=0.5; and (ii) as β 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 ( β=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 θ[x,t] in terms of an orthonormal basis defined in space, τi[x] and the corresponding temporal basis functions φi[t], such that,
θ[x, t]= i=1 nφi[t] τi[x],
in which n corresponds to the dimension of matrices K and C.
Bearing in mind that the orthonormal basis τi[x] satisfy
ij, τiTKτ j= τi T C τj=0 ,
and substituting Eq. (19) to Eq. (12), we obtain
τjTCτ j φ˙[t]+ τjTKτ jφ[ t]= τjT f.
Denoting
cj=τ jTC τj, kj =τ jTK τj, and fj= τjT f,
we have
cj φ ˙j[t] +kjφ j[t ]=f j[t ],
in which cj, kj, and fj[t] 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 θ[x,t] in the form of θ[x,t]=τ[x]eαt with e as the Napier’s constant (e2.71828), and ignoring the load vector (external excitation) f, the system of equations in Eq. (12) becomes
Cθ˙+Kθ=( Kα C)τ eαt =0.
For the global capacitance matrix C, there exists a lower triangular matrix L such that C=LL T. From Eq. (24), we obtain the following eigen problem
( K¯α I ) LT τ=0, with K¯= L1K LT.
The n×n matrix K¯ has n real non-negative eigenvalues ( α1, α2,...,αn) that correspond to n eigenvectors τ¯i=L T τk. One can easily derive that K τi =α iC τi. Referring to Eq. (22), it can be obtained that
aj= kjcj.
Eventually, Eq. (23) can be re-written as
φ˙j[t]+α jφj[ t]=fj[t] cj.
Substituting the time differentiating scheme in Eq. (13) into Eq. (27),
φjt+Δt φjt+αjΔt( βφjt+Δt+(1β)φjt)=Δtf j cj,
which can be rearranged to give
φjt+Δt= 1(1β)αjΔ t1+βαjΔtφ jt+Δt1+β αjΔt f jcj.
To ensure that φjt+Δt remains bounded in time, the recurrence factor should satisfy
| 1 (1β) αjΔt1+βα jΔ t|1,j.
Note that the eigenvalues of a positive matrix are positive. The stability condition can be obtained as
(1 2β)α jΔ t2,j,
which indicates that the recurrence in Eq. (29) is unconditionally stable if 0.5β 1, and conditionally stable if 0β<0.5. For the conditionally stable cases, the critical time-step for stability is given by
Δtsta = 2(1 2β)α¯,
in which α¯=min( αj ), 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., 1(1β) αjΔt1+βα jΔ t< 0, 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 αj is 1(1β )α jΔ t>0, i.e.,
ΔtΔtosc = 1(1 β) αj.
In literature such as Refs. [64,65], it is suggested to use a time-step size smaller than 1(1β)α¯ with α¯ αj, 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 α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 α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. [66]) and it is recommended to choose a time-step size using
Δt Δx2ρc k,
with Δ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 Γ2 for the same amount material used. The optimal solution for the minimum boundary problem is to have a circular outer boundary Γ2. The problem can be expressed as
J= 0T Γ2 h(θ θ e) dΓdt,
with the volume fixed.
Fig.2 The initial plate design and the NURBS parameterization (values in m) [16,20]. Problem parameters are: θ0[x] =100°C, xΩs, θe =0°C, ρ =7800 kg/m3, c =420 J/(kg·°C), k =20 W/(m·°C) and h =50 W/(m2·°C).

Full size|PPT slide

The design time is chosen to be T =300 s. For simplicity, only a quarter of the plate is parametrized, as depicted in Fig. 2(b), using NURBS with knot vectors ξ=[000 13 1 223111] and η=[000111], and control points shown in Table 1. The locations of six control points, denoted as CI, I = 1, 2, …, 6, are chosen as design variables. Due to symmetry, control point‚ C1 is restricted to move only horizontally while C6 only vertically.

Stability and oscillations in the transient analysis

The analysis model is refined with standard k-refinement approach using knot vectors [ 1 20,220,...,1920] and [ 1 10,210,...,910] in the two orthonormal index directions, respectively. This eventually produce an isogeometric model with 288 control points, leading to matrices C and K with a dimension of 288×288. Following Eq. (25), the largest eigenvalue for matrix K¯ is 422.22. The critical time-step size of the stability condition in Eq. (32) versus the time integration scheme coefficient β is plotted in Fig. 3, which shows that the time-step size needs to be very small to ensure the analysis stability for β<0.5. When the time-step size Δt>0.00948 s with β =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 β.

Full size|PPT slide

The 288 eigenvalues of matrix K¯ are depicted in Fig. 4(a) with the corresponding critical time-step sizes for the oscillatory conditions in Eq. (33) of β =0.5 and β =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 β =0.75 is twice of that with β =0.5. When time-step size is Δt=1 s, the non-oscillatory eigenmodes correspond to α <2 for β=0.5 and α <4 for β=0.75. When time-step size is Δt=0.1 s, the non-oscillatory eigenmodes corresponds to α <20 for β=0.5 and α <40 for β=0.75. The temperature histories of initial time-steps at points A and C1 with different time-step sizes are plotted in Fig. 5 for β =0.5 and Fig. 6 for β =0.75, respectively. It can be observed that (i) the oscillations of temperature for β =0.5 are much more noticeable than those for β =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 β =0.5 and 0.75, respectively.

Full size|PPT slide

Fig.5 Temperature oscillations at point A and C1 with different time-step sizes and β=0.5 for the first few iterative steps.

Full size|PPT slide

Fig.6 Temperature oscillations at point A and C1 with different time-step sizes and β=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 [16,20]
I (i, j) Location Weight ‚I (i, j) Location Weight
x1I x2I x1I x2I
(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 Δt=0.01 s, the sensitivity analysis, termed GfI for design control point xI, is computed using finite difference (FD) method for β=0.5, 0.75, and 1, respectively:
GfI:=J[ xI+δx] J[ xI ]δ x,
where δx is the perturbation of the locations of design control points. The FD computation for β <0.5 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:
G¯fI= 13( G f1I+ Gf2I+G f3I),
where Gf1I, Gf2I, and Gf3I are the FD gradient for β=0.5, 0.75, and 1, respectively. The relative difference of Gfi comparing to the referential sensitivity analysis is calculated using
Dfi= Gfi G¯fmax G¯f.
The L2 norm of Dfi shows the magnitude of the difference between calculated sensitivity and the referential sensitivity. The L2 norms of the relative difference for Gf1, Gf2, and Gf3 are 1.2924× 10 5, 1.6588× 105, and 2.5184× 10 5, respectively, which provides a quantification of the close match between the three cases.
Tab.2 Sensitivity analysis using FD with Δt=0.01 s for different β and the referential sensitivity of the minimum boundary problem
CI Component FD Referential FD
β =0.5 β =0.75 β =1
C1 1 135874.5067 135876.4566 135887.5961 135879.5198
C2 1 243335.5767 243333.5103 243350.4778 243339.8549
2 –37170.7865 –37162.2082 –37166.9958 –37166.6635
C3 1 897323.9019 897326.7759 897308.0330 897319.5703
2 306189.0384 306194.4444 306182.1844 306188.5557
C4 1 306189.1548 306194.2916 306176.7857 306186.7440
2 897323.3635 897312.5441 897306.3377 897314.0818
C5 1 –37165.4896 –37167.1631 –37163.8816 –37165.5115
2 243336.8791 243328.4244 243347.0800 243337.4611
C6 2 135881.9864 135871.1306 135886.7594 135879.9588

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

For β=0.5, 0.75, and 1, the adjoint shape sensitivity is computed for time-step sizes Δ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 L2 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 β =0.75 has an overall better performance than other two cases. The larger error associated with β =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 L2 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 β<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 β at a given time interval discretization

To further investigate the numerical error of adjoint shape sensitivity with respect to β, the adjoint shape sensitivity is computed for Δt=10, 3, and 0.3 s (corresponding to 30, 100, and 1000 time-steps), respectively, with different values of β ranging from 0.5 to 1. The numerical error versus coefficient β 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 β increases from 0.52 to 1. When β =0.5, the numerical errors for the two coarser time discretizations are relatively big due to the oscillatory solutions. A larger value of β can help to reduce the error of the shape sensitivity analysis caused by the oscillatory solutions.
Fig.8 The L2 norm of the relative difference of the adjoint sensitivity analysis versus β 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 [16,20,67]. The model is parametrized with knot arrays ξ=[0,0,0,0.2,0.4,0.5,0.7,1,1,1] and η=[0, 0 ,0,1,1,1], and control points listed in Table 4. The heat convection coefficient are h1=3.15 ×104 W/(mm2·°C) on Γ1 and h3=2.88 ×104 W/(mm2·°C) on Γ3. Other related parameters are k =27.52× 10 3 W/(mm2·°C), ρ c=2.288 ×103 J/(mm3·°C) and T =500 s.
Fig.9 NURBS parameterization of the initial plunger model (values in mm) [16,20].

Full size|PPT slide

Tab.4 Initial locations of the design control points for the plunger design problem [16,20]
I (i,j) Location Weight ‚I (i,j) Location Weight
x1I x2I x1I x2I
(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 Γ3, which contact the molten glass, is assumed to be θe3=1000°C, while the ambient temperature of boundary Γ1, which contact the cooling fluid, is assumed to be θe1=0°C. Temperature difference along the fixed boundary Γ3 affects the quality of the television bulb, which makes it necessary to design boundary Γ 1 such that the temperature difference along boundary Γ3 can be minimized. For this problem, an objective functional is introduced as
J= 0T Γ3 (θ θ ˜) 2dΓdt,
where the average temperature along Γ3, θ˜[t], is computed using
θ˜=1|Γ3| Γ3 θdΓ.

Stability and oscillations in the transient analysis

The analysis model is refined with standard k-refinement approach using knot vectors [ 1 50,250,...,4950] and [ 1 8,28,...,78] in the two orthonormal index directions, respectively. This eventually produces an isogeometric model with 520 control points and matrices C and K with a dimension of 520×520. Following Eq. (25), the maximum eigenvalue for matrix K¯ is 62.67. The critical time-step size of the stability condition in Eq. (32) versus the time integration scheme coefficient β 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 β <0.5. When the time-step size Δt>0.01778 s with β =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 β.

Full size|PPT slide

The 520 eigenvalues of matrix K¯ are depicted in Fig. 11(a) with the corresponding critical time-step sizes for the oscillatory conditions Eq. (32) of β =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 β =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 δ x i=10 6 and time-step size of Δt=0.05 s. The FD computation for β <0.5 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 L2 norm of the relative difference for Gf1, Gf2, and Gf3 are 0.9862× 10 4, 0.2142× 104, and 1.0319× 10 4, respectively, which conform the close match of the three cases and the average gradient G¯f can be used as the referential gradient.
Tab.5 Sensitivity analysis using FD with Δ t=0.05 s for different β and the referential sensitivity of the plunger design problem
CI Component FD Referential FD
β =0.5 β =0.75 β =1
C1 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
C2 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
C3 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
C4 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
C5 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 β

For β= 0.5, 0.75, and 1, the adjoint shape sensitivity is computed for time-step sizes Δ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 L2 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 β =0.5 has an overall better performance than other two cases.
Fig.12 The L2norm 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 β<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 β at a given time interval discretization

To further investigate the numerical error of adjoint shape sensitivity with respect to β, the adjoint shape sensitivity is computed for Δt=10, 1, and 0.1 s (corresponding to 50, 500, and 5000 time-steps), respectively, with different values of β ranging from 0.48 to 1. The numerical error versus coefficient β 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 β increases from 0.5 to 1.
Fig.13 The L2 norm of the relative difference of the adjoint sensitivity analysis versus β 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 β <0.5 are not recommended due to numerical stability concerns; Crank–Nicolson approach with β =0.5 may induce large error because of oscillatory solutions; semi-implicit approaches with β >0.5 are preferred; and fully implicit approach with β =1 has a lower accuracy than the semi-implicit approaches. Hence, a value around of β 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/.
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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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)

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

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

DOI

Outlines

/