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

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.


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. [3][4][5][6][7] with topology optimizations and in Refs. [8][9][10] 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,[21][22][23][24][25][26], 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.
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 in which the time characteristic function &½t is defined as and q½x, t; s is the heat flux on the boundary. The transient heat conduction problem within a time interval T ¼ ½0, T is governed by 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 ¼ kr 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 r and r 2 are the gradient and the Laplacian operators, respectively. The problems Fig. 1 Schematics of initial design at s ¼ 0 (left) and updated design at s (right). considered in this work have an essential boundary condition defined over Γ s , a natural boundary condition q defined over Γ s q and a Robin boundary condition defined over Γ s e . 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 in which # ¼ #½x, t is the adjoint temperature field with the equality constraint nested in With the constraint satisfied, we have J ¼ e J . The transient adjoint temperature field can be solved with the following adjoint problem in which τ ¼ Tt is adjoint time, Q * is adjoint volumetric heat supply,# andq * 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 dJ ds 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] dJ ds j^ # , q#, rq⋅n, where with 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 R I , i.e., 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 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 β 2 ½0, 1, the temperature discretization of Eq. (12) in the time domain can be written as 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 Using Taylor series expansion, we have and From Eq. (15), it can be obtained that The truncation error can thus be obtained from the difference between Eqs. (16) and (17) as This indicates that (i) the numerical error of the temperature rate _ is in the order of Δt when β≠0:5 and Δt 2 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 f i ½t, such that, in which n corresponds to the dimension of matrices K and C.
Bearing in mind that the orthonormal basis τ i ½x satisfy and substituting Eq. (19) to Eq. (12), we obtain Denoting we have in which c j , k j , and f j ½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 ¼ τ½xe -αt with e as the Napier's constant (e % 2.71828), and ignoring the load vector (external excitation) f , the system of equations in Eq. (12) becomes 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 The n Â n matrix K has n real non-negative eigenvalues (α 1 , α 2 , :::, α n ) that correspond to n eigenvectors Referring to Eq. (22), it can be obtained that Eventually, Eq. (23) can be re-written as Substituting the time differentiating scheme in Eq. (13) into Eq. (27), which can be rearranged to give To ensure that f tþΔt j remains bounded in time, the recurrence factor should satisfy Note that the eigenvalues of a positive matrix are positive. The stability condition can be obtained as 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 in which α ¼ minðα j Þ, 8j.

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 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., In literature such as Refs. [64,65], it is suggested to use a time-step size smaller than 1 ð1 -βÞα with α³α j , 8 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 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. 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 with the volume fixed.
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 and η ¼ ½0 0 0 1 1 1, 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 1 is restricted to move only horizontally while C 6 only vertically.

Stability and oscillations in the transient analysis
The analysis model is refined with standard k-refinement approach using knot vectors 1 20 , 2 20 , :::, 19 20 ! and 1 10 , 2 10 , :::, 9 10 ! in the two orthonormal index direc-tions, 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.
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 C 1 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.

Referential sensitivity calculation
Using a relatively small time-step Δt = 0.01 s, the  sensitivity analysis, termed G I f for design control point x I , is computed using finite difference (FD) method for β = 0.5, 0.75, and 1, respectively: 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: where G I f 1 , G I f 2 , and G I f 3 are the FD gradient for β = 0.5, 0.75, and 1, respectively. The relative difference of G fi comparing to the referential sensitivity analysis is calculated using The L 2 norm of D fi shows the magnitude of the difference between calculated sensitivity and the referential sensitivity. The L 2 norms of the relative difference for G f 1 , G f 2 , and G f 3 are 1:2924 Â 10 -5 , 1:6588 Â 10 -5 , and 2:5184 Â 10 -5 , respectively, which provides a quantification of the close match between the three cases.  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.
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.

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. 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 h 1 ¼ 3:15 Â 10 -4 W/(mm 2 $°C) on Γ 1 and h 3 ¼ 2:88 Â 10 -4 W/(mm 2 $°C) on Γ 3 . Other related parameters are k ¼ 27:52 Â 10 -3 W/(mm 2 $°C), c ¼ 2:288 Â 10 -3 J/(mm 3 $°C) and T ¼ 500 s.
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 where the average temperature along Γ 3 ,½t, is computed using   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 timestep 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  Fig. 7 The L 2 norm of the relative difference of the adjoint sensitivity analysis versus number of time-steps for the minimum boundary problem.  unbounded, which leads to failure of the sensitivity analysis.
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.

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 L 2 norm of the relative difference for G f 1 , G f 2 , and G f 3 are 0:9862 Â 10 -4 , Fig. 9 NURBS parameterization of the initial plunger model (values in mm) [16,20].  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.
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. 5.2.5 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.

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. Fig. 12 The L 2 norm of the relative difference of the adjoint sensitivity analysis versus number of time-steps for the plunger design problem.  Fig. 13 The L 2 norm of the relative difference of the adjoint sensitivity analysis versus β for the plunger design problem.