Scaled boundary finite element method with exact defining curves for two-dimensional linear multi-field media

Abstract

This paper presents an efficient and accurate numerical technique based upon the scaled boundary finite element method for the analysis of two-dimensional, linear, second-order, boundary value problems with a domain completely described by a circular defining curve. The scaled boundary finite element formulation is established in a general framework allowing single-field and multi-field problems, bounded and unbounded bodies, distributed body source, and general boundary conditions to be treated in a unified fashion. The conventional polar coordinates together with a properly selected scaling center are utilized to achieve the exact description of the circular defining curve, exact geometry of the domain, and exact spatial differential operators. Standard finite element shape functions are employed in the discretization of both trial and test functions in the circumferential direction and the resulting eigenproblem is solved by a selected efficient algorithm. The computational performance of the implemented procedure is then fully investigated for various scenarios to demonstrate the accuracy in comparison with standard linear elements.

Keywords

multi-field problems / defining curve / exact geometry / general boundary conditions / SBFEM

Introduction

A scaled boundary finite element method (SBFEM) has been well recognized, in the past three decades, as one of accurate and powerful solution procedures for analysis of various problems in applied mechanics, especially when involving domains are unbounded [¹,²]. In this particular technique, the solution is mainly assumed in a separable form in terms of functions described by two sets of spatial coordinates commonly termed the scaling and boundary coordinates. A part of the solution depending on the boundary coordinate is approximated using available standard discretization procedures whereas that depending on the scaling coordinate is determined directly by solving a certain set of ordinary differential equations [¹,³]. Due to the reduction of one spatial dimension in the discretization, degrees of freedom are equivalently located at the domain boundary, therefore, the method can be categorized as a particular type of boundary element methods. It is worth noting, however, that the formulation of the SBFEM contains no fundamental solution and, as a direct consequence, the non-trivial treatment of singular integrals is not required in comparison with boundary integral equation methods [⁴–⁶]. Such positive and attractive features render this particular technique suitable and computationally efficient when applied to various scenarios including problems related to potential fields, linear elasticity, electrostatics, fractures, elasto-plastic, wave propagation, and geo-mechanics applications [¹,³,⁷–¹⁶]. The collection of earlier investigations indicating the history and background, the current advances, and the enhanced capability of the SBFEM is briefly summarized below.

The primitive formulation of the SBFEM was initially proposed by Wolf and Song [²] in the dynamic analysis of linear elastic unbounded media. The complexity of their original formulation based upon the mechanical-based approach was later improved by applying either the standard weighted residual technique [¹⁷,¹⁸] or the well-known principle of virtual work [³,¹⁹]. While the SBFEM, in its early age, was developed specifically to treat special classes of problems involving unbounded media, its capability has been continuously enhanced and is recently far superior to its original version [⁷–¹¹,¹⁶,²⁰–²⁵]. For instance, use of various families of interpolation functions, including both element-based and element-free bases, in the solution discretization and the corresponding computational performance was fully investigated [²¹–²³,²⁶]. Both h- and p-hierarchical adaptive procedures together with various refinement criteria were also integrated into the SBFEM to enhance the solution accuracy with optimal computational costs [¹,²⁴,²⁵]. The capability of the SBFEM to solve linear elasticity problems associated with applied concentrated forces was achieved by directly incorporating the known fundamental solutions into the solution procedure [⁹]. Extension of the technique to investigate linear boundary value problems involving coupled-field materials such as piezoelectric solids has been also successfully developed [⁷,⁸,¹⁶,²⁷,²⁸]. In addition, applications of the SBFEM to solve problems involving embedded static and moving singularities such as the displacement discontinuities and their propagation has been increasingly found and its success and attractiveness stem directly from the fundamental characteristic of the technique to directly capture the order of singularity and behavior of the surrounding fields [⁸,¹⁶,²⁹–³¹]. Besides linear boundary value problems, the SBFEM along with the polygon formulation was also applied successfully in the stress analysis of bodies made of elasto-plastic materials [¹⁰]. It is known that the accuracy of the SBFEM relies mainly upon the discretization of both the geometry and solution. For relatively complex domains, it is generally required sufficiently fine meshes to achieve the accurate representation of the geometry and this is a non-trivial issue in the computer-aided designs. The concept of isogeometric analysis (IGA) has been recently employed to enhance the approximation of the geometry in the SBFEM (e.g., Ref. [³²]). In addition, IGA has also been successfully integrated in the boundary integral equation methods for the analysis of both two- and three-dimensional problems in linear elasticity (e.g., Refs. [³³,³⁴]).

All aforementioned studies have demonstrated the significant progress of the SBFEM in the analysis of numerous engineering problems. However, the key formulation and implementations of the solution procedure were mainly established in a problem-dependent style. Treatment of various types of boundary value problems and all essential ingredients such as body loads, side-faces, and prescribed conditions based upon a unified or single formulation has not been well established. In addition, use of the exact description of the defining curve to further improve the approximation of the domain geometry and its corresponding boundary has not been recognized. Standard C⁰ linear finite elements have been commonly employed in the SBFEM and found to provide reasonably accurate results with relatively low number of degrees of freedom in the discretization when the domain boundary consists of piecewise linear segments. However, such low-order elements can be inefficient, in terms of degrees of freedom required, when applied to problems with circular or curved boundaries. Relatively fine meshes are generally required to accurately capture the domain geometry and the prescribed boundary conditions. In the present study, the SBFEM is formulated in a general framework allowing two-dimensional, linear, second-order, multi-field boundary value problems to be treated in a unified fashion. The exact defining curve for a circular arc is also established and integrated into the current implementation to investigate the quality of approximate solutions in comparison with those generated by simple linear elements. Results from the present study should, at least, provide some useful information and guidelines for choosing either exact or approximate defining curves in the SBFEM. Remaining sections of this paper are organized to cover the problem formulation, scaled boundary formulation, solution procedure, numerical results and discussion, and conclusion and remarks.

Problem formulation

Consider a two-dimensional body occupying a region \mathrm{\Omega }\in {ℝ}^{\mathrm{2}} as shown schematically in Fig. 1. The region is assumed sufficiently smooth such that all involved mathematical operators such as differentiations and integrations are well defined on this region. In addition, the boundary of the body \mathrm{\Omega }, denoted by \partial \mathrm{\Omega }, is assumed piecewise smooth and an outward unit normal vector at any smooth point on \partial \mathrm{\Omega } is denoted by n={\{n_{\mathrm{1}}{n}_{\mathrm{2}}\}}^{\mathrm{T}}. The state of the body \mathrm{\Omega } is assumed to be completely described by the following three basic field quantities: the state variable denoted by a \mathrm{\Lambda }-component vector u=u(x), the state-variable gradient denoted by a \mathrm{2}\mathrm{\Lambda }-component vector \bar{\epsilon }=\bar{\epsilon }(x) and the body flux denoted by a \mathrm{2}\mathrm{\Lambda }-component vector \sigma =\sigma (x) where\mathrm{\Lambda }is any positive integer. The three field quantities u(x), \bar{\epsilon }(x) and \sigma (x) are governed by three linear field equations including the fundamental law of conservation, the constitutive law of materials, and the relation between the state variable and its measure of variation:

(1)L^T\sigma +b=\mathrm{0},

(2)\sigma =D\overline{\epsilon },

(3)\overline{\epsilon }=Lu,

where L represents the linear differential operator defined by

