|Computer Modeling in Engineering & Sciences|
Isogeometric Boundary Element Method for Two-Dimensional Steady-State Non-Homogeneous Heat Conduction Problem
1School of Architectural Engineering, Huanghuai University, Zhumadian, 463000, China
2College of Intelligent Construction, Wuchang University of Technology, Wuhan, 430223, China
*Corresponding Author: Yanming Xu. Email: email@example.com
Received: 10 November 2021; Accepted: 08 January 2022
Abstract: The isogeometric boundary element technique (IGABEM) is presented in this study for steady-state inhomogeneous heat conduction analysis. The physical unknowns in the boundary integral formulations of the governing equations are discretized using non-uniform rational B-spline (NURBS) basis functions, which are utilized to build the geometry of the structures. To speed up the assessment of NURBS basis functions, the Bézier extraction approach is used. To solve the extra domain integrals, we use a radial integration approach. The numerical examples show the potential of IGABEM for dimension reduction and smooth integration of CAD and numerical analysis.
Keywords: Isogeometric analysis; NURBS; boundary element method; heat conduction; non-homogeneous; radial integration method
Isogeometric analysis(IGA) has garnered considerable attention in engineering research since Hughes’ pioneering study. IGA employs spline functions for geometric expressions in CAD, such as non-uniform rational B-spline(NURBS) as basis functions to approximate unknown physics fields as an alternative to the usual numerical techniques based on Lagrange polynomials. IGA offers high-order consistency and versatile refinement strategies, which are very appealing in numerical simulations. IGA has been effectively implemented in numerous sectors owing to various advantages as described next[2–6]. The Bézier extraction procedure has been used to increase the simulation performance and make it easier to apply IGA. Although the notion of IGA was initially established in the context of the finite element technique, its implementation in engineering practice is limited owing to its reliance on volume parameterization, which is incompatible with the boundary representation of CAD. However, the boundary element method (BEM) requires boundary parameterization and is inherently compatible with CAD models. The isogeometric boundary element method (IGABEM) has been effectively applied to a variety of domains, including linear elasticity[8–11], potential problems[12–14], fracture mechanics[15–18], electromagnetics, acoustics[20–29], and structural optimization[30–35], and others.
In this study, we use the NURBS method to construct a two-dimensional isogeometric model. The Bézier curve, proposed by Pierre Bézier in 1962, paved the path for the application of spline functions for modeling curved surfaces. The application of spline functions in curved surface modeling has been initiated, and NURBS has become the most widely used modeling tool in geometric modeling in the past few decades. NURBS has outstanding advantages for curved surface modeling. It can accurately represent quadratic curved surfaces, thus providing a uniform mathematical representation of regular and free surfaces. It includes a weight factor that affects the shape of curved surfaces, making it easier to control and achieve the required shape.
In this study, we focus on the heat conduction problem. Many researchers have used traditional methods to solve heat conduction problems, such as the singular boundary method for steady-state nonlinear heat conduction problems with temperature-dependent thermal conductivity, the combination of the BEM and radial basis function method (RBF) for the evaluation of domain integrals with boundary-only discretization[37,38], IGABEM based on NURBS for 3D steady-state thermal conduction problems for a medium with non-homogeneous inclusions[39,40], and for solving transient heat conduction problems with heat sources[41–46]. The subtraction of singularity technique is applied to compute the singular integrals in the BEM. Subsequently, we have coupled RBF and IGABEM for problems of steady heat conduction with variable coefficients. RBF is used to transform the domain integral generated by variable thermal conductivity parameters and the source term into a boundary integral. The operation enables full utilization of the advantages of IGABEM with NURBS.
The remainder of this paper is organized as follows. Section 2 describes the NURBS basis functions with B-spline functions. In Section 3, the basic theory of steady-state heat conduction is described, and in Section 4, the specific operations of the domain integration treatment and the formulation of the system equations for the iterative solution are detailed. In Section 5, several numerical examples of geometric models constructed by NURBS are discussed to verify the correctness of the algorithm. Finally, the conclusions of this study are summarized in Section 6.
The principles of NURBS are briefly discussed in this section for completeness. The basis functions of NURBS are defined over a knot vector, which is a monotonic growing real number sequence, represented by , where ui is the i-th knot in vector U. Iteratively, the B-spline basis function is evaluated as Eq. (1):
When , we have the expression for Ni, p, as Eq. (2):
In Eqs. (1) and (2), p denotes the rank of the polynomial in the B-spline basis functions.
The derivative of the basis function can be written as Eq. (3):
The linear combination of NURBS basis functions and the related control points describes the B-spline curve as Eq. (4):
where, are the Cartesian coordinates of control points.
The weight coefficient is used to rationalize the B-spline basis, resulting in non-uniform rational B-spline (NURBS) basis functions as Eq. (5):
where we have the expression for W as shown in Eq. (6):
Similar to B-spline, the NURBS curve is expressed as Eq. (7):
The rational form of the NURBS basis function, which consists of the B-spline basis function plus the weight factor , allows the precise description of conic curves such as ellipses and circles.
3 Integral Equations for Steady-State Inhomogeneous Heat Conduction Problems
For most materials, the thermal conductivity varies with temperature. For example, the thermal conductivity of a certain composite material at room temperature is 10.85 , at 200 K and 400 K, its thermal conductivity could be 17.73 and 20.16 , respectively. If the thermal conductivity is a function of temperature, the control equation is a nonlinear partial differential equation, and we study the BEM for steady-state inhomogeneous heat conduction problems.
The steady-state heat transfer equation for an isotropic medium containing an internal heat source with temperature-dependent thermal conductivity is expressed as Eq. (8):
In Eq. (8), the thermal conductivity is a function of temperature , is a function of the coordinates of point , indicates the value of the heat source at point , and xi is the coordinate component of point . The boundary conditions are shown in Eqs.(9), (10) and (11):
In Eq. (9), indicates a constant temperature boundary; in Eq. (10) denotes a constant temperature flux boundary, and is the direction vector of the unit outer normal at the boundary point. Eq. (11) is a mixed boundary condition. Now, we introduce the weight function and the weighted residual operation and divisional integration operation, which can be expressed as Eq. (12):
Applying the divisional integration method and Gaussian scattering theorem to the first term, we get Eq. (13):
Considering the case where the source point p is on the boundary containing the field point q, we may write Eq. (14):
4 Analytic Expressions for the Conversion of Domain Integrals to Boundary Integrals
Two domain integrals and two boundary integrals appear in the unified equation. The boundary integrals can be calculated directly by discretizing the Gauss-Legendre quadrature. However, the latter two domain integrals are due to the nonlinear thermal conductivity and heat source terms, respectively. For the domain integrals from the nonlinear thermal conductivity, as the integrals contain , interpolation approximation substitution is needed using radial basis functions and primary polynomials. For integration due to the heat source term, we can use the radial integration method to convert the coordinate representation. Therefore, we focus on the domain integrals due to thermal conductivity nonlinearity and the heat source term.
4.1 Treatment of Domain Integration Due to Thermal Conductivity Nonlinearity
We use radial basis functions and primary polynomials to represent k, iT, which is , as shown in Eqs.(15) and (16):
Here, ; and are the coefficients to be determined, is the radial basis function, which is a most commonly used function, as shown in Table 1.
By applying the radial basis function in Eqs.(15) and (16) to all nodes, we obtain a set of algebraic equations written in matrix form as Eq. (17):
where is the column vector consisting of the coefficients. Matrix f is written as Eq. (18):
The matrix is invertible under the condition that no nodes overlap, and this is obtained by the chain rule shown in Eq. (19):
where is the value of the derivative of the unknown temperature at the node, which can be solved by the integral equation method or by performing an augmented basis function expression approximation. We have the expression for T(x), as shown in Eqs.(20) and (21):
The temperature gradient can be expressed as Eq. (22):
in which denotes the derivative of the radial basis function with respect to R, as shown in Table 2. We can write the expressions for R and its partial derivative with respect to xi, as shown in Eqs.(23) and (24):
In Eqs.(20) and (21), the coefficients to be determined can be obtained from the collocation of the boundary and interior points, and the matrix equation can be written as Eq. (25):
In Eq. (25), the column vector of the temperature at each point can be calculated using the value of the previous iteration, and the last term of the domain integral of Eq. (14) is approximated by substitution using the radial integration method. Hence, we may write Eq. (26):
and Eq. (27):
We continue the simplification of the equation of C in Eq. (27), as shown in Eq. (28):
The original equation is thus obtained by simplifying and combining similar terms as Eq. (29):
in which, we have Eq. (30):
The Gaussian product formula is used to solve the radial integral analytically for the inverse complex quadratic radial basis function, which is usually in the form of a tight branch, as shown in Eq. (31):
It can be noted that is invalid for points outside the compactly supported radius. Eq. (31) can be solved radially by substituting into Eq. (30). The computation of radial integrals by numerical methods is time consuming and can be improved by the analytical solution of Eq. (30).
For Eq. (31), we use the radial integration method for the analytical solution, as shown in Eq. (32):
Fig. 1 shows the radial integration distance diagram. For a few special cases, the source point, field point, and action point are in a straight line as shown in Fig. 2. In such cases, we can solve Eq. (32) analytically. In the radial basis division method, the use of numerical integration requires considerable time, while the use of analytical solutions (as in Eq. (30)), can improve the efficiency of the calculation by several times.
4.2 Processing of Domain Integrals Caused by Source Terms
We first consider the domain integral for the heat source, which is the third term on the right hand side of Eq. (14). By using the radial basis function method, the following expression can be obtained because the boundary integral can be calculated directly by Gaussian integration. We only need to consider the treatment of the domain integral. First, we consider the domain integral owing to the heat source, i.e., the second term on the right-hand side of Eq. (13) is treated using the radial basis method, and the following expression can be obtained as Eq. (33):
where is the basis function obtained from the heat source function after the radial basis method and is expressed as Eq. (34):
We have considered a three-dimensional problem here. Therefore, the power of r is 3-1=2. A specific derivation method can be observed in the radial basis partition. For simple heat source distribution functions, the basis function B can be derived analytically, and for complex heat source distribution functions, it can be approximated by Gaussian numerical integration.
4.3 The Set of System Equations
The unified integral boundary is discretized into Ne boundary cells, the number of boundary nodes is Nb, and the domain has Ni points. The total number of nodes is NA = Nb + Ne, NA nodes are considered as the source p in turn, and the domain integral matrix equation can be calculated, as Eq. (35):
Substituting Eq. (25) into Eq. (35), we obtain Eq. (36):
Finally, a unified system of algebraic equations is obtained, as shown in Eq. (37):
indicates that the elements of the matrix are temperature dependent. The system equations are then obtained by imposing the boundary conditions in Eqs.(9) and (10), and we obtain the system of equations as Eq. (38):
In Eq. (38), A denotes the elements in matrix as a function of temperature T, is a column vector consisting of all unknown node temperatures and unknown boundary fluxes, matrix has values only at the node positions corresponding to the thermal radiation boundary conditions, and the elements at other positions are 0.
4.4 Iterative Solution of the System Equations
Solving the system of linear equations using the Newton-Raphson iteration, we assume that after the nth iteration, the residual of Eq. (38) is . By iterating continuously until convergence, we obtain as Eq. (39):
For the iteration of the (n + 1)th time step and assuming the residual to be 0, we use Taylor series expansion to get Eq. (40):
We need to derive Eq. (39), and the formula can be obtained from Eq. (41):
Substituting Eq. (41) into Eq. (40), we obtain Eq. (42):
We can find the correction value of the unknown quantity , which can then be expressed as Eq. (43):
where is the relaxation factor that changes the convergence, and 0 1. The convergence of iterations can be influenced by changing the relaxation factor.
5 Numerical Examples
5.1 2D Plate Model
We consider a two-dimensional model diagram consisting of a square of 1m side length on the left and a semicircle of radius 0.5m on the right side. Additionally, there are adiabatic left and right boundaries, a constant wall temperature of 100K for the upper boundary, and three boundary conditions of 200K for the lower boundary with convective heat transfer and convective radiative heat transfer. The material density is , specific heat capacity cp is , and thermal emissivity h is . The thermal conductivity of the material as a function of temperature is k(T) = 15+0.01T2 , and the boundary element calculation model and the configuration of the interior points are shown in Fig. 3.
The original NURBS curve and control point positions are shown in Fig. 4a. The parameter space vector of the boundary element may be produced by normalizing the knot vector, which is written as intervals of [0, 0.25), [0.25, 0.5), [0.5, 0.75), and [0.75, 1), respectively. The new control point sequence, knot vector, and weight factor vector may be produced after refinement by dividing each NURBS element evenly and inserting new knots.
Fig. 5a shows the temperature distribution under convective boundary conditions. In the figure, near the heat source at high temperature, the temperature from the lower to upper region increases in accordance with the law of heat transfer. Fig. 5b shows the temperature distribution cloud under temperature boundary conditions. Fig. 5c shows the temperature distribution cloud under convective radiation boundary conditions and it is similar to Fig. 5a, and Fig. 6 shows the temperature of the inner point. When the boundary condition is only temperature-based, the magnitude of the temperature value is proportional to the proximity of the heat source, and when the boundary condition is based on convective or convective radiation, the temperature rises more slowly and the difference between the two is insignificant.
Table 3 shows the calculated temperatures of the internal points under different boundary conditions. We provide the calculated temperatures from the IGABEM and FLUENT software on the flat plate boundary. From Table3 and Fig. 6, we can observe that the IGABEM and FLUENT results have an error of less than 0.7%. Hence, the radial integration method has a high accuracy for solving two-dimensional steady-state nonlinear heat conduction problems. In Fig. 6, we refer to the temperature boundary conditions as TBC, CHT is the convective heat transfer, and CRHT is the convective radiative heat transfer.
Now, we consider the temperature boundary condition as an example to study the influence of the number of internal nodes on the accuracy of the calculation results. The number of internal points is considered as 14, 42, 66, and 86 for comparison, as shown in Fig. 7, and the relative error is shown in Fig. 8 and Table 4. It is observed that compared with 86 internal points, the errors of 14 and 42 internal points are larger, with a maximum error of 2.8%, but the calculation result of 66 internal points is very close. Therefore, we can conclude that for smaller number of interior points, there is greater uncertainty in the result, and for larger number of interior points, the temperature value tends to be more stable, which verifies the correctness of the algorithm. Increasing the number of interior points can enhance computation accuracy while simultaneously increasing calculation difficulty. A topic worth investigating is how to choose the number and position of internal points. It’s tough to present the determined internal points utilizing a theoretical formulation for this problem’s examination. However, uniform distribution works effectively in most cases.
5.2 Rectangular and Semi-Elliptical Models
Let us consider a two-dimensional model diagram of rectangular and semi-elliptical models. The upper region is a rectangle of length 2m and width 1m, the lower region is a semi-ellipse and the analytical equation is . The domain has adiabatic left and right boundaries, a constant wall temperature of 150K for the upper boundary, and 100K for the lower boundary. The material density is , and the specific heat capacity cp is . The thermal conductivity of the material as a function of temperature is . The boundary element calculation model is shown in Fig. 9a, and the configuration of the interior points is shown in Fig. 9b. We choose different coefficients of T2, and verify the correctness of the algorithm.
Fig. 9a shows the model of the rectangular and semi-ellipse region. We assign the coordinates (1, 0), (1, 1), ( −1, 1), ( −1, 0), and (0, −0.5) to describe the silhouette of the plate, in order to make the calculation effective. The inner points are shown in Fig. 9b, and the coordinates of the inner points are ( −0.9, 0.1), ( −0.8, 0.1) ... (0.9, 0.1). There are a total of 19 inner points.
Fig. 10a indicates that the thermal conductivity is a constant value . In the diagram, the right side is a constant 150K heat source, the left side is a constant 100K heat source and heat flows from higher to lower temperature values. The left and right temperature trends are basically symmetrical, which conforms to the thermodynamics theorem. As the inhomogeneous coefficient increases, the temperature symmetry changes. For higher temperature value, the thermal conductivity of the materialis higher, which leads to an increase in the heat transfer rate and relatively high temperature values. This can be observed in Fig. 10b and Fig. 10c. In order to better visualize the effect of inhomogeneous thermal conductivity on the variation in temperature values, we present Fig. 11. With constant thermal conductivity, the temperature is lower than that with variable coefficients, and for larger coefficient, the temperature value is larger.
In this study, the isogeometric boundary element technique (IGABEM) is applied to solve the two-dimensional non-homogeneous steady state heat conduction problem, and NURBS is used to construct a smooth geometric model. The radial integration method is used to transform the domain integral caused by the variable coefficient into the boundary integral. The numerical results show that the developed algorithm can effectively solve non-homogeneous steady-state heat conduction problems with variable coefficients. In the future, we will extend the algorithm to transient analysis.
Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.
Funding Statement: This work is supported by Key Scientific Research Projects of Universities and Key Scientific and Technological Projects in Henan Province, which numbers are 21A440015, 22A570007 and 212102310601, respectively.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|