Numerical investigation of circle defining curve for two-dimensional problem with general boundaries using the scaled boundary finite element method

Abstract

The scaled boundary finite element method (SBFEM) is applied to the static analysis of two dimensional elasticity problem, boundary value problems domain with the domain completely described by a circular defining curve. The scaled boundary finite element equations is formulated within a general framework integrating the influence of the distributed body force, general boundary conditions, and bounded and unbounded domain. This paper investigates the possibility of using exact geometry to form the exact description of the circular defining curve and the standard finite element shape function to approximate the defining curve. Three linear elasticity problems are presented to verify the proposed method with the analytical solution. Numerical examples show the accuracy and efficiency of the proposed method, and the performance is found to be better than using standard linear element for the approximation defining curve on the scaled boundary method.

Keywords

exact geometry / circular defining curve / general boundaries / SBFEM

Introduction

The scaled boundary finite element method is recognized as a semi-analytical technique combining features of both analytical schemes and the finite element approximation. The SBFEM is achieved in two purposes such with regards to the analytical and numerical method and to the standard procedure of the finite element and boundary element method within the numerical procedures [¹]. In last two decades, the SBFEM has been developed for unbounded and bounded domains in two and three-dimensional media. The method was originally derived to compute the dynamic stiffness of the unbounded domain [²]. The SBFEM has proved to be more general than initially investigated, then developments has allowed analysis of incompressible material and bounded domain [³], and inclusion of body loads [⁴]. The complexity of the original derivation of this technique led to develop weighted residual formulation [⁵,⁶]. Then Deeks and Wolf [⁷] and Deeks [⁸] commenced with virtual work and novel semi-analytical approach of the scaled boundary finite element method to derive the standard finite element method for two dimensional problems in solid mechanics accessibly.

Deeks and Wolf [⁹] used an h-hierachical adaptive procedure in the SBFEM. This technique took the ability of the SBFEM to model stress singularities at the scaling centre and to avoid discretization of certain adjacent segments of the boundary. Vu and Deeks [¹⁰] investigated high-order elements in the SBFEM. The spectral element and hierarchical approach were employed in this study. They found that the spectral element approach was better than the hierarchical approach. Doherty and Deeks [¹¹] developed a meshless scaled boundary method to model the far field and the conventional meshless local Petrov-Galerkin modelling. This combining was general and could be employed to other techniques of modelling the far field. Although, the SBFEM has demonstrated many advantages in the approach method, it also has had disadvantaged in solving problems involving an unbounded domain or stress singularities. When the number of degrees of freedom became too large, the computational expense was a trouble. So, Vu and Deeks [¹²] developed a p-adaptive in the SBFEM for the two dimensional problem. These authors investigated an alternative set of refinement criteria. This led to be maximize the solution accuracy and minimizing the cost. Furthermore, He et al. [¹³] presented a new Element-free Galerkin scaled boundary method to approximate in the circumferential direction. This work was applied to a number of standard linear elasticity problems, and the technique was found to offer higher and better convergence than the original SBFEM. Additionally, He et al. [¹⁴] investigated the possibility of using the Fourier shape functions in the SBFEM to approximate in the circumferential direction. This research used to solve three elastostactic and steady-state heat transfer problems. They found that the accuracy and convergence were better than using polynomial elements or using an element-free Galerkin to approximate on the circumferential direction in the SBFEM. These published papers focused on the shape functions and applied in individual particular problem to show the advantages of the SBFEM. They may not provide the advantages of the using other kinds of shape function for general problem. They only used typical function to apply in the SBFEM for each particular problem such as using linear shape functions for heat problem, two dimensional problems, etc. In nearly years, Vu and Deeks [¹⁵] used fundamental solutions in the scaled boundary finite element method to solve problems with concentrated loads. Liu and Lin [¹⁶] extended the SBFEM to treat electrostatic problems. He et al. [¹⁷] presented an approach method to develop for numerical analysis of 2D elastic systems with rotationally periodic symmetry under arbitrary load conditions. Ooi et al. [¹⁸,¹⁹] developed an efficient methodology for automatic dynamic crack propagation simulations using scaled boundary polygon elements. In the other hand, the accuracy of the SBFEM relies on the meshing of geometric and approximate state variables. If mesh refinement is required, the geometry representations in Computer Aided Design is difficult. Recently, isogeometric analysis (IGA) has become popular in the finite element method, boundary element method. The most feature of this variant is that uses non-uniform rational B-splines as the basic functions to describe the geometry and approximate state variables [²⁰–²²]. IGA has been successfully applied to e.g. three-dimensional elasticity problems and two-dimensional crack problem [²³,²⁴].

The aforementioned works have shown various important progresses to implement the SBFEM in analysis of engineering problems. However, the existing methods are nearly all focused on the analysis of the structural response through the SBFEM. The alternative computational procedure was only solved for individual particular boundary conditions and presented each advantage of the SBFEM, while less work has been conducted for its advantages to apply for general boundary conditions. In the other hand, they used shape functions in the circumferential direction to investigate the influence on the accuracy of the resulting solutions. In this present study, the SBFEM is developed for solving two-dimensional elasticity problem, general boundary conditions. Basic field equations governing are formulated in a general setting allowing various types of shape to be treated in the same manner. In addition, the conventional polar coordinates are investigated to discretize on the defining curve.

Problem and scaled boundary formulation

Consider a two-dimensional elasticity body occupying a domain \mathbf{\Omega } in {\mathrm{ℝ}}^{\mathrm{2}}as shown schematically in Fig. 1. The domain is assumed smooth in the sense that all involved mathematical operators (e.g., integrations and differentiations) can be performed over this domain. In addition, the boundary of the body\mathbf{\Omega }, denoted by\partial \mathbf{\Omega }, is assumed piecewise smooth and an outward unit normal vector at any smooth point on \partial \mathbf{\Omega } is denoted by\mathit{n}={\{n_{\mathrm{1}}{n}_{\mathrm{2}}\}}^{\mathrm{T}}. The interior of the body is denoted by\mathrm{int}\mathbf{\Omega }.

Three basic 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, which relate the three field quantitiesu(x),\overline{\epsilon }(\mathit{x})and \sigma (\mathit{x}), are given explicitly by

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

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

(3)\overline{\epsilon }=\mathit{L}\mathit{u},

where L is a linear differential operator defined, in terms of a \mathrm{4}\times \mathrm{2}-matrix, by

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

with I and 0 denoting a \mathrm{2}\times \mathrm{2}-identity matrix and a \mathrm{2}\times \mathrm{2}-zero matrix, respectively. By applying the law of conservation at any smooth point \mathit{x} on the boundary\partial \mathbf{\Omega }, the surface flux \mathit{t}(\mathit{x})can be related to the body flux\sigma (\mathit{x})and the outward unit normal vector\mathit{n}(\mathit{x})by\mathit{t}=\left[\begin{array}{ll}{n}_{\mathrm{1}}\mathit{I}& n_{\mathrm{2}}\mathit{I}\end{array}\right]\sigma ,where n_{\mathrm{1}} and n_{\mathrm{2}}are components of\mathit{n}(\mathit{x}).

By applying the standard weighted residual technique to the law of conservation (1), then integrating certain integral by parts via Gauss-divergence theorem, and finally employing the relations (2) and (3), the weak-form equation in terms of the state variable is given by

(5){\displaystyle \underset{\mathrm{\Omega }}{\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{\Omega }}{\int }{\mathit{w}}^{\mathrm{T}}\mathit{b}}\mathrm{d}A,

where\mathit{w}is a weight function.

Scaled boundary finite element equations

Let{\mathit{x}}_{\mathrm{0}}=(x_{\mathrm{10}},x_{\mathrm{20}})be a point in {\mathrm{ℝ}}^{\mathrm{2}} and C be a simple, piecewise smooth curve in{\mathrm{ℝ}}^{\mathrm{2}} parameterized by a function\left(x_{\mathrm{10}}+{\widehat{x}}_{\mathrm{1}}(s),x_{\mathrm{20}}+{\widehat{x}}_{\mathrm{2}}(s)\right)\in {\mathrm{ℝ}}^{\mathrm{2}}\mathit{r}:s\in [a,b]\toas shown in Fig. 2. Now, let us introduce the following coordinate transformation

(6)x_{\alpha }=x_{\alpha \mathrm{0}}+\xi {\widehat{x}}_{\alpha }(s).

In this present study, two types of elements used to describe the domain geometry and to discretize a 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. For the 2-node isoparametric linear element, the approximate defining curve is parameterized by

(7){\widehat{x}}_{\mathrm{1}}(s)=x_{\mathrm{11}}{\displaystyle \frac{(\mathrm{1}-s)}{\mathrm{2}}}+x_{\mathrm{12}}{\displaystyle \frac{(\mathrm{1}+s)}{\mathrm{2}}},

(8){\widehat{x}}_{\mathrm{2}}(s)=x_{\mathrm{21}}{\displaystyle \frac{(\mathrm{1}-s)}{\mathrm{2}}}+x_{\mathrm{22}}{\displaystyle \frac{(\mathrm{1}+s)}{\mathrm{2}}},

wheres\in [-\mathrm{1},\mathrm{1}] denotes the local or element boundary coordinate. For the 2-node circular-arc element, the defining curve is parameterized by

(9){\widehat{x}}_{\mathrm{1}}(s)=\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\left({\mathit{\theta }}_a{\displaystyle \frac{(\mathrm{1}-s)}{\mathrm{2}}}+{\mathit{\theta }}_b{\displaystyle \frac{(\mathrm{1}+s)}{\mathrm{2}}}\right),

(10){\widehat{x}}_{\mathrm{2}}(s)=\mathrm{r}\sin \left({\mathit{\theta }}_a{\displaystyle \frac{(\mathrm{1}-s)}{\mathrm{2}}}+{\mathit{\theta }}_b{\displaystyle \frac{(\mathrm{1}+s)}{\mathrm{2}}}\right).

The linear differential operator \mathit{L} given by Eq. (4) can now be expressed in terms of partial derivatives with respect to\xi andsby

(11)\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{4}\times \mathrm{2}-matrices defined by

(12){\mathit{b}}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{J}}\left[\begin{array}{l}-{\widehat{x}}_{\mathrm{2}}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\\ {\widehat{x}}_{\mathrm{1}}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right];J={\widehat{x}}_{\mathrm{1}}{\displaystyle \frac{\mathrm{d}{\widehat{x}}_{\mathrm{2}}}{\mathrm{d}s}}-{\widehat{x}}_{\mathrm{2}}{\displaystyle \frac{\mathrm{d}{\widehat{x}}_{\mathrm{1}}}{\mathrm{d}s}}.{\mathit{b}}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{J}}\left[\begin{array}{l}{\displaystyle \frac{\mathrm{d}{\widehat{x}}_{\mathrm{2}}}{\mathrm{d}s}}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\\ -{\displaystyle \frac{\mathrm{d}{\widehat{x}}_{\mathrm{1}}}{\mathrm{d}s}}{\mathit{I}}_{\mathrm{2}\times \mathrm{2}}\end{array}\right];

From the coordinate transformation along with the approximation, the state variable\mathit{u}is now approximated by{\mathit{u}}^hin a form

(13){\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,

where{\mathit{u}}_{(i)}^h(\xi )denotes the value of the state variable along the lines=s_{(i)},{\mathit{N}}^Sis a \mathrm{2}\times m\mathrm{2}-matrix containing all basis functions, and {\mathit{U}}^his a vector containing all functions{\mathit{u}}_{(i)}^h(\xi ).The approximation of the body flux\mathit{\sigma }is given by

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

where B₁ and B₂ are given by

(15){\mathit{B}}_{\mathrm{1}}={\mathit{b}}_{\mathrm{1}}^{}{\mathit{N}}^S;{\mathit{B}}_{\mathrm{2}}={\mathit{b}}_{\mathrm{2}}^{}\mathrm{d}{\mathit{N}}^S/\mathrm{d}s.

Similarly, the weight function \mathit{w} and its derivatives \mathit{L}\mathit{w} can be approximated, in a similar fashion, by

(16){\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,

where{\mathit{w}}_{(i)}^h(\xi )denotes an arbitrary function of the coordinate \xi along the line s=s_{(i)}and {\mathit{W}}^h is a vector containing all functions{\mathit{w}}_{(i)}^h(\xi ).

A set of scaled boundary finite element equations is established for a generic, two-dimensional body\mathbf{\Omega }as shown in Fig. 3. The boundary of the domain \partial \mathbf{\Omega } is assumed consisting of four parts resulting from the scale boundary coordinate transformation with the scaling center x_{\mathrm{0}} and defining curveC: the inner boundary\partial {\mathbf{\Omega }}_{\mathrm{1}}, the outer boundary\partial {\mathbf{\Omega }}_{\mathrm{2}}, the side-face-1\partial {\mathbf{\Omega }}_{\mathrm{1}}^s, and the side-face-2\partial {\mathbf{\Omega }}_{\mathrm{2}}^s. The body is considered in this general setting to ensure that the resulting formulation is applicable to various cases. As a result of the boundary partition\partial \mathbf{\Omega }=\partial {\mathrm{\Omega }}_{\mathrm{1}}\cup \partial {\mathrm{\Omega }}_{\mathrm{2}}\cup \partial {\mathrm{\Omega }}_{\mathrm{1}}^s\cup \partial {\mathrm{\Omega }}_{\mathrm{1}}^s, by changing to the \xi ,s-coordinates via the transformation, the weak-form Eq. (5) 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 }}{(\mathit{L}\mathit{w})}^T\mathit{D}(\mathit{L}\mathit{u})}}J\xi \mathrm{d}\xi \mathrm{d}s=\\ {\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{w}}_{\mathrm{1}}^T{\mathit{t}}_{\mathrm{1}}}(s)J^s(s){\xi }_{\mathrm{1}}\mathrm{d}s+{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\mathit{w}}_{\mathrm{2}}^T{\mathit{t}}_{\mathrm{2}}(s)}{J}^s(s){\xi }_{\mathrm{2}}\mathrm{d}s\\ +{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{({\mathit{w}}_{\mathrm{1}}^s)}^T{\mathit{t}}_{\mathrm{1}}^s(\xi )}{J}_{\mathrm{1}}^{\xi }\mathrm{d}\xi +{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{({\mathit{w}}_{\mathrm{2}}^s)}^T{\mathit{t}}_{\mathrm{2}}^s(\xi )}{J}_{\mathrm{2}}^{\xi }\mathrm{d}\xi \\ +{\displaystyle \underset{s_{\mathrm{1}}}{\overset{s_{\mathrm{2}}}{\int }}{\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{\mathit{w}}^T\mathit{b}}}J\xi \mathrm{d}\xi \mathrm{d}s.\end{array}

By manipulating the involved matrix algebra, integrating the first two integrals by parts with respect to the coordinate

\xi

, the weak-form Eq. (17) is approximated by

(18)$$\begin{gathered} \begin{array}{c}{\displaystyle {\displaystyle \underset{{\xi }_{\mathrm{1}}}{\overset{{\xi }_{\mathrm{2}}}{\int }}{({\mathit{W}}^h)}^{\mathrm{T}}}}[-\xi {\mathit{E}}_0{\mathit{U}}_{,\xi \xi }^h+({\mathit{E}}_1-{\mathit{E}}_1^{\mathrm{T}}-{\mathit{E}}_0){\mathit{U}}_{,\xi }^h+\frac{1}{\xi }{\mathit{E}}_2{\mathit{U}}^h-{\mathit{F}}^t-\xi {\mathit{F}}^b]\mathrm{d}\xi \end{array} \\ +{({\mathit{W}}_2^h)}^{\mathrm{T}}[{\{\xi {\mathit{E}}_0{\mathit{U}}_{,\xi }^h+{\mathit{E}}_1^{\mathrm{T}}{\mathit{U}}^h\}}_{\xi ={\xi }_2}-{\mathit{P}}_2]- \\ {({\mathit{W}}_1^h)}^{\mathrm{T}}[{\{\xi {\mathit{E}}_0{\mathit{U}}_{,\xi }^h+{\mathit{E}}_1^{\mathrm{T}}{\mathit{U}}^h\}}_{\xi ={\xi }_1}+{\mathit{P}}_1]=0, \end{gathered}$$

where the matrices {\mathit{E}}_{\mathrm{0}}, {\mathit{E}}_{\mathrm{1}}, and {\mathit{E}}_{\mathrm{2}} , and the following quantities are defined by

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

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

(21)$$\begin{gathered} {\mathit{F}}_1^t={\mathit{F}}_1^t(\xi )={({\mathit{N}}_1^S)}^{\mathrm{T}}{\mathit{t}}_1^s(\xi )J_1^{\xi } \\ {\mathit{F}}_2^t={\mathit{F}}_2^t(\xi )={({\mathit{N}}_2^S)}^{\mathrm{T}}{\mathit{t}}_2^s(\xi )J_2^{\xi } \\ {\mathit{F}}^t={\mathit{F}}_1^t+{\mathit{F}}_2^t, \end{gathered}$$

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

From the arbitrariness of the weight function{\mathit{W}}^h, it can be deduced that

(23){\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}{\mathit{U}}_{,\xi \xi }^h+\xi ({\mathit{E}}_{\mathrm{0}}+{\mathit{E}}_{\mathrm{1}}^{\mathrm{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=0\forall \xi \in ({\xi }_{\mathrm{1}},{\xi }_{\mathrm{2}}),

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

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

where the vector {\mathit{Q}}^h={\mathit{Q}}^h(\xi ) commonly 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}}^{\mathrm{T}}{\mathit{U}}^h.

Eqs. (23)–(25) form a set of the so-called scaled boundary finite element equations governing the function{\mathit{U}}^h={\mathit{U}}^h(\xi ). It can be remarked that Eq. (23) forms a system of linear, second-order, nonhomogeneous, ordinary differential equations with respect to the coordinate\xiwhereas Eqs. (24)–(25) pose the boundary conditions on the inner and outer boundaries of the body. It should be evident from Eqs. (23)–(25) that the information associated with the prescribed distributed body source and the prescribed boundary conditions on both inner and outer boundaries can be integrated into the formulation via the term {\mathit{F}}^b and the conditions Eqs. (24)‒(25), respectively. Consistent with the partition of the vector {\mathit{U}}^h, the vector {\mathit{F}}^t can also be partitioned into{\mathit{F}}^t={\left\{\begin{array}{ll}{\mathit{F}}^{tu}& {\mathit{F}}^{tc}\end{array}\right\}}^{\mathrm{T}}where {\mathit{F}}^{tu}={\mathit{F}}^{tu}(\xi ) contains many0functions and known functions associated with prescribed surface flux on the side face and has the same dimension as that of {\mathit{U}}^{hu} and {\mathit{F}}^{tc}={\mathit{F}}^{tc}(\xi )contains unknown functions associated with the unknown surface flux on the side face and has the same dimension as that of {\mathit{U}}^{hc}. According to the this partition, the system of differential Eq. (23) can be expressed as

(27)\begin{array}{l}{\xi }^{\mathrm{2}}\left[\begin{array}{ll}{\mathit{E}}_{\mathrm{0}}^{uu}& {\mathit{E}}_{\mathrm{0}}^{uc}\\ {({\mathit{E}}_{\mathrm{0}}^{uc})}^{\mathrm{T}}& {\mathit{E}}_{\mathrm{0}}^{cc}\end{array}\right]\left\{\begin{array}{l}{\mathit{U}}_{,\xi \xi }^{hu}\\ {\mathit{U}}_{,\xi \xi }^{hc}\end{array}\right\}+\\ \xi \left[\begin{array}{ll}{\mathit{E}}_{\mathrm{0}}^{uu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{uu}& {\mathit{E}}_{\mathrm{0}}^{uc}+{({\mathit{E}}_{\mathrm{1}}^{cu})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{uc}\\ {({\mathit{E}}_{\mathrm{0}}^{uc})}^T+{({\mathit{E}}_{\mathrm{1}}^{uc})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{cu}& {\mathit{E}}_{\mathrm{0}}^{cc}+{({\mathit{E}}_{\mathrm{1}}^{cc})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{cc}\end{array}\right]\left\{\begin{array}{l}{\mathit{U}}_{,\xi }^{hu}\\ {\mathit{U}}_{,\xi }^{hc}\end{array}\right\}\\ -\left[\begin{array}{ll}{\mathit{E}}_{\mathrm{2}}^{uu}& {\mathit{E}}_{\mathrm{2}}^{uc}\\ {({\mathit{E}}_{\mathrm{2}}^{uc})}^{\mathrm{T}}& {\mathit{E}}_{\mathrm{2}}^{cc}\end{array}\right]\left\{\begin{array}{l}{\mathit{U}}^{hu}\\ {\mathit{U}}^{hc}\end{array}\right\}+\xi \left\{\begin{array}{l}{\mathit{F}}^{tu}\\ {\mathit{F}}^{tc}\end{array}\right\}+{\xi }^{\mathrm{2}}\left\{\begin{array}{l}{\mathit{F}}^{bu}\\ {\mathit{F}}^{bc}\end{array}\right\}=0.\end{array}

Similarly, Eq. (26) for the nodal internal flux can also be expressed, in a partitioned form, as

(28)\left\{\begin{array}{l}{\mathit{Q}}^{hu}\\ {\mathit{Q}}^{hc}\end{array}\right\}=\xi \left[\begin{array}{ll}{\mathit{E}}_{\mathrm{0}}^{uu}& {\mathit{E}}_{\mathrm{0}}^{uc}\\ {({\mathit{E}}_{\mathrm{0}}^{uc})}^{\mathrm{T}}& {\mathit{E}}_{\mathrm{0}}^{cc}\end{array}\right]\left\{\begin{array}{l}{\mathit{U}}_{,\xi }^{hu}\\ {\mathit{U}}_{,\xi }^{hc}\end{array}\right\}+\left[\begin{array}{ll}{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}& {({\mathit{E}}_{\mathrm{1}}^{cu})}^{\mathrm{T}}\\ {({\mathit{E}}_{\mathrm{1}}^{uc})}^{\mathrm{T}}& {({\mathit{E}}_{\mathrm{1}}^{cc})}^{\mathrm{T}}\end{array}\right]\left\{\begin{array}{l}{\mathit{U}}^{hu}\\ {\mathit{U}}^{hc}\end{array}\right\}.

Eq. (27) can be separated into two systems:

(29){\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}^{uu}{\mathit{U}}_{,\xi \xi }^{hu}+\xi [{\mathit{E}}_{\mathrm{0}}^{uu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{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},

(30)\xi {\mathit{F}}^{tc}=-{\xi }^{\mathrm{2}}{\mathit{F}}^{bc}-{\mathit{F}}^{suc}-{\mathit{F}}^{scc},

where the vectors{\mathit{F}}^{suu},{\mathit{F}}^{suc}, and {\mathit{F}}^{scc}are defined by

(31){\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})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{uc}){\mathit{U}}_{,\xi }^{hc}-{\mathit{E}}_{\mathrm{2}}^{uc}{\mathit{U}}^{hc},

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

(33){\mathit{F}}^{scc}={\xi }^{\mathrm{2}}{\mathit{E}}_{\mathrm{0}}^{cc}{\mathit{U}}_{,\xi \xi }^{hc}+\xi [{\mathit{E}}_{\mathrm{0}}^{cc}+{({\mathit{E}}_{\mathrm{1}}^{cc})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{cc}]{\mathit{U}}_{,\xi }^{hc}-{\mathit{E}}_{\mathrm{2}}^{cc}{\mathit{U}}^{hc}.

By following the same procedure, the partitioned Eq. (28) can also be separated into two systems:

(34){\mathit{Q}}^{hu}(\xi )=\xi {\mathit{E}}_{\mathrm{0}}^{uu}{\mathit{U}}_{,\xi }^{hu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}{\mathit{U}}^{hu}+{\mathit{Q}}^{huc}(\xi ),

(35){\mathit{Q}}^{hc}(\xi )=\xi {({\mathit{E}}_{\mathrm{0}}^{uc})}^{\mathrm{T}}{\mathit{U}}_{,\xi }^{hu}+{({\mathit{E}}_{\mathrm{1}}^{uc})}^{\mathrm{T}}{\mathit{U}}^{hu}+{\mathit{Q}}^{hcc}(\xi ).

where the known vectors {\mathit{Q}}^{huc}(\xi )and{\mathit{Q}}^{hcc}(\xi ) are defined by

(36){\mathit{Q}}^{huc}(\xi )=\xi {\mathit{E}}_{\mathrm{0}}^{uc}{\mathit{U}}_{,\xi }^{hc}+{({\mathit{E}}_{\mathrm{1}}^{uc})}^{\mathrm{T}}{\mathit{U}}^{hc},

(37){\mathit{Q}}^{hcc}(\xi )=\xi {\mathit{E}}_{\mathrm{0}}^{cc}{\mathit{U}}_{,\xi }^{hc}+{({\mathit{E}}_{\mathrm{1}}^{cc})}^{\mathrm{T}}{\mathit{U}}^{hc}.

Now, a system of differential Eq. (29) along with the following two boundary conditions on the inner and outer boundaries:

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

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

A homogeneous solution of the system of linear differential Eq. (29), denoted by{\mathit{U}}_{\mathrm{0}}^{hu}, is derived following standard procedure from the theory of differential equations. The homogeneous solution {\mathit{U}}_{\mathrm{0}}^{hu} must satisfy

(40){\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})}^{\mathrm{T}}-{\mathit{E}}_{\mathrm{1}}^{uu}]{\mathit{U}}_{\mathrm{0},\xi }^{hu}-{\mathit{E}}_{\mathrm{2}}^{uu}{\mathit{U}}_{\mathrm{0}}^{hu}=0,

and the corresponding nodal internal flux , denoted by {\mathit{Q}}_{\mathrm{0}}^{hu}(\xi ), is given by

(41){\mathit{Q}}_{\mathrm{0}}^{hu}(\xi )=\xi {\mathit{E}}_{\mathrm{0}}^{uu}{\mathit{U}}_{\mathrm{0},\xi }^{hu}+{({\mathit{E}}_{\mathrm{1}}^{uu})}^{\mathrm{T}}{\mathit{U}}_{\mathrm{0}}^{hu}.

Since Eq. (40) is a set of (m-p)\mathrm{2} linear, second-order, Euler-Cauchy differential equations, the solution {\mathit{U}}_{\mathrm{0}}^{hu} takes the following form

(42){\mathit{U}}_{\mathrm{0}}^{hu}(\xi )={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}(m-p)\mathrm{2}}{c}_i{\xi }^{{\lambda }_i}{\mathit{\psi }}_i^u},

where a constant {\lambda }_i is termed the modal scaling factor, {\mathit{\psi }}_i is the (m-p)\mathrm{2}-component vector representing the i^th mode of the state variable, and c_iare arbitrary constants denoting the contribution of each mode to the solution. By substituting Eq. (42) into Eq. (40) and Eq. (41), then introducing a \mathrm{2}(m-p)\mathrm{2}-component vector {\mathit{X}}_i such that{\mathit{X}}_i={\left\{\begin{array}{ll}{\mathit{\psi }}_i^u& {\mathit{q}}_i^u\end{array}\right\}}^T, Eqs. (40) and (41) can be combined into a system of linear algebraic equations

(43)\mathit{H}{\mathit{X}}_i={\lambda }_i{\mathit{X}}_i.

where the matrix \mathit{H} is given by

(44)\mathit{H}=\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].

Determination of all \mathrm{2}(m-p)\mathrm{2} pairs \{{\lambda }_i,{\mathit{X}}_i\} is achieved by solving the eigenvalue problem Eq. (41) where {\lambda }_i denote the eigenvalues and {\mathit{X}}_i are associated eigenvectors. In fact, only a half of the eigenvalues has the positive real part whereas the other half has negative real part. Let {\lambda }^{+}and {\lambda }^{-} be(m-p)\mathrm{2}\times (m-p)\mathrm{2}diagonal matrices containing eigenvalues with the positive real part and the negative real part, respectively. Also, let {\mathit{\Phi }}^{\psi +} and {\mathit{\Phi }}^{q+} be matrices whose columns containing, respectively, all vectors {\mathit{\psi }}_i^u and {\mathit{q}}_i^u obtained from the eigenvectors{\mathit{X}}_i={\left\{\begin{array}{ll}{\mathit{\psi }}_i^u& {\mathit{q}}_i^u\end{array}\right\}}^{\mathrm{T}} associated with all eigenvalues contained in {\mathit{\lambda }}^{+} and let {\mathit{\Phi }}^{\psi -} and {\mathit{\Phi }}^{q-} be matrices whose columns containing, respectively, all vectors {\mathit{\psi }}_i^u and {\mathit{q}}_i^u obtained from the eigenvectors {\mathit{X}}_i={\left\{\begin{array}{ll}{\mathit{\psi }}_i^u& {\mathit{q}}_i^u\end{array}\right\}}^{\mathrm{T}} associated with all eigenvalues contained in{\lambda }^{-}. Now, the homogeneous solutions {\mathit{U}}_{\mathrm{0}}^{hu} and {\mathit{Q}}_{\mathrm{0}}^{hu}(\xi ) are given by

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

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

where {\mathit{\Pi }}^{+} and {\mathit{\Pi }}^{-} are diagonal matrices obtained by simply replacing the diagonal entries {\lambda }_i of the matrices {\lambda }^{+}and {\lambda }^{-}by the a function {\xi }^{{\lambda }_i}, respectively; and {\mathit{C}}^{+}and {\mathit{C}}^{-}are vectors containing arbitrary constants representing the contribution of each mode. It is apparent that the diagonal entries of {\mathit{\Pi }}^{+}become infinite when \xi \to \infty whereas those of {\mathit{\Pi }}^{-}is unbounded when \xi \to \mathrm{0}. As a result, {\mathit{C}}^{+}is taken to 0 to ensure the boundedness of the solution for unbounded bodies and, similarly, the condition {\mathit{C}}^{-}=0is enforced for bodies containing the scaling center.

A particular solution of Eq. (29), denoted by{\mathit{U}}_{\mathrm{1}}^{hu}, associated with the distributed body source, the surface flux on the side face and the prescribed state variable on the side face can also be obtained from a standard procedure in the theory of differential equations such as the method of undetermined coefficient. Once the particular solution{\mathit{U}}_{\mathrm{1}}^{hu}is obtained, the corresponding particular nodal internal flux {\mathit{Q}}_{\mathrm{1}}^{hu} can be calculated. Finally, the general solution of Eq. (29) and the corresponding nodal internal flux are then given by

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

(48){\mathit{Q}}^{hu}(\xi )={\mathit{\Phi }}^{q+}{\mathit{\Pi }}^{+}(\xi ){\mathit{C}}^{+}+{\mathit{\Phi }}^{q-}{\mathit{\Pi }}^{-}(\xi ){\mathit{C}}^{-}+{\mathit{Q}}_{\mathrm{1}}^{hu}(\xi ).

To determine the constants contained in {\mathit{C}}^{+}and{\mathit{C}}^{-}, the boundary conditions on both inner and outer boundaries are enforced. By enforcing the conditions Eqs. (38)‒(39), it gives rise to

(49)\left\{\begin{array}{l}{\mathit{C}}^{+}\\ {\mathit{C}}^{-}\end{array}\right\}={\left[\begin{array}{ll}{\mathit{\Phi }}^{q+}{\mathit{\Pi }}^{+}({\xi }_{\mathrm{1}})& {\mathit{\Phi }}^{q-}{\mathit{\Pi }}^{-}({\xi }_{\mathrm{1}})\\ {\mathit{\Phi }}^{q+}{\mathit{\Pi }}^{+}({\xi }_{\mathrm{2}})& {\mathit{\Phi }}^{q-}{\mathit{\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).

From Eq. (49), it can readily be obtained and substituting Eq. (45) into its yields

(50)\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\mathit{K}, commonly termed the stiffness matrix, is given by

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

Then applying the prescribed surface flux and the state variable on both inner and outer boundaries, a system of linear algebraic Eq. (29) is sufficient for determining all involved unknowns. Once the unknowns on both the inner and outer boundaries are solved, the approximate field quantities such as the state variable and the surface flux within the body can readily be post-processed.

(52){\mathit{u}}^h(\xi ,s)={\mathit{N}}^{Su}(s){\mathit{U}}^{hu}(\xi )+{\mathit{N}}^{Sc}(s){\mathit{U}}^{hc}(\xi ),

where{\mathit{N}}^{Su}and {\mathit{N}}^{Sc}are matrices resulting from the partition of {\mathit{N}}^S. Similarly, the approximated body flux can be computed from Eq. (14) as

(53){\sigma }^h(\xi ,s)=\mathit{D}\left[{\mathit{B}}_{\mathrm{1}}^u(s){\mathit{U}}_{,\xi }^{hu}(\xi )+{\displaystyle \frac{\mathrm{1}}{\xi }}{\mathit{B}}_{\mathrm{2}}^u(s){\mathit{U}}^{hu}(\xi )\right]+\mathit{D}\left[{\mathit{B}}_{\mathrm{1}}^c(s){\mathit{U}}_{,\xi }^{hc}(\xi )+{\displaystyle \frac{\mathrm{1}}{\xi }}{\mathit{B}}_{\mathrm{2}}^c(s){\mathit{U}}^{hc}(\xi )\right],

where{\mathit{B}}_{\mathrm{1}}^u,{\mathit{B}}_{\mathrm{1}}^cand{\mathit{B}}_{\mathrm{2}}^u,{\mathit{B}}_{\mathrm{2}}^c are matrices resulting from the partition of the matrices B₁ and B₂, respectively. It is emphasized here again that the solutions Eq. (52) and Eq. (53) also apply to the special cases of bounded and unbounded bodies. For bounded bodies containing the scaling center, {\mathit{C}}^{-}simply vanishes and, for unbounded bodies,{\mathit{C}}^{+}=\mathrm{0}. To investigate the error of the scaled boundary finite element approximation, an error vector-valued functioner=er(x)is first defined by

(54)\mathit{e}r(x)={\mathit{u}}^{\mathrm{exact}}(\mathit{x})-{\mathit{u}}^h(\mathit{x}),

where {\mathit{u}}^{\mathrm{exact}}and {\mathit{u}}^hdenotes the exact solution and scaled boundary finite element solution of the state variables, respectively. The following standardL^{\mathrm{2}}-norm is employed to measure the magnitude of the error function, i.e.,

(55){\Vert \mathit{e}r\Vert }_{L^{\mathrm{2}}}=\sqrt{\left({\displaystyle \underset{\mathrm{\Omega }}{\int }\mathit{e}{r}^{\mathrm{T}}(x)\mathit{e}r(x)}\right)}.

Note that the evaluation of the above integral can be achieved efficiently in the\xi -sspace. Finally, the relative prevent error, denoted byerror, is obtained by normalizing the L²-norm of the error function by the L²-norm of the exact state variable, i.e.,

(56)error={\displaystyle \frac{{\Vert \mathit{e}r\Vert }_{L^{\mathrm{2}}}}{{\Vert {\mathit{u}}^{\mathrm{exact}}\Vert }_{L^{\mathrm{2}}}}}\times \mathrm{100},

where

{\Vert {\mathit{u}}^{\mathrm{exact}}\Vert }_{L^{\mathrm{2}}}=\sqrt{\left({\displaystyle \underset{\Omega }{\int }{[{\mathit{u}}^{\mathrm{exact}}(x)]}^{\mathrm{T}}{\mathit{u}}^{\mathrm{exact}}(x)}\right)}.

It is noticed that for problems without the exact solution, the converged solution obtained from a particular, sufficiently fine mesh can be used to estimate the error for any level of meshes used in the approximation.

Performance of the method

Some numerical examples to verify the proposed technique and demonstrate its performance and capabilities. To demonstrate its capability to treat a variety of boundary value problems, general boundary conditions, and prescribed data on the side faces, the types of problems associated with linear elasticity under various scenarios are considered. The conventional polar coordinates are exploited to achieve the exact description of the circular defining curve, exact geometry of domain, and the standard finite element shape functions are employed to discretize both the defining curve and the trial and test functions. The accuracy and convergence of numerical solutions are confirmed by benchmarking with available analytical solutions and carrying out the analysis via a series of meshes.

Pressurized circular hole in an infinite domain

Consider a pressurized circular hole of radius R₁ in an infinite domain as shown in Fig. 4(a). The medium is made of a homogeneous, linearly elastic, isotropic material with Young’s modulus E and Poisson’s ratio \nu and is subjected to a plane-strain condition and uniform pressure p₁ on the boundary, respectively. Due to the symmetry, it is sufficient to model this problem using only a quarter of the cylinder (see Fig. 4(b)) with appropriate conditions on both side faces (i.e., the normal displacement and tangential traction on the side faces vanish).

To describe the geometry, the scaling center is chosen at the center of the cylinder whereas the inner boundary is treated as the defining curve. In a numerical study, the Poisson’s ratio \nu =\mathrm{0.3} and meshes with N identical linear elements are employed. Both exact geometry and approximation geometry are investigated in the modelling. Results for the relative error norm versus number of elements are reported along with existing analytical solutions [²⁵] in Fig. 5.

As can be seen exact defining curve exhibits more accuracy compared to approximate defining cure as the same number of elements is employed for both cases. Results for the normalized radial displacement (u_r/(p_{\mathrm{1}}{R}_{\mathrm{1}}/E)), normalized radial stress ({\sigma }_{rr}/p_{\mathrm{1}}) and normalized hoop stress ({\sigma }_{\theta \theta }/p_{\mathrm{1}}) are reported along with existing analytical solutions [²⁵] in Figs. 6, 7, respectively, for four meshes of exact defining curve. It is seen that numerical solutions generated by the proposed technique converge and exhibit excellent agreement with the benchmark solution. It is worth noting that the discretization with only few linear elements can capture numerical solution with the sufficient accuracy.

Thick circular disk

Another problem, consider a thick circular disk of radius R as shown in Fig. 8(a). The medium is made of a homogeneous, linearly elastic, isotropic material with Young’s modulus E and Poisson’s ratio v and subjected to radial body force b_r=\rho w{r}^{\mathrm{2}}whereas its boundary is restrained against the movement in the radial direction. The exact solution for this particular under plane strain condition can be readily obtained from the classical theory of linear elasticity.

Due to the symmetry, it is sufficient to model this problem using only a quarter of the circular disk (see Fig. 8(b)) with appropriate condition on side face (i.e., the normal displacement and tangential traction on the side faces vanish). To describe the geometry, the scaling centre is chosen at the centre of the circular disk whereas the boundary is treated as the defining curve. In a numerical study, the Poisson’s ratio \nu =\mathrm{0.3} and meshes with Nidentical linear elements are employed both exact and approximation defining curve. The rate of convergence of the approximation is investigated with linear elements in both cases of defining curve. The rate of convergence of the approximation is investigated with linear elements in both cases of defining curve. Plots of the relative percent errors versus the number of elements are reported in the Fig. 9. It is seen that the rate of convergence from two cases of defining curves is similar. However, if the number of elements at the same mesh, a significant increase in accuracy can be obtained higher by using exact defining curve. Resulting from the normalized radial displacement{\underline{\mathrm{u}}}_r=(\mathrm{15}\mu (\kappa +\mathrm{1}){\mathit{u}}_r/((\kappa -\mathrm{1})\rho w{R}^{\mathrm{4}})),normalized radial stress {\underset{\_}{\sigma }}_{rr}=(\mathrm{15}(\kappa +\mathrm{1}){\sigma }_{rr}/(\rho w{R}^{\mathrm{3}})), and normalized hoop stress {\underset{\_}{\sigma }}_{\theta \theta }=(\mathrm{15}(\kappa +\mathrm{1}){\sigma }_{\theta \theta }/(\rho w{R}^{\mathrm{3}}))(exact defining curve) are reported with exact solutions in Figs. 10, 11. It can be seen that numerical solutions converges to the exact solutions as the number of elements N used to discretize in the defining curve increase and, in addition, only few number of degrees of freedom is sufficient to obtain accurate displacements. Similar to displacements, the present method also yields highly accurate stress components and the good behaviour.

Concentrated force in an infinite domain

As the last example, a representative boundary value problem associated with an infinite domain subjected a concentrated force P to investigate the capability of the proposed technique as shown in Fig. 12(a). The medium is made of a homogeneous, linearly elastic, isotropic material with Young’s modulus E and Poisson’s ratio\nu, modulus matrix D for the plane strain condition. With concentrated forcePis applied in an infinite domain, the analytical solution is available [²⁶]. The external load at radial r₀ of an infinite domain is considered as the traction on surface matching the circular hole in an infinite domain

(57)p_{x_{\mathrm{1}}}={\displaystyle \frac{P}{\mathrm{2}\pi (\mathrm{1}-v)r_{\mathrm{0}}}}\cos (\theta )\sin (\varphi )

(58)p_{x_{\mathrm{2}}}={\displaystyle \frac{P}{\mathrm{4}\pi (\mathrm{1}-v)r_{\mathrm{0}}}}((\mathrm{1}-\mathrm{2}v)+\mathrm{2}{\sin }^{\mathrm{2}}(\varphi )).

In the geometry modeling, the scaling center is chosen at the center of an infinite domain whereas the r_{\mathrm{0}} boundary is treated as defining curve (see Fig. 12(b)). In the analysis, the Poisson’s ratio is taken as \nu =\mathrm{0.2} and defining curve is discretized byNidentical linear elements. The exact geometry is used in defining curve and linear elements are approximated in displacement fields. The normalized stress components{\sigma }_{\mathrm{11}}/(P/r_{\mathrm{0}})and {\sigma }_{\mathrm{22}}/(P/r_{\mathrm{0}})are reported along radial (angle θ=90°) in Fig. 13 and Fig. 14. It can be seen that the discretization with only few linear elements can capture numerical solution to exact solution with the sufficient accuracy.

Conclusions

A numerical technique based on the scaled boundary finite element method has been successfully developed for solving two-dimensional elasticity, boundary value problems with the domain completely described by a circular defining curve. Both the formulation and implementations have been established in an associated data (such as the domain geometry, the prescribed distributed body source, boundary conditions, and contribution of the side face) to be treated in a single, unified manner. Results from several numerical study have indicated that the proposed SBFEM yields highly accurate numerical solutions with the percent error weakly dependent on the level of mesh refinement. The results also show that it is advantageous to use circular defining curve, and that higher convergence can be obtained.

References

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

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

[3] Wolf J P, Song C. Finite-element modelling of undounded media. In: Proceedings of Eleventh World Conference on Earthquake Engineering, 1996: Paper No. 70

[4] 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

[5] 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

[6] 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

[7] 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

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

[9] 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

[10] 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

[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, 2005, 32(6): 436–444

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

[13] 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

[14] 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

[15] 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

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

[17] He Y, Yang H, Xu M, Deeks A J. A scaled boundary finite element method for cyclically symmetric two-dimensional elastic analysis. Computers & Structures, 2013, 120: 1–8

[18] 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

[19] 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

[20] Chan C L, Anitescu C, Rabczuk T. Volumetric parametrization from a level set boundary representation with PHT-splines. Computer Aided Design, 2017, 82: 29–41

[21] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: an overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[22] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

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

[24] 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

[25] Sadd M H. Elasticity: Theory, Application, and Numerics. Elsevier Academic Press, 2005

[26] Karasudhi P. Foundation of Solid Mechanics. Kluwer Academic Publishers, 1991

RIGHTS & PERMISSIONS

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