(4)\mathit{L}={\mathit{L}}_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}+{\mathit{L}}_{\mathrm{2}}{\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}},{\mathit{L}}_{\mathrm{1}}=\left[\begin{array}{l}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ {\mathbf{0}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],L_{\mathrm{2}}=\left[\begin{array}{l}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ {\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],

with \partial /\partial x_{\mathrm{1}} and \partial /\partial x_{\mathrm{2}} denoting the partial derivatives with respect to the coordinates x_{\mathrm{1}} and x_{\mathrm{2}}, respectively, and I_{\mathrm{\Lambda }\times \mathrm{\Lambda }} and {\mathrm{0}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }} standing for \mathrm{\Lambda }\times \mathrm{\Lambda } identity- and zero-matrices, respectively, the superscript “T” denotes the matrix transpose operator, b is a \mathrm{\Lambda }-component vector denoting the prescribed distributed body source, and D is a \mathrm{2}\mathrm{\Lambda }\times \mathrm{2}\mathrm{\Lambda }modulus matrix containing prescribed constants used for completely characterizing the behaviour of a constituting material of the body. At any smooth point x\in \partial \mathrm{\Omega }, the surface flux, denoted by a \mathrm{\Lambda }-component vector t=t(x), is linearly related to the body flux \sigma (x) by

(5)\mathit{t}=\left[\begin{array}{ll}{n}_{\mathrm{1}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}& n_{\mathrm{2}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right]\sigma \mathrm{.}

Based on the prescribed information on the boundary of the body, \partial \Omega can be decomposed into a portion \partial {\Omega }^u on which the state variable is fully prescribed (i.e., u=u^{\mathrm{0}}(x) where u^{\mathrm{0}} is a given function) and a (complementary) portion \partial {\Omega }^t on which the surface flux is prescribed (i.e., t=t^{\mathrm{0}}(x) where t^{\mathrm{0}} is a given function). It should be clear from the form of the basic field Eqs. (1)–(3) and the body-surface-flux relation Eq. (5) that various types of boundary value problems can be readily treated as a special case by appropriately setting the value of the integer \mathrm{\Lambda }and defining the modulus matrix D, for instance, Laplace’s and Poisson’s equation (\Lambda =\mathrm{1}), linear elasticity (\Lambda =\mathrm{2}), linear piezoelectric and linear piezomagnetic problems (\Lambda =\mathrm{3}), and linear piezoelectromagnetic problems (\Lambda =\mathrm{4}) (see also the work of Chung [²⁷]).

The weak formulation of the above governing field equations can be readily established via the following standard procedure: The standard weighted residual technique is first applied to the law of conservation Eq. (1); the term involving the operator L is then integrated by parts via Gauss-divergence theorem; and finally, the constitutive law Eq. (2) and the relation Eq. (3) are directly employed upon substitution. The final weak-form equation in terms of the state variable is given by

(6){\displaystyle \underset{\mathrm{d}}{\int }{(\mathit{L}\mathit{w})}^{\mathrm{T}}\mathit{D}(\mathit{L}\mathit{u})}\mathrm{d}A={\displaystyle \underset{\partial \mathrm{\Omega }}{\int }{\mathit{w}}^{\mathrm{T}}\mathit{t}}\mathrm{d}l+{\displaystyle \underset{\mathrm{d}}{\int }{\mathit{w}}^{\mathrm{T}}\mathit{b}}\mathrm{d}A,

where w is a \Lambda-component vector of sufficiently smooth weight or test functions. To ensure the integrability of all involved integrals in Eq. (6), the test functions w must satisfy the following condition

(7)

Scaled boundary formulation

To introduce the scaled boundary coordinate transformation in {ℝ}^{\mathrm{2}}, a pair (x_{\mathrm{0}},C) is first selected where x_{\mathrm{0}}=(x_{\mathrm{10}},x_{\mathrm{20}})\in {ℝ}^{\mathrm{2}} is termed a scaling center and C\in {ℝ}^{\mathrm{2}} is termed a defining curve. The defining curve can be parameterized by (x_{\mathrm{10}}+{\widehat{x}}_{\mathrm{1}}(s),x_{\mathrm{20}}+{\widehat{x}}_{\mathrm{2}}(s)), s\in [s_{\mathrm{1}},s_{\mathrm{2}}] where s is termed the boundary coordinate and {\widehat{x}}_{\mathrm{1}} and {\widehat{x}}_{\mathrm{2}} are prescribed functions of s. Any point x=(x_{\mathrm{1}},x_{\mathrm{2}})\in {ℝ}^{\mathrm{2}} can be now related to the scaling center x_{\mathrm{0}} and the defining curve C via the following coordinate transformation

(8)x_{\alpha }=x_{\alpha \mathrm{0}}+\xi {\widehat{x}}_{\alpha }(s),\alpha =\mathrm{1},\mathrm{2.}

A choice of the scaling center and the defining curve must ensure that the transformation Eq. (8) is one-to-one. By applying the coordinate transformation Eq. (8), the first-order linear differential operator L given by Eq. (4) can be expressed in terms of partial derivatives with respect to the scaled boundary coordinates \xi and s by

(9)\mathit{L}={\mathit{b}}_{\mathrm{1}}{\displaystyle \frac{\partial }{\partial \xi }}+{\mathit{b}}_{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\xi }}{\displaystyle \frac{\partial }{\partial s}},

where b₁ and b₂ are\mathrm{2}\mathrm{\Lambda }\times \mathrm{\Lambda }-matrices defined by

(10){\mathit{b}}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{J}}\left[\begin{array}{l}{\displaystyle \frac{d{\widehat{x}}_{\mathrm{2}}}{ds}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ -{\displaystyle \frac{d{\widehat{x}}_{\mathrm{1}}}{ds}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],b_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{J}}\left[\begin{array}{l}-{\widehat{x}}_{\mathrm{2}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ {\widehat{x}}_{\mathrm{1}}{\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],

with

J={\widehat{x}}_{\mathrm{1}}\left(d{\widehat{x}}_{\mathrm{2}}\mathrm{/d}s\right)-{\widehat{x}}_{\mathrm{2}}\left(d{\widehat{x}}_{\mathrm{1}}\mathrm{/d}s\right)

. To provide the complete description of the body \mathrm{\Omega } in terms of the scaled boundary coordinates \xi and s, pairs of scaling center and defining curves must be properly chosen and such choices are generally not unique and depends primarily on the geometry of the body itself. In the present study, we focus mainly on the body \mathrm{\Omega } that can be fully described by a single pair of (x_{\mathrm{0}},C). The treatment of bodies with multiple pairs of (x_{\mathrm{0}},C) can be readily extended using the sub-domain technique (e.g., Refs. [³⁵,³⁶]). Now, let any point x\in \mathrm{\Omega } be fully described by the coordinate transformation Eq. (8) for (\xi ,s)\in [{\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}}]\times [s_{\mathrm{1}},s_{\mathrm{2}}]. Portions of the boundaries associated with \xi ={\xi }_{\mathrm{1}} and \xi ={\xi }_{\mathrm{2}} are termed the inner and outer boundaries, respectively. If the defining curve C is not a closed curve (i.e., \widehat{x}(s_{\mathrm{1}})\ne \widehat{x}(s_{\mathrm{2}})), the body is said to be opened and portions of its boundary associated with s=s_{\mathrm{1}} and s=s_{\mathrm{2}} are termed the side-faces; otherwise, the body is closed and contains no side-face. If {\xi }_{\mathrm{2}} is finite, the body is bounded; otherwise, it is unbounded (i.e., {\xi }_{\mathrm{2}}\to \infty). The body contains the scaling center if and only if {\xi }_{\mathrm{1}}=\mathrm{0}.

Solution discretization

The defining curve C is first discretized into n finite elements and m nodes. The state variable \mathit{u} and the test function w can be approximated by

(11){\mathit{u}}^h={\mathit{u}}^h(\xi ,s)={\displaystyle \sum _{i=\mathrm{1}}^m{\varphi }_{(i)}(s){\mathit{u}}_{(i)}^h(\xi )}=\mathit{N}(s){\mathit{U}}^h(\xi ),

(12){\mathit{w}}^h={\mathit{w}}^h(\xi ,s)={\displaystyle \sum _{i=\mathrm{1}}^m{\varphi }_{(i)}(s){\mathit{w}}_{(i)}^h(\xi )}=\mathit{N}(s){\mathit{W}}^h(\xi ),

where the superscript “h” is employed, here and in what follows, to designate approximate quantities, {\phi }_{(i)}(s) denotes the ith nodal basis function, u_{(i)}^h(\xi ) denotes the ith unknown nodal function of the state variable, w_{(i)}^h(\xi ) denotes the ith arbitrary nodal function, N(s) denotes a matrix containing {\phi }_{(i)}(s), and U^h(\xi ) and W^h(\xi ) denotes vector containing nodal functions u_{(i)}^h(\xi ) and w_{(i)}^h(\xi ), respectively. By employing Eqs. (2), (3), (9), (11) and (12), the body flux \sigma and \mathit{L}\mathit{w} can be approximated by

(13){\mathit{\sigma }}^h={\mathit{\sigma }}^h(\xi ,s)=\mathit{D}({\mathit{L}}^h{\mathit{u}}^h)=\mathit{D}\left({\mathit{B}}_{\mathrm{1}}{\mathit{U}}_{,\xi }^h+{\displaystyle \frac{\mathrm{1}}{\xi }}{\mathit{B}}_{\mathrm{2}}{\mathit{U}}^h\right),

(14){\mathit{L}}^h{\mathit{L}}^h={\mathit{L}}^h{\mathit{w}}^h(\xi ,s)=\left({\mathit{B}}_{\mathrm{1}}{\mathit{W}}_{,\xi }^h+{\displaystyle \frac{\mathrm{1}}{\xi }}{\mathit{B}}_{\mathrm{2}}{\mathit{W}}^h\right),

where {(\cdot )}_{,\xi } denotes the derivative with respect to the scaling coordinate \xi and the matrices B₁ and B₂ are given by

(15){\mathit{B}}_{\mathrm{1}}={\mathit{b}}_{\mathrm{1}}^h\mathit{N,}

(16){\mathit{B}}_{\mathrm{2}}={\mathit{b}}_{\mathrm{2}}^h\mathit{B}={\mathit{b}}_{\mathrm{2}}^h\mathrm{d}\mathit{N}/\mathrm{d}\mathit{s}\mathrm{.}

Scaled boundary finite element equations

To establish a set of scaled boundary finite element equations for a general, two-dimensional body \mathrm{\Omega } as shown schematically in Fig. 2, let us assume that the boundary of the domain \partial \mathrm{\Omega } consists of four parts resulting from the scaled boundary coordinate transformation with the scaling center x₀, the defining curve C, and (\xi ,s)\in [{\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}}]\times [s_{\mathrm{1}},s_{\mathrm{2}}]: The inner boundary \partial {\mathrm{\Omega }}_{\mathrm{1}} (associated with \xi ={\xi }_{\mathrm{1}}), the outer boundary \partial {\mathrm{\Omega }}_{\mathrm{2}} (associated with \xi ={\xi }_{\mathrm{2}}), the side-face-1 \partial {\mathrm{\Omega }}_{s\mathrm{1}} (associated with s=s_{\mathrm{1}}), and the side-face-2 \partial {\mathrm{\Omega }}_{s\mathrm{2}} (associated with s=s_{\mathrm{2}}). As a result of such boundary partition \partial \mathrm{\Omega }=\partial {\mathrm{\Omega }}_{\mathrm{1}}\cup \partial {\mathrm{\Omega }}_{\mathrm{2}}\cup \partial {\mathrm{\Omega }}_{s\mathrm{1}}\cup \partial {\mathrm{\Omega }}_{s\mathrm{2}} and the coordinate transformation Eq. (8), the weak-form Eq. (6) becomes

(17)\begin{array}{l}{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{(Lw)}^{T}D(Lu)}}J\xi d\xi ds={\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{w}_{\mathrm{1}}^T{t}_{\mathrm{1}}}(s)J_s(s){\xi }_{\mathrm{1}}ds+{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{w}_{\mathrm{2}}^T{t}_{\mathrm{2}}(s)}{J}_s(s){\xi }_{\mathrm{2}}ds+{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{(w_{s\mathrm{1}})}^T{t}_{s\mathrm{1}}(\xi )}{J}_{\xi \mathrm{1}}d\xi \\ \begin{array}{llll}\begin{array}{lll}\begin{array}{lll}& & \end{array}& & \end{array}& & & \end{array}+{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{(w_{s\mathrm{2}})}^T{t}_{s\mathrm{2}}(\xi )}{J}_{\xi \mathrm{2}}d\xi +{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{w}^{T}b}}J\xi d\xi ds\end{array}

where {\mathit{w}}_{\mathrm{1}}, {\mathit{w}}_{\mathrm{2}},{\mathit{w}}_{s\mathrm{1}},{\mathit{w}}_{s\mathrm{2}} are restrictions of the test function on the boundaries \partial {\mathrm{\Omega }}_{\mathrm{1}}, \partial {\mathrm{\Omega }}_{\mathrm{2}}, \partial {\mathrm{\Omega }}_{s\mathrm{1}}, \partial {\mathrm{\Omega }}_{s\mathrm{2}}, respectively, t_{\mathrm{1}}, t_{\mathrm{2}}, t_{s\mathrm{1}}, t_{s\mathrm{2}} are surface fluxes on the boundaries \partial {\mathrm{\Omega }}_{\mathrm{1}}, \partial {\mathrm{\Omega }}_{\mathrm{2}}, \partial {\mathrm{\Omega }}_{s\mathrm{1}}, \partial {\mathrm{\Omega }}_{s\mathrm{2}}, respectively, and J_{\xi \mathrm{1}}=\sqrt{{\widehat{x}}_{\mathrm{1}}^{\mathrm{2}}(s_{\mathrm{1}})+{\widehat{x}}_{\mathrm{2}}^{\mathrm{2}}(s_{\mathrm{1}})}, J_{\xi \mathrm{2}}=\sqrt{{\widehat{x}}_{\mathrm{1}}^{\mathrm{2}}(s_{\mathrm{2}})+{\widehat{x}}_{\mathrm{2}}^{\mathrm{2}}(s_{\mathrm{2}})} and J_s(s)=\sqrt{{\left(d{\widehat{x}}_{\mathrm{1}}\mathrm{/d}s\right)}^{\mathrm{2}}+{\left(d{\widehat{x}}_{\mathrm{2}}\mathrm{/d}s\right)}^{\mathrm{2}}}. By substituting the approximations Eqs. (11)–(14) into Eq. (17), then integrating the integral on the left hand side by parts, and finally carrying out all involved matrix algebra, it leads to

(18)$$\begin{gathered} \underset{{\xi }_1}{\overset{{\xi }_2}{\int }}{\left(W^h\right)}^{\mathrm{T}}[-\xi E_0{U}_{,\xi \xi }^h+(E_1-E_1^{\mathrm{T}}-E_0)U_{,\xi \xi }^h+\frac{1}{\xi }{E}_2{U}^h-F^t-\xi F^b]\mathrm{d}\xi \\ +{\left(W^h\left({\xi }_2\right)\right)}^{\mathrm{T}}[{(\xi E_0{U}_{,\xi }^h+E_1^{\mathrm{T}}{U}^h)}_{\xi ={\xi }_2}-P_2] \\ -{\left(W^h\left({\xi }_1\right)\right)}^{\mathrm{T}}[{(\xi E_0{U}_{,\xi }^h+E_1^{\mathrm{T}}{U}^h)}_{\xi ={\xi }_1}+P_1]=0, \end{gathered}$$

where the matrices E₀, E₁, E₂, and all involved quantities are defined by

(19)\begin{array}{lll}{\mathit{E}}_{\mathrm{0}}={\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{B}}_{\mathrm{1}}^T\mathit{D}{\mathit{B}}_{\mathrm{1}}}{J}^h\mathrm{d}s,& {\mathit{E}}_{\mathrm{1}}={\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{B}}_{\mathrm{2}}^T\mathit{D}{\mathit{B}}_{\mathrm{1}}}{J}^h\mathrm{d}s,& {\mathit{E}}_{\mathrm{2}}={\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{B}}_{\mathrm{2}}^T\mathit{D}{\mathit{B}}_{\mathrm{2}}}{J}^h\mathrm{d}s\end{array},

(20)\begin{array}{ll}{\mathit{P}}_{\mathrm{1}}={\xi }_{\mathrm{1}}{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{({\mathit{N}}^S)}^T{\mathit{t}}_{\mathrm{1}}}(s)J_s^h(s)\mathrm{d}s,& {\mathit{P}}_{\mathrm{2}}={\xi }_{\mathrm{2}}{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{({\mathit{N}}^S)}^T{\mathit{t}}_{\mathrm{2}}}(s)J_s^h(s)\mathrm{d}s\end{array},

(21)\begin{array}{ll}{\mathit{F}}_{\mathrm{1}}^t={({\mathit{N}}_{\mathrm{1}}^S)}^T{\mathit{t}}_{s\mathrm{1}}(\xi )J_{\xi \mathrm{1}}^h,& {\mathit{F}}_{\mathrm{2}}^t={({\mathit{N}}_{\mathrm{2}}^S)}^T{\mathit{t}}_{s\mathrm{2}}(\xi )J_{\xi \mathrm{2}}^h\end{array},

(22){\mathit{F}}^b={\mathit{F}}^b(\xi )={\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{N}}^T(s)\mathit{b}(\xi ,s)}{J}^h\mathrm{d}s,

From the arbitrariness of the nodal functions W^h, it can be deduced from Eq. (18) that

(23){\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}{\mathit{U}}_{,\xi \xi }^h+\xi ({\mathit{E}}_{\mathrm{0}}+{\mathit{E}}_{\mathrm{1}}^T-{\mathit{E}}_{\mathrm{1}}){\mathit{U}}_{,\xi }^h-{\mathit{E}}_{\mathrm{2}}{\mathit{U}}^h+\xi {\mathit{F}}^t+{\xi }^{\mathrm{2}}{\mathit{F}}^b=\mathrm{0},\forall \xi \in ({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}}),

(24){\mathit{Q}}^h({\xi }_{\mathrm{1}})=-{\mathit{P}}_{\mathrm{1}},

(25){\mathit{Q}}^h({\xi }_{\mathrm{2}})={\mathit{P}}_{\mathrm{2}},

where the vector {\mathit{Q}}^h={\mathit{Q}}^h(\xi ) known as the nodal internal flux is defined by

(26){\mathit{Q}}^h(\xi )=\xi {\mathit{E}}_{\mathrm{0}}{\mathit{U}}_{,\xi }^h+{\mathit{E}}_{\mathrm{1}}^T{\mathit{U}}^h

Equations (23)–(25) form a set of the so-called scaled boundary finite element equations governing the function {\mathit{U}}^h={\mathit{U}}^h(\xi ). In particular, Eq. (23) is a system of linear, second-order, nonhomogeneous, ordinary differential equations with respect to the scaling coordinate \xi whereas Eqs. (24) and (25) correspond to the boundary conditions prescribed on the inner and outer boundaries of the body.

Treatment of prescribed conditions on side-faces

For bodies containing side-faces, the treatment of prescribed conditions on those boundaries can be achieved by partitioning the vector {\mathit{U}}^h into {\left\{\begin{array}{ll}{\mathit{U}}^{hu}& {\mathit{U}}^{hc}\end{array}\right\}}^{\mathrm{T}} where {\mathit{U}}^{hu} denotes a vector containing only the unknown nodal functions and {\mathit{U}}^{hc}denotes a vector containing known nodal functions associated with the prescribed state variable on the side-faces. By partitioning all other involved quantities in a consistent fashion and then partitioning the system of Eq. (23) and the boundary conditions Eq. (24) and (25), it finally leads to the scaled boundary finite element equations for the unknown nodal functions {\mathit{U}}^{hu}:

(27){\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}^{uu}{\mathit{U}}_{,\xi \xi }^{hu}+\xi [{\mathit{E}}_{\mathrm{0}}^{uu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^T-{\mathit{E}}_{\mathrm{1}}^{uu}]{\mathit{U}}_{,\xi }^{hu}-{\mathit{E}}_{\mathrm{2}}^{uu}{\mathit{U}}^{hu}=-\xi {\mathit{F}}^{tu}-{\xi }^{\mathrm{2}}{\mathit{F}}^{bu}-{\mathit{F}}^{suu},

(28){\mathit{Q}}^{hu}({\xi }_{\mathrm{1}})=-{\mathit{P}}_{\mathrm{1}}^u,

(29){\mathit{Q}}^{hu}({\xi }_{\mathrm{2}})={\mathit{P}}_{\mathrm{2}}^u,

where the function F^{suu} is defined by

(30){\mathit{F}}^{suu}={\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}^{uc}{\mathit{U}}_{,\xi \xi }^{hc}+\xi [{\mathit{E}}_{\mathrm{0}}^{uc}+{({\mathit{E}}_{\mathrm{1}}^{cu})}^T-{\mathit{E}}_{\mathrm{1}}^{uc}]{\mathit{U}}_{,\xi }^{hc}-{\mathit{E}}_{\mathrm{2}}^{uc}{\mathit{U}}^{hc},

and \{{\mathit{E}}_{\mathrm{0}}^{uu},{\mathit{E}}_{\mathrm{0}}^{uc}\}, \{{\mathit{E}}_{\mathrm{1}}^{uu},{\mathit{E}}_{\mathrm{1}}^{uc}\}, \{{\mathit{E}}_{\mathrm{2}}^{uu},{\mathit{E}}_{\mathrm{2}}^{uc}\}, {\mathit{F}}^{tu}, {\mathit{F}}^{bu}, {\mathit{Q}}^{hu}, {\mathit{P}}_{\mathrm{1}}^u, {\mathit{P}}_{\mathrm{2}}^u result from the partition of {\mathit{E}}_{\mathrm{0}}, {\mathit{E}}_{\mathrm{1}}, {\mathit{E}}_{\mathrm{2}}, {\mathit{F}}^t, {\mathit{F}}^b, {\mathit{Q}}^h, {\mathit{P}}_{\mathrm{1}}, {\mathit{P}}_{\mathrm{2}}, respectively (also see the work of Chung [²⁷]). It should be remarked that the right-hand side of Eq. (27) now involves only the prescribed conditions on the side-faces and the distributed body source.

Solution procedure

A homogeneous solution of Eq. (27), denoted by {\mathit{U}}_{\mathrm{0}}^{hu}, can be readily obtained by solving a system of linear, second-order, Euler-Cauchy differential equations:

(31){\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}^{uu}{\mathit{U}}_{\mathrm{0},\xi \xi }^{hu}+\xi [{\mathit{E}}_{\mathrm{0}}^{uu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^T-{\mathit{E}}_{\mathrm{1}}^{uu}]{\mathit{U}}_{\mathrm{0},\xi }^{hu}-{\mathit{E}}_{\mathrm{2}}^{uu}{\mathit{U}}_{\mathrm{0}}^{hu}=\mathrm{0},

and it takes the following form

(32){\mathit{U}}_{\mathrm{0}}^{hu}(\xi )={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}N}{c}_i{\xi }^{{\lambda }_i}{\psi }_i},

where N is the number of equations in the system Eq. (31), ({\lambda }_i,{\psi }_i) denote a scaling factor and a vector representing the state variable of the ith mode, and c_i is an unknown constant denoting the contribution of the ith mode to the solution. The nodal internal flux generated by {\mathit{U}}_{\mathrm{0}}^{hu}, denoted by {\mathit{Q}}_{\mathrm{0}}^{hu}, is obtained by substituting Eq. (32) into {\mathit{Q}}^{hu} as

(33){\mathit{Q}}_{\mathrm{0}}^{hu}(\xi )={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}N}{c}_i{\xi }^{{\lambda }_i}{\mathit{q}}_i},

where {\mathit{q}}_i=[{\lambda }_i{\mathit{E}}_{\mathrm{0}}^{uu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}]{\psi }_i is termed the ith modal internal flux. A pair ({\lambda }_i,X_i) with {\mathit{X}}_i={\left\{\begin{array}{ll}{\psi }_i& {\mathit{q}}_i\end{array}\right\}}^{\mathrm{T}} can be obtained by solving the corresponding linear eigenvalue problem

(34)\mathit{A}{\mathit{X}}_i={\lambda }_i{\mathit{X}}_i

where the coefficient matrix A is defined by

(35)\mathit{A}=\left[\begin{array}{ll}-{({\mathit{E}}_{\mathrm{0}}^{uu})}^{-\mathrm{1}}{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}& {({\mathit{E}}_{\mathrm{0}}^{uu})}^{-\mathrm{1}}\\ {\mathit{E}}_{\mathrm{2}}^{uu}-{\mathit{E}}_{\mathrm{1}}^{uu}{({\mathit{E}}_{\mathrm{0}}^{uu})}^{-\mathrm{1}}{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}& {\mathit{E}}_{\mathrm{1}}^{uu}{({\mathit{E}}_{\mathrm{0}}^{uu})}^{-\mathrm{1}}\end{array}\right],

For further use in the solution procedure of different types of bodies, the homogeneous solution for both {\mathit{U}}_{\mathrm{0}}^{hu} and {\mathit{Q}}_{\mathrm{0}}^{hu} can be re-expressed as

(36){\mathit{U}}_{\mathrm{0}}^{hu}(\xi )={\mathbf{\Phi }}^{\psi +}{\mathrm{\Pi }}^{+}(\xi ){\mathit{C}}^{+}+{\mathbf{\Phi }}^{\psi -}{\mathrm{\Pi }}^{-}(\xi ){\mathit{C}}^{-},

(37){\mathit{Q}}_{\mathrm{0}}^{hu}(\xi )={\mathbf{\Phi }}^{q+}{\mathrm{\Pi }}^{+}(\xi ){\mathit{C}}^{+}+{\mathrm{\Phi }}^{q-}{\mathrm{\Pi }}^{-}(\xi ){\mathit{C}}^{-},

where the superscripts “+” and “–” are used to designate quantities associated with {\lambda }_ipossessing positive and negative real parts, respectively, {\mathbf{\Phi }}^{\psi +}and {\mathbf{\Phi }}^{\psi -} are matrices with each column containing a vector {\psi }_i, {\mathbf{\Phi }}^{q+} and {\mathbf{\Phi }}^{q-} are matrices with each column containing a vector q_i; {\mathrm{\Pi }}^{+} and {\mathrm{\Pi }}^{-} are diagonal matrices with each diagonal entries containing a function {\xi }^{{\lambda }_i}; and {\mathit{C}}^{+} and {\mathit{C}}^{-} are vectors containing unknown constants representing the contribution of each mode. Note in particular that the order of columns and entries of involved matrices and vectors in Eqs. (36) and (37) follows consistently the eigenvalues{\lambda }_i. For bodies containing the scaling center, it should be apparent that {\mathrm{\Pi }}^{-} becomes infinite at \xi =\mathrm{0}and the vector {\mathit{C}}^{-} must vanish to maintain the boundedness of the solution at the scaling center. Similarly, for unbounded bodies, the vector {\mathit{C}}^{+} must be chosen equal to zero if the solution is expected to be bounded at infinity.

A particular solution of the system Eq. (27) associated with the distributed body source and the prescribed boundary conditions on the side-faces, denoted by {\mathit{U}}_{\mathrm{1}}^{hu}, can be readily established from existing standard procedures such as the method of undetermined coefficients. Once the particular solution {\mathit{U}}_{\mathrm{1}}^{hu} is obtained, the corresponding particular nodal internal flux {\mathit{Q}}_{\mathrm{1}}^{hu} can be directly calculated in the same fashion as {\mathit{Q}}_{\mathrm{0}}^{hu}. Finally, the general solution of Eq. (27) and the corresponding nodal internal flux are given by

(38){\mathit{U}}^{hu}(\xi )={\mathit{U}}_{\mathrm{0}}^{hu}(\xi )+{\mathit{U}}_{\mathrm{1}}^{hu}(\xi )={\mathbf{\Phi }}^{\psi +}{\mathrm{\Pi }}^{+}(\xi ){\mathit{C}}^{+}+{\mathbf{\Phi }}^{\psi -}{\mathrm{\Pi }}^{-}(\xi ){\mathit{C}}^{-}+{\mathit{U}}_{\mathrm{1}}^{hu}(\xi ),

(39){\mathit{Q}}^{hu}(\xi )={\mathit{Q}}_{\mathrm{0}}^{hu}(\xi )+{\mathit{Q}}_{\mathrm{1}}^{hu}(\xi )={\mathbf{\Phi }}^{q+}{\mathrm{\Pi }}^{+}(\xi ){\mathit{C}}^{+}+{\mathbf{\Phi }}^{q-}{\mathrm{\Pi }}^{-}(\xi ){\mathit{C}}^{-}+{\mathit{Q}}_{\mathrm{1}}^{hu}(\xi ),

By enforcing the boundary conditions on both inner and outer boundaries along with Eqs. (28), (29) and (39), the unknown constants {\mathit{C}}^{+} and {\mathit{C}}^{-} can be obtained as

(40)\left\{\begin{array}{l}{\mathit{C}}^{+}\\ {\mathit{C}}^{-}\end{array}\right\}={\left[\begin{array}{ll}{\mathbf{\Phi }}^{q+}{\mathrm{\Pi }}^{+}({\xi }_{\mathrm{1}})& {\mathbf{\Phi }}^{q-}{\mathrm{\Pi }}^{-}({\xi }_{\mathrm{1}})\\ {\mathbf{\Phi }}^{q+}{\mathrm{\Pi }}^{+}({\xi }_{\mathrm{2}})& {\mathbf{\Phi }}^{q-}{\mathrm{\Pi }}^{-}({\xi }_{\mathrm{2}})\end{array}\right]}^{-\mathrm{1}}\left(\left\{\begin{array}{l}-{\mathit{P}}_{\mathrm{1}}^u\\ {\mathit{P}}_{\mathrm{2}}^u\end{array}\right\}-\left\{\begin{array}{l}{\mathit{Q}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{1}})\\ {\mathit{Q}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{2}})\end{array}\right\}\right)\mathrm{.}

By substituting Eq. (40) into Eq. (38) and then evaluating the result at \xi ={\xi }_{\mathrm{1}} and \xi ={\xi }_{\mathrm{2}}, it finally yields

(41)\mathit{K}\left\{\begin{array}{l}{\mathit{U}}^{hu}({\xi }_{\mathrm{1}})\\ {\mathit{U}}^{hu}({\xi }_{\mathrm{2}})\end{array}\right\}=\left\{\begin{array}{l}-{\mathit{P}}_{\mathrm{1}}^u\\ {\mathit{P}}_{\mathrm{2}}^u\end{array}\right\}+\mathit{K}\left\{\begin{array}{l}{\mathit{U}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{1}})\\ {\mathit{U}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{2}})\end{array}\right\}-\left\{\begin{array}{l}{\mathit{Q}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{1}})\\ {\mathit{Q}}_{\mathrm{1}}^{hu}({\xi }_{\mathrm{2}})\end{array}\right\},

where the coefficient matrix K is given by

(42)K=\left[\begin{array}{ll}{\Phi }^{q+}{\Pi }^{+}({\xi }_{\mathrm{1}})& {\Phi }^{q-}{\Pi }^{-}({\xi }_{\mathrm{1}})\\ {\Phi }^{q+}{\Pi }^{+}({\xi }_{\mathrm{2}})& {\Phi }^{q-}{\Pi }^{-}({\xi }_{\mathrm{2}})\end{array}\right]{\left[\begin{array}{ll}{\Phi }^{\psi +}{\Pi }^{+}({\xi }_{\mathrm{1}})& {\Phi }^{\psi -}{\Pi }^{-}({\xi }_{\mathrm{1}})\\ {\Phi }^{\psi +}{\Pi }^{+}({\xi }_{\mathrm{2}})& {\Phi }^{\psi -}{\Pi }^{-}({\xi }_{\mathrm{2}})\end{array}\right]}^{-\mathrm{1}}\mathrm{.}

By applying the prescribed conditions on both inner and outer boundaries to Eq. (41), all involved unknowns on those boundaries can be determined. Once such unknowns are solved, the constants {\mathit{C}}^{+} and {\mathit{C}}^{-} are computed from Eq. (40), the unknown nodal functions contained in {\mathit{U}}^{hu} can be determined from Eq. (38), and the state variable and the body flux at any point within the body can be obtained from Eqs. (11) and (13), respectively. The relative error of the scaled boundary finite element solution is defined by

(43)Error={\displaystyle \frac{\left|\right|\mathit{e}|{|}_{L^{\mathrm{2}}}}{\left|\right|{\mathit{u}}^{\mathrm{exact}}|{|}_{L^{\mathrm{2}}}}}\times \mathrm{100}\%,

where \mathit{e}={\mathit{u}}^{\mathrm{exact}}-{\mathit{u}}^h denotes the error function, {\mathit{u}}^{\mathrm{exact}} is an exact or reference solution, and {\Vert \mathit{f}\Vert }_{L^{\mathrm{2}}} denotes the L_{\mathrm{2}}-norm of a vector-value function \mathit{f} given by

(44){\Vert \mathit{f}\Vert }_{L^{\mathrm{2}}}=\sqrt{{\displaystyle \underset{\Omega }{\int }{f}^{T}fdA}}\mathrm{.}

Description of linear straight and circular arc elements

In this section, the description of two types of elements used in the present study to describe the domain geometry and to discretize the solution in the boundary direction, one associated with a 2-node isoparametric linear element and the other corresponding to a 2-node circular-arc element (Fig. 3), is briefly presented. For the 2-node isoparametric linear element, the shape functions used in the discretization of the defining curve and the solution are given explicitly by

(45){\varphi }_{\mathrm{1}}(s)=(\mathrm{1}-s)/\mathrm{2},{\varphi }_{\mathrm{2}}(s)=(\mathrm{1}+s)/\mathrm{2},

,

where

s\in [-\mathrm{1},\mathrm{1}]

denotes the local or element boundary coordinate. The approximate defining curve of this particular element is parameterized by

(46)x_{\alpha }^h(s)=x_{\alpha \mathrm{0}}+{\widehat{x}}_{\alpha \mathrm{1}}{\varphi }_{\mathrm{1}}(s)+{\widehat{x}}_{\alpha \mathrm{2}}{\varphi }_{\mathrm{2}}(s),\alpha =\mathrm{1},\mathrm{2}

where x_{\mathrm{1}}=(x_{\mathrm{10}}+{\widehat{x}}_{\mathrm{11}},x_{\mathrm{20}}+{\widehat{x}}_{\mathrm{21}}) and x_{\mathrm{2}}=(x_{\mathrm{10}}+{\widehat{x}}_{\mathrm{12}},x_{\mathrm{20}}+{\widehat{x}}_{\mathrm{22}}) are physical coordinates of the two nodes of the element on the defining curve. By using the parameterization Eq. (46), the approximation of the matrices b₁ and b₂ and J in Eq. (10) can be obtained explicitly by

(47){\mathit{b}}_{\mathrm{1}}^h={\displaystyle \frac{\mathrm{1}}{J^h}}\left[\begin{array}{l}{\mathrm{\Delta }}_{\mathrm{2}}{\mathit{I}}_{\mathrm{\Delta }\times \mathrm{\Delta }}\\ -{\mathrm{\Delta }}_{\mathrm{1}}{\mathit{I}}_{\mathrm{\Delta }\times \mathrm{\Delta }}\end{array}\right],{\mathit{b}}_{\mathrm{2}}^h={\displaystyle \frac{\mathrm{1}}{J^h}}\left[\begin{array}{l}-({\Sigma }_{\mathrm{2}}+{\mathrm{\Delta }}_{\mathrm{2}}s){\mathit{I}}_{\mathrm{\Delta }\times \mathrm{\Delta }}\\ ({\Sigma }_{\mathrm{1}}+{\mathrm{\Delta }}_{\mathrm{1}}s){\mathit{I}}_{\mathrm{\Delta }\times \mathrm{\Delta }}\end{array}\right],J^h={\Sigma }_{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{2}}-{\Sigma }_{\mathrm{2}}{\mathrm{\Delta }}_{\mathrm{1}},

where {\Sigma }_{\mathrm{1}}=\left({\widehat{x}}_{\mathrm{11}}+{\widehat{x}}_{\mathrm{12}}\right)/\mathrm{2}, {\Sigma }_{\mathrm{2}}=\left({\widehat{x}}_{\mathrm{21}}+{\widehat{x}}_{\mathrm{22}}\right)/\mathrm{2}, {\mathrm{\Delta }}_{\mathrm{1}}=\left({\widehat{x}}_{\mathrm{12}}-{\widehat{x}}_{\mathrm{11}}\right)/\mathrm{2} and {\mathrm{\Delta }}_{\mathrm{2}}=\left({\widehat{x}}_{\mathrm{22}}-{\widehat{x}}_{\mathrm{21}}\right)/\mathrm{2}. It is apparent from Eq. (47) that both J^h and b_{\mathrm{1}}^h are constant whereas b_{\mathrm{2}}^h is a linear function of s. Using Eq. (47) along with Eqs. (15) and (16) and the fact that the solution on the element is also discretized by the same linear shape functions Eq. (45), all entries of the matrices B₁ and B₂ are linear functions of s and, as a result, the matrices E₀, E₁, E₂ for the 2-node linear element can be readily obtained in a closed form. It is noted that the representation of the geometry of the body and its associated boundaries using the 2-node linear elements can be exact if the defining curve is piecewise linear. However, when the defining curve used in describing the body becomes more complex or non-straight, a significant number of linear elements is required, in general, to achieve the accurate description of the geometry.

For the 2-node circular-arc element, the defining curve on this element is parameterized in terms of the local boundary coordinate s\in [-\mathrm{1},\mathrm{1}] by

(48)x_{\mathrm{1}}=x_{\mathrm{10}}+R\cos \theta (s),x_{\mathrm{2}}=x_{\mathrm{20}}+R\sin \theta (s),

where R denotes the radius of the circular arc, \theta (s)={\theta }_{\mathrm{1}}{\phi }_{\mathrm{1}}(s)+{\theta }_{\mathrm{2}}{\phi }_{\mathrm{2}}(s), {\varphi }_{\mathrm{1}} and {\varphi }_{\mathrm{2}} are given by Eq. (45), and {\theta }_{\mathrm{1}} and {\theta }_{\mathrm{2}} are angles between the x_{\mathrm{1}}-axis and a straight line passing through the scaling center and Nodes 1 and 2, respectively (Fig. 3(b)). The approximate solution along the boundary direction for this particular element is obtained by interpolating the nodal functions at the two nodes using the linear shape functions given by Eq. (45). It should be evident that the 2-node circular-arc elements can be used to describe the domain geometry exactly when the defining curve consists of only circular arcs. By employing the parameterization Eq. (48) along with choosing the origin of the x_{\mathrm{1}}-x_{\mathrm{2}} coordinate system coincident with the scaling center, the matrices b_{\mathrm{1}}^h and b_{\mathrm{2}}^h and J^h for this case become

(49){\mathit{b}}_{\mathrm{1}}^h={\displaystyle \frac{\mathrm{1}}{R}}\left[\begin{array}{l}\cos \theta {\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ \sin \theta {\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],{\mathit{b}}_{\mathrm{2}}^h={\displaystyle \frac{\mathrm{2}}{R({\theta }_{\mathrm{2}}-{\theta }_{\mathrm{1}})}}\left[\begin{array}{l}-\sin \theta {\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\\ \cos \theta {\mathit{I}}_{\mathrm{\Lambda }\times \mathrm{\Lambda }}\end{array}\right],J^h={\displaystyle \frac{R^{\mathrm{2}}({\theta }_{\mathrm{2}}-{\theta }_{\mathrm{1}})}{\mathrm{2}}}\mathrm{.}

As evident from Eq. (49), J^h is constant whereas b_{\mathrm{1}}^h and b_{\mathrm{2}}^h involve simple elementary functions. This together with the linear approximation along the boundary direction allows the matrices E₀, E₁ and E₂ to be constructed analytically. To further enhance the approximation of the solution along the boundary direction, higher-order Lagrange basis functions can be employed by simply adding interior nodes while the geometry is still described by the exact representation Eq. (48).

Numerical results and discussion

In this section, results of three examples are reported not only to verify the proposed technique but also to demonstrate its computational performance when the discretization of the defining curve is exact. Representative problems of heat conduction (\mathrm{\Lambda }=\mathrm{1}), linear elasticity (\mathrm{\Lambda }=\mathrm{2}) and linear piezoelectricity (\mathrm{\Lambda }=\mathrm{3}) with bodies possessing circular boundaries and subjected to various scenarios are chosen in the numerical study to mainly serve the latter purpose and to additionally exhibit the capability of the technique to treat general boundary conditions and prescribed data on the side-faces. Two types of elements described in the previous section are employed in the discretization and results are then compared to investigate the quality of approximate solutions when the exact description of the defining curve is utilized. For brevity, the symbols “Type-1” and “Type-2” are later used to indicate 2-node linear elements and 2-node circular-arc elements, respectively. Accuracy and convergence of computed numerical solutions are confirmed by benchmarking with available reference solutions and carrying out the analysis via a series of meshes with different levels of refinement, respectively.

Heat conduction problem

First, a representative problem associated with the steady-state heat conduction (\mathrm{\Lambda }=\mathrm{1}) is investigated. Consider a two-dimensional domain which is a quarter of a ring of inner radius R₁ and outer radius R₂ and occupies a region in {ℝ}^{\mathrm{2}} as shown schematically in Fig. 4. The domain is made of a medium with an isotropic heat conductivity k_{\mathrm{0}} (or, equivalently, D_{\mathrm{11}}=D_{\mathrm{22}}=k_{\mathrm{0}},D_{\mathrm{12}}=D_{\mathrm{21}}=\mathrm{0}) and subjected to a non-uniform heat source b_{\mathrm{1}}=Q_{\mathrm{0}}(x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}\mathrm{)/}{R}_{\mathrm{1}}^{\mathrm{2}}, where Q_{\mathrm{0}} is a constant, and mixed boundary conditions as indicated below.

Side AB: u_{\mathrm{1}}^{AB}=-Q_{\mathrm{0}}{x}_{\mathrm{1}}^{\mathrm{4}}\mathrm{/16}{k}_{\mathrm{0}}{R}_{\mathrm{1}}^{\mathrm{2}},

Side BC: u_{\mathrm{1}}^{BC}=-Q_{\mathrm{0}}{R}_{\mathrm{2}}^{\mathrm{4}}\mathrm{/16}{k}_{\mathrm{0}}{R}_{\mathrm{1}}^{\mathrm{2}},

Side CD: t_{\mathrm{1}}^{CD}=\mathrm{0},

Side DA: t_{\mathrm{1}}^{DA}=Q_{\mathrm{0}}{R}_{\mathrm{1}}/\mathrm{4}.

For this particular case, the exact solutions for the temperature field u_{\mathrm{1}}(x_{\mathrm{1}},x_{\mathrm{2}}) and the body heat flux {\sigma }_{\mathrm{11}}(x_{\mathrm{1}},x_{\mathrm{2}}) and {\sigma }_{\mathrm{21}}(x_{\mathrm{1}},x_{\mathrm{2}}) are given by

(50)u_{\mathrm{1}}=-Q_{\mathrm{0}}{(x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}})}^{\mathrm{2}}\mathrm{/16}{k}_{\mathrm{0}}{R}_{\mathrm{1}}^{\mathrm{2}},

(51){\sigma }_{\mathrm{11}}=-Q_{\mathrm{0}}(x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}})x_{\mathrm{1}}/\mathrm{4}{R}_{\mathrm{1}}^{\mathrm{2}},{\sigma }_{\mathrm{21}}=-Q_{\mathrm{0}}(x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}})x_{\mathrm{2}}/\mathrm{4}{R}_{\mathrm{1}}^{\mathrm{2}}\mathrm{.}

The domain geometry is described by the scaling center at Point O and the defining curve AD; as a result, the boundaries AB and CD become the side faces (Fig. 4). In the numerical study, R_{\mathrm{2}}/R_{\mathrm{1}}=\mathrm{2}is chosen and the defining curve is discretized by a uniform mesh containing n identical elements of either Type-1 or Type-2. Results for the normalized temperature along the circular arc between AD and BC (i.e., along \sqrt{x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}}=(R_{\mathrm{1}}+R_{\mathrm{2}})/\mathrm{2}) are normalized by the existing analytical solution and then reported in Table 1 at different values of normalized angle \theta \mathrm{/\pi }={\tan }^{-\mathrm{1}}(x_{\mathrm{2}}/x_{\mathrm{1}})\mathrm{/\pi } and for four different meshes and both types of elements. It is seen that numerical solutions generated by the proposed technique with Type-1 elements converge to the reference solution as the mesh is refined whereas those obtained from the Type-2 elements exhibit excellent agreement with the exact solution for all levels of refinement considered. It is worth noting that numerical solutions generated from the Type-2 elements are much more accurate than those from the Type-1 elements for the same level of refinement. To additionally confirm this argument, the relative errors of scaled boundary finite element solutions obtained for both types of elements are computed for various levels of discretization and results are shown in Fig. 5. As indicated by this set of results, using the exact defining curve in the solution procedure yields significantly less solution error for all levels of refinement and obtained approximate solutions are therefore nearly identical to the exact solution. This is due mainly to that the geometry of the body, prescribed data and all involved differential operators are treated exactly. As the direct result, using circular-arc elements with quadratic element shape functions in the discretization of solution along the boundary direction yields nearly identical relative error to that for the case of Type-2 elements.

Linear elasticity problem

Consider, next, a representative problem in linear elasticity (\mathrm{\Lambda }=\mathrm{2}) corresponding to a plane-strain, elastic, hollowed disk with the inner radius and outer radius denoted by R₁ and R₂, respectively, as shown in Fig. 6(a). The body is made of a homogeneous, isotropic, linearly elastic material with Young’s modulus E and Poisson’s ratio \nu and free of the body force field. On the inner boundary, the normal traction vanishes and the uniform shear traction {\tau }_{\mathrm{0}} is fully prescribed whereas on the outer boundary, the body is fully restrained against the movement. For this particular problem, the exact solution for the displacement and stress fields can be readily obtained from the theory of linear elasticity and they are given explicitly by

(52)u_{\mathrm{1}}=-{\displaystyle \frac{{\tau }_{\mathrm{0}}}{\mathrm{2}\mu }}{\overline{R}}_{\mathrm{1}}^{\mathrm{2}}{R}_{\mathrm{2}}\left({\overline{x}}_{\mathrm{2}}-{\displaystyle \frac{{\overline{x}}_{\mathrm{2}}}{{\overline{x}}_{\mathrm{1}}^{\mathrm{2}}+{\overline{x}}_{\mathrm{2}}^{\mathrm{2}}}}\right),u_{\mathrm{2}}={\displaystyle \frac{{\tau }_{\mathrm{0}}}{\mathrm{2}\mu }}{\overline{R}}_{\mathrm{1}}^{\mathrm{2}}{R}_{\mathrm{2}}\left({\overline{x}}_{\mathrm{1}}-{\displaystyle \frac{{\overline{x}}_{\mathrm{1}}}{{\overline{x}}_{\mathrm{1}}^{\mathrm{2}}+{\overline{x}}_{\mathrm{2}}^{\mathrm{2}}}}\right),

(53){\sigma }_{\mathrm{11}}=-{\sigma }_{\mathrm{22}}=-\mathrm{2}{\tau }_{\mathrm{0}}{\overline{R}}_{\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{{\overline{x}}_{\mathrm{1}}{\overline{x}}_{\mathrm{2}}}{{({\overline{x}}_{\mathrm{1}}^{\mathrm{2}}+{\overline{x}}_{\mathrm{2}}^{\mathrm{2}})}^{\mathrm{2}}}},{\sigma }_{\mathrm{12}}={\sigma }_{\mathrm{21}}={\tau }_{\mathrm{0}}{\overline{R}}_{\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{{\overline{x}}_{\mathrm{1}}^{\mathrm{2}}-{\overline{x}}_{\mathrm{2}}^{\mathrm{2}}}{{({\overline{x}}_{\mathrm{1}}^{\mathrm{2}}+{\overline{x}}_{\mathrm{2}}^{\mathrm{2}})}^{\mathrm{2}}}},

where \mu is the shear modulus, {\bar{R}}_{\mathrm{1}}=R_{\mathrm{1}}/R_{\mathrm{2}}, {\bar{x}}_{\mathrm{1}}=x_{\mathrm{1}}/R_{\mathrm{2}} and {\bar{x}}_{\mathrm{2}}=x_{\mathrm{2}}/R_{\mathrm{2}}. Due to the rotational symmetry and to demonstrate the capability of the present method to treat mixed boundary conditions and prescribed non-zero data on the side-faces, a quarter of the body, as illustrated in Fig. 6(b), is considered. On four sides of the reduced domain, following non-uniform mixed boundary conditions are prescribed:

Side AB:

t_{\mathrm{1}}^{AB}=-{\tau }_{\mathrm{0}}{\displaystyle \frac{{\bar{R}}_{\mathrm{1}}^{\mathrm{2}}}{{\bar{x}}_{\mathrm{1}}^{\mathrm{2}}}},u_{\mathrm{2}}^{AB}={\displaystyle \frac{{\tau }_{\mathrm{0}}}{\mathrm{2}\mu }}{\bar{R}}_{\mathrm{1}}^{\mathrm{2}}{R}_{\mathrm{2}}\left({\bar{x}}_{\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{{\bar{x}}_{\mathrm{1}}}}\right)

,

Side BC: u_{\mathrm{1}}^{BC}=u_{\mathrm{2}}^{BC}=\mathrm{0},

Side CD:

u_{\mathrm{1}}^{CD}=-{\displaystyle \frac{{\tau }_{\mathrm{0}}}{\mathrm{2}\mu }}{\bar{R}}_{\mathrm{1}}^{\mathrm{2}}{R}_{\mathrm{2}}\left({\bar{x}}_{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{{\bar{x}}_{\mathrm{2}}}}\right),t_{\mathrm{2}}^{CD}={\tau }_{\mathrm{0}}{\displaystyle \frac{{\bar{R}}_{\mathrm{1}}^{\mathrm{2}}}{{\bar{x}}_{\mathrm{2}}^{\mathrm{2}}}}

,

Side AD: t_{\mathrm{1}}^{AD}={\tau }_{\mathrm{0}}{\displaystyle \frac{x_{\mathrm{2}}}{R_{\mathrm{1}}}},t_{\mathrm{2}}^{AD}=-{\tau }_{\mathrm{0}}{\displaystyle \frac{x_{\mathrm{1}}}{R_{\mathrm{1}}}}.

To describe the geometry, the scaling center is chosen at the center of the hollowed disk whereas its inner boundary is treated as the defining curve. In numerical simulations, R_{\mathrm{2}}/R_{\mathrm{1}}=\mathrm{2}is chosen and uniform meshes containing n elements of either Type-1 or Type-2 are employed to discretize the defining curve of the reduced domain. Numerical solutions for the displacement components along the circular arc between AD and BC (i.e., along \sqrt{x_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}^{\mathrm{2}}}=(R_{\mathrm{1}}+R_{\mathrm{2}})/\mathrm{2}) are shown in Tables 2 and 3 at different values of normalized angle \theta \mathrm{/\pi }={\tan }^{-\mathrm{1}}(x_{\mathrm{2}}/x_{\mathrm{1}})\mathrm{/\pi } and for different levels of refinements and two types of elements used in the discretization of the defining curve. Besides the good agreement between the scaled boundary finite element and benchmark solutions and the good convergence behaviour upon the mesh refinement, using the 2-node circular-arc elements in the discretization of the defining curve yields more accurate solutions. The improvement of the solution accuracy can also be confirmed from the relationship between the relative error and the degrees of freedom used in the discretization shown in Fig. 7. Similar to the previous case, the use of the exact defining curve does not alter the rate of convergence as the mesh is refined but significantly reducing the solution error. Note also that if the quadratic circular-arc elements are employed in the solution discretization, the rate of convergence is apparently enhanced and, for the same number of degrees of freedom, the relative error is less than that for the case of Type-2 elements.

Linear piezoelectricity problem

As a final example, a representative boundary value problem associated with the linear piezoelectricity (

\Lambda =\mathrm{3}

) is investigated to further emphasize the capability of the proposed technique to treat linear multi-field problems. Let us consider a quarter of hollowed circular plate with the inner radius and outer radius denoted by R₁ and R₂, respectively, and made of a linear piezoelectric solid with all material constants taken from PZT-4 [⁸] as shown in Fig. 8. In particular, all non-zero entries of the modulus matrix D are given explicitly by D_{\mathrm{11}}=\mathrm{139}\mathrm{GPa}, D_{\mathrm{15}}=D_{\mathrm{51}}=\mathrm{74.3}\mathrm{GPa}, D_{\mathrm{16}}=D_{\mathrm{61}}=-\mathrm{6.98}{\mathrm{C/m}}^2, D_{\mathrm{22}}=D_{\mathrm{24}}=\mathrm{25.6}\mathrm{GPa}, D_{\mathrm{42}}=D_{\mathrm{44}}=\mathrm{25.6}\mathrm{GPa}, D_{\mathrm{55}}=\mathrm{113}\mathrm{GPa}, D_{\mathrm{33}}=\mathrm{6.0}\times {\mathrm{10}}^{-\mathrm{9}}\mathrm{C/(V}\cdot \mathrm{m)}, D_{\mathrm{34}}=D_{\mathrm{43}}=\mathrm{13.44}{\mathrm{C/m}}^2, D_{\mathrm{56}}=D_{\mathrm{65}}=\mathrm{13.84}{\mathrm{C/m}}^2, and D_{\mathrm{66}}=-\mathrm{5.47}\times {\mathrm{10}}^{-\mathrm{9}}\mathrm{C/(V}\cdot \mathrm{m)}. For the purpose of verification, a set of prescribed data such as the distributed body source and boundary conditions on four sides is chosen such that the exact solution for the displacements (u_{\mathrm{1}},u_{\mathrm{2}}) and the electric potential (u_{\mathrm{3}}) takes the following form:

(54)u_{\mathrm{1}}=u_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}^{\mathrm{2}},u_{\mathrm{2}}=u_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}^{\mathrm{2}},u_{\mathrm{3}}={\varphi }_{\mathrm{0}}{\overline{x}}_{\mathrm{2}},

where u_{\mathrm{10}}, u_{\mathrm{20}}, {\phi }_{\mathrm{0}} are given constants, {\bar{x}}_{\mathrm{1}}=x_{\mathrm{1}}/R and {\bar{x}}_{\mathrm{2}}=x_{\mathrm{2}}/R with R taken as R_{\mathrm{1}}. The corresponding exact stress field ({\sigma }_{\mathrm{11}},{\sigma }_{\mathrm{12}},{\sigma }_{\mathrm{21}},{\sigma }_{\mathrm{22}}) and electrical induction vector ({\sigma }_{\mathrm{13}},{\sigma }_{\mathrm{23}}) are given by

(55){\sigma }_{\mathrm{11}}=(\mathrm{2}{D}_{\mathrm{11}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{15}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{16}}{\varphi }_{\mathrm{0}})/R,{\sigma }_{\mathrm{22}}=(\mathrm{2}{D}_{\mathrm{15}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{55}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{56}}{\varphi }_{\mathrm{0}})/R,{\sigma }_{\mathrm{12}}={\sigma }_{\mathrm{21}}=\mathrm{0},

(56){\sigma }_{\mathrm{13}}=\mathrm{0},{\sigma }_{\mathrm{23}}=(\mathrm{2}{D}_{\mathrm{16}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{56}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{66}}{\varphi }_{\mathrm{0}})/R\mathrm{.}

The distributed body source required to ensure the equilibrium state is given by

(57)b_{\mathrm{1}}=-\mathrm{2}{D}_{\mathrm{11}}{u}_{\mathrm{10}}/R^{\mathrm{2}},b_{\mathrm{2}}=-\mathrm{2}{D}_{\mathrm{55}}{u}_{\mathrm{20}}/R^{\mathrm{2}},b_{\mathrm{3}}=-\mathrm{2}{D}_{\mathrm{56}}{u}_{\mathrm{20}}/R^{\mathrm{2}}\mathrm{.}

In the analysis, the scaling center is chosen at a Point O and the prescribed conditions on the side faces AB and DC and the boundaries AD and BC are given below:

Side AB: t_{\mathrm{1}}^{AB}=\mathrm{0},u_{\mathrm{2}}^{AB}=\mathrm{0},u_{\mathrm{3}}^{AB}=\mathrm{0},

Side CD: u_{\mathrm{1}}^{CD}=\mathrm{0},u_{\mathrm{2}}^{CD}=u_{\mathrm{20}}{\bar{x}}_{\mathrm{2}}^{\mathrm{2}},t_{\mathrm{3}}^{CD}=\mathrm{0},

Side DA:

t_{\mathrm{1}}^{DA}=-(\mathrm{2}{D}_{\mathrm{11}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{15}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{16}}{\phi }_{\mathrm{0}}){\overline{x}}_{\mathrm{1}}/R,

t_{\mathrm{2}}^{DA}=-(\mathrm{2}{D}_{\mathrm{15}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{55}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{56}}{\phi }_{\mathrm{0}}){\overline{x}}_{\mathrm{2}}/R,t_{\mathrm{3}}^{DA}=-(\mathrm{2}{D}_{\mathrm{16}}{u}_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}+\mathrm{2}{D}_{\mathrm{56}}{u}_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}+D_{\mathrm{66}}{\phi }_{\mathrm{0}}){\overline{x}}_{\mathrm{2}}/R,

Side BC: u_{\mathrm{1}}^{BC}=u_{\mathrm{10}}{\overline{x}}_{\mathrm{1}}^{\mathrm{2}},u_{\mathrm{2}}^{BC}=u_{\mathrm{20}}{\overline{x}}_{\mathrm{2}}^{\mathrm{2}},u_{\mathrm{3}}^{BC}={\phi }_{\mathrm{0}}{\overline{x}}_{\mathrm{2}}.

The prescribed boundary conditions indicated above are chosen mainly to demonstrate the capability of the present technique to deal with general prescribed data on the side faces and the boundary of the domain. In numerical simulations, uniform meshes containing n elements of either Type-1 or Type-2 are employed to discretize the defining curve (i.e., the boundary DA) and R_{\mathrm{2}}/R_{\mathrm{1}}=\mathrm{2}, u_{\mathrm{20}}/u_{\mathrm{10}}=\mathrm{2} and D_{\mathrm{16}}{\phi }_{\mathrm{0}}/D_{\mathrm{11}}{u}_{\mathrm{10}}=\mathrm{1} are chosen. Computed displacements, electric potentials along the line x_{\mathrm{1}}=x_{\mathrm{2}}, normalized by the analytical solution, are reported at different values of {\overline{x}}_{\mathrm{1}}=x_1/Rin Tables 4–6 for various levels of discretization and both types of elements. It is seen, for this particular problem, that the proposed technique yields highly accurate results even when relatively coarse meshes containing only few degrees of freedom are employed. Again, use of the exact defining curve to describe the domain geometry and its boundary can significantly improve the accuracy of the numerical solutions for the entire domain and allows coarser meshes to be utilized to obtain the same level of accuracy. This is confirmed by the plot of the relative error versus the number of degrees of freedom used in the discretization of the defining curve shown in Fig. 8. Similar to the previous problem, the rate of convergence of the numerical solution can be improved and the relative error can be significantly reduced by using the quadratic circular-arc elements.

Conclusions and remarks

An efficient numerical technique based on the scaled boundary finite element method has been successfully implemented for the analysis of two-dimensional, linear, boundary value problems. Both the formulation and numerical implementations have been established in a general framework allowing a variety of linear field problems (e.g., potential and Laplace-type problems, linear elasticity, linear coupled-field problems, etc.) and the general associated data (e.g., domain geometry, prescribed distributed body source, boundary conditions, contribution of prescribed conditions on the side-faces, etc.) to be treated in a single, unified fashion. The exact treatment of the defining curve used for describing the geometry of the domain and the corresponding boundaries and its influence on the improvement of solution accuracy has also been fully investigated via a special class of boundary value problems concerning domains with a circular boundary. A two-node circular arc element has been developed and implemented to serve that purpose and obtained results have been benchmarked with those generated by standard 2-node linear elements. Numerical results obtained from extensive numerical experiments have indicated that the proposed SBFEM yields highly accurate solutions for various types of boundary value problems considered and, in particular, only few degrees of freedom are required along the defining curve to attain such high-quality results. In addition, it has been found that use of the exact description of the defining curve in the solution procedure can significantly reduce the solution error and, as a result, reduce the number of degrees of freedom required to achieve the target accuracy in comparison with standard linear elements.

References

[1] Wolf J P. The Scaled Boundary Finite Element Method. Chichester: John Wiley & Sons, 2003

[2] Wolf J P, Song C. Finite-Element Modelling of Unbounded Domain. Chichester: John Wiley & Sons, 1996

[3] Deeks J A, Wolf J P. A virtual work derivation of the scaled boundary finite-element method for elastostatics. Computational Mechanics, 2002, 28(6): 489–504

[4] Cruse T A. Boundary Element Analysis in Computational Fracture Mechanics. Dordrecht: Kluwer Academic Publishers, 1988

[5] Brebbia C A, Dominguez J. Boundary Elements: An Introductory Course. 2nd ed. New York: McGraw-Hill, 1989

[6] Bonnet M, Maier G, Polizzotto C. Symmetric Galerkin Boundary Element Methods. Applied Mechanics Reviews, 1998, 51(11): 669–703

[7] Liu J, Lin G A. A scaled boundary finite element method applied to electrostatic problems. Engineering Analysis with Boundary Elements, 2012, 36(12): 1721–1732

[8] Li C, Man H, Song C, Gao W. Fracture analysis of piezoelectric materials using the scaled boundary finite element method. Engineering Fracture Mechanics, 2013, 97: 52–71

[9] Vu T H, Deeks A J. Using fundamental solutions in the scaled boundary finite element method to solve problems with concentrated loads. Computational Mechanics, 2014, 53(4): 641–657

[10] Ooi E T, Song C, Tin-Loi F. A scaled boundary polygon formulation for elasto-plastic analyses. Computer Methods in Applied Mechanics and Engineering, 2005, 268: 905–937

[11] Doherty J P, Deeks A J. Adaptive coupling of the finite-element and scaled boundary finite-element methods for non-linear analysis of unbounded media. Computers and Geotechnics, 2015, 32(6): 436–444

[12] Li F, Ren P. A novel solution for heat conduction problems by extending scaled boundary finite element method. International Journal of Heat and Mass Transfer, 2016, 95: 678–688

[13] Li M, Zhang H, Guan H. Study of offshore monopole behavior due to ocean waves. Ocean Engineering, 2011, 38(17–18): 1946–1956

[14] Meng X N, Zou Z J. Radiation and diffraction of water waves by an infinite horizontal structure with a sidewall using SBFEM. Ocean Engineering, 2013, 60: 193–199

[15] Gravenkamp H, Birk C, Song C. The computation of dispersion relations for axisymmetric waveguides using the scaled boundary finite element method. Ultrasonics, 2014, 54(5): 1373–1385

[16] Li C, Ooi E T, Song C, Natarajan S. SBFEM for fracture analysis of piezoelectric composites under thermal load. International Journal of Solids and Structures, 2015, 52: 114–129

[17] Song C, Wolf J P. The scaled boundary finite-element method—alias consistent infinitesimal finite-element cell method—for elastodynamics. Computer Methods in Applied Mechanics and Engineering, 1997, 147(3–4): 329–355

[18] Wolf J P, Song C. The scaled boundary finite-element method: A fundamental solution-less boundary-element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190(42): 5551–5568

[19] Deeks A J. Prescribed side-face displacements in the scaled boundary finite-element method. Computers & Structures, 2004, 82(15–16): 1153–1165

[20] Song C, Wolf J P. Body loads in scaled boundary finite-element method. Computer Methods in Applied Mechanics and Engineering, 1999, 180(1–2): 117–135

[21] He Y, Yang H, Deeks A J. An element-free Galerkin (EFG) scaled boundary method. Finite Elements in Analysis and Design, 2012, 62: 28–36

[22] Vu T H, Deeks A J. Use of higher-order shape functions in the scaled boundary finite element method. International Journal for Numerical Methods in Engineering, 2006, 65(10): 1714–1733

[23] He Y, Yang H, Deeks A J. Use of Fourier shape functions in the scaled boundary method. Engineering Analysis with Boundary Elements, 2014, 41: 152–159

[24] Deeks A J, Wolf J P. An h-hierarchical adaptive procedure for the scaled boundary finite-element method. International Journal for Numerical Methods in Engineering, 2002, 54(4): 585–605

[25] Vu T H, Deeks A J. A p-adaptive scaled boundary finite element method based on maximization of the error decrease rate. Computational Mechanics, 2007, 41(3): 441–455

[26] Deeks A J, Augarde C E. A meshless local Petrov-Galerkin scaled boundary method. Computational Mechanics, 2005, 36(3): 159–170

[27] Chung N V. Analysis of two-dimensional linear field problems by scaled boundary finite element method. Dissertation for the Doctoral Degree. Bangkok: Chulalongkorn University, 2016

[28] Chung N V, Jaroon R, Phoonsak P. Scaled boundary finite element method for two-dimensional linear multi-field media. Engineering Journal (Thailand), 2017, 21(7): 334–360

[29] Ooi E T, Song C, Tin-Loi F, Yang Z J. Automatic modelling of cohesive crack propagation in concrete using polygon scaled boundary finite elements. Engineering Fracture Mechanics, 2012, 93: 13–33

[30] Ooi E T, Shi C, Song C, Tin-Loi F, Yang Z J. Dynamic crack propagation simulation with scaled boundary polygon elements and automatic remeshing technique. Engineering Fracture Mechanics, 2013, 106: 1–21

[31] Dieringer R, Becker W. A new scaled boundary finite element formulation for the computation of singularity orders at cracks and notches in arbitrarily laminated composites. Composite Structures, 2015, 123: 263–270

[32] Natarajan S, Wang J, Song C, Birk C. Isogeometric analysis enhanced by the scaled boundary finite element method. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 733–762

[33] Nguyen B H, Tran H D, Anitescu C, Zhuang X, Rabczuk T. Isogeometric symmetric Galerkin boundary element method for three-dimensional elasticity problems. Computer Methods in Applied Mechanics and Engineering, 2017, 323: 132–150

[34] Nguyen B H, Tran H D, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 252–275

[35] Li P, Liu J, Lin G, Zhang P, Xu B. A combination of isogeometric technique and scaled boundary method for solution of the steady-state heat transfer problems in arbitrary plane domain with Robin boundary. Engineering Analysis with Boundary Elements, 2017, 82: 43–56

[36] Li F, Qiang T. The scaled boundary finite element analysis of seepage problems in multi-material regions. International Journal of Computational Methods, 2012, 9(1): 1240008

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature