Table of Contents Author Guidelines Submit a Manuscript
Advances in Acoustics and Vibration
Volume 2011, Article ID 973591, 11 pages
http://dx.doi.org/10.1155/2011/973591
Research Article

Dynamic Analysis of Wind Turbine Blades Using Radial Basis Functions

Department of Electrical Engineering, National Penghu University of Science and Technology, Penghu 880, Taiwan

Received 30 October 2010; Revised 11 April 2011; Accepted 11 April 2011

Academic Editor: K. M. Liew

Copyright © 2011 Ming-Hung Hsu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Wind turbine blades play important roles in wind energy generation. The dynamic problems associated with wind turbine blades are formulated using radial basis functions. The radial basis function procedure is used to transform partial differential equations, which represent the dynamic behavior of wind turbine blades, into a discrete eigenvalue problem. Numerical results demonstrate that rotational speed significantly impacts the first frequency of a wind turbine blade. Moreover, the pitch angle does not markedly affect wind turbine blade frequencies. This work examines the radial basis functions for dynamic problems of wind turbine blade.

1. Introduction

Wind energy will likely play an important role as a future energy source in countries worldwide. Energy is essential to economic growth. In response to environmental concerns, the amount of energy generated using renewable resources is increasing [13]. Herbert et al. [2] reviewed assessment models for wind resources, site selection models, and aerodynamic models, including the wake effect model. Varol et al. [3] demonstrated the positive effect of aerofoils surrounding wind blades. Their notion of surrounding aerofoils resembles that of steering blades in a water turbine. Spee et al. [4] created a novel control strategy for a brushless doubly fed machine applied to variable-speed wind-power generation systems. Vitale and Rossi [5] designed low-power, horizontal-axis wind turbine blades, using an iterative algorithm. With their software, one can determine the optimum blade shape for a wind turbine to satisfy the energy requirements of an electrical system with optimum rotor efficiency. Song and Dhinakaran [6] developed a nonlinear wind turbine control method. Storti and Aboelnaga [7] analyzed transverse deflections of a straight tapered symmetrical beam attached to a rotating hub as a model for bending vibration of blades in turbomachinery. The natural frequencies of a single tapered and pretwisted turboblade were calculated by Rao [8, 9], Abrate [10], Hodges et al. [11], and Dawson and Carneige [12, 13] using the Rayleigh-Ritz scheme. Gupta and Rao [14] applied the finite element method to identify the frequencies of natural vibration of doubly tapered and twisted beams. Swaminathan and Rao [15] described the vibrations of a rotating pretwisted and tapered blade. Subrahmanyam et al. [16, 17] derived the natural frequencies and mode shapes of a uniform pretwisted cantilever blade using the Reissner approach. Chen and Keer [18] examined the transverse vibrations of a rotating pretwisted Timoshenko beam under axial loading. Storti and Aboelnaga [19] examined transverse deflections of a straight tapered symmetrical beam attached to a rotating hub as a model for the bending vibration of blades. Wagner [20] derived the forced vibration response of subsystems with different natural frequencies and damping attached to a foundation with a finite stiffness and mass. Griffin and Sinha [2123] determined the dynamic responses of frictionally damped turbine blades. Determining the dynamic features of wind turbine blades is of significant importance. Radial basis functions are important elements of approaches generally called meshless methods. In this study, dynamic problems associated with wind turbine blades are formulated using radial basis functions.

2. Radial Basis Function

A radial basis function is a real-value function whose value depends on distance from an origin. Kansa [24, 25] investigated a given function or partial derivatives of a function with respect to a coordinate direction expressed as a linear weighted sum of all functional values at all mesh points along the direction initiated based on the concept of the radial basis function. In their algorithm, node distribution was entirely unstructured. Wang and Liu [26] developed a point interpolation meshless method based on radial basis functions that incorporated the Galerkin weak form for solving partial differential equations. Elfelsoufi [27] investigated buckling, flutter, and vibration of beams using radial basis functions. Hon et al. [28] utilized radial basis functions for function fitting and solving partial differential equations using global nodes and collocation procedures. Liu et al. [29] constructed shape functions with the delta function property using radial and polynomial basis functions. Devi and Pepper [30] generated a simplified radial basis function approach for calculating coupled heat transfer of a fluid flow using a local pressure correction scheme. In their study, shape functions were constructed using radial basis functions. Vrankar et al. [31] applied a relatively new approach for modeling radionuclide migration through the geosphere using a radial basis function scheme. Qiao and Ernst [32] applied a nonlinear approach for constructing color conversions based on radial basis functions. Their experimental results demonstrated that color conversion can be achieved effectively using the radial basis function approach when constructing nonlinear data maps. A radial basis function can be expressed as follows [33, 34]: 𝐵𝑖(𝑟)=𝑟𝑟𝑖2+𝑐2,(1) where 𝑐 is a shape parameter. The radial basis function is generally utilized to develop functional approximations in the following forms [33, 34]:𝑢𝑘(𝑟,𝑡)=𝑁𝑖=1𝑎𝑘𝑖(𝑡)𝐵𝑖𝑣(𝑟),for𝑘=1,2,3,(2)𝑘(𝑟,𝑡)=𝑁𝑖=1𝑎𝑘𝑖(𝑡)𝐵𝑖𝑈(𝑟),for𝑘=1,2,3,(3)𝑘(𝑟)=𝑁𝑖=1𝑏𝑘𝑖𝐵𝑖𝑉(𝑟),for𝑘=1,2,3,(4)𝑘(𝑟)=𝑁𝑖=1𝑏𝑘𝑖𝐵𝑖(𝑟),for𝑘=1,2,3,(5) where 𝑎𝑘𝑖, 𝑎𝑘𝑖, 𝑏𝑘𝑖, and 𝑏𝑘𝑖 are coefficients to be determined. Blade deflection 𝑢𝑘(𝑟,𝑡) is the sum of 𝑁 radial basis functions, each associated with a different center 𝑟𝑖. Blade deflection 𝑣𝑘(𝑟,𝑡) denotes the sum of 𝑁 radial basis functions, each associated with a different center 𝑟𝑖. Blade deflection 𝑈𝑘(𝑟) is the sum of 𝑁 radial basis functions, each associated with a different center 𝑟𝑖. Blade deflection 𝑉𝑘(𝑟) denotes the sum of 𝑁 radial basis functions, each associated with a different center 𝑟𝑖. The domain contains 𝑁 collocation points.

3. Dynamic Analysis of Wind Turbine Blades

The kinetic energy of rotating wind turbine blades can be derived as𝑇𝑒=3𝑘=112𝐿0𝜌𝐴𝜕𝑢𝑘𝜕𝑡2+𝜕𝑣𝑘𝜕𝑡2+Ω𝑣𝑘cos𝜃2+Ω𝑢𝑘sin𝜃2𝑑𝑟,(6) where Ω is hub rotational speed, 𝐴 is cross-sectional area, 𝜃 is pitch angle, 𝜌 is blade density, and 𝐿 is blade length. The corresponding strain energy of rotating blades is𝑈𝑒=3𝑘=112𝐿0𝐸𝐼𝑦𝑦𝜕𝑢𝑘𝜕𝑟2+2𝐼𝑥𝑦𝜕𝑢𝑘𝜕𝑟𝜕𝑣𝑘𝜕𝑟+𝐼𝑥𝑥𝜕𝑣𝑘𝜕𝑟2+1𝑑𝑟2𝐿0𝐿𝑟𝜌𝐴Ω2𝑟+𝑟0×𝑑𝑟𝜕𝑢𝑘𝜕𝑟2+𝜕𝑣𝑘𝜕𝑟2,𝑑𝑟(7) where 𝑟0 is hub radius, 𝑢1(𝑟,𝑡) is the displacement of the first blade on the 𝑥-axis, 𝑣1(𝑟,𝑡) is the displacement of the first blade on the 𝑦-axis, 𝑢2(𝑟,𝑡) is the displacement of the second blade on the 𝑥-axis, 𝑣2(𝑟,𝑡) is the displacement of the second blade on the 𝑦-axis, 𝑢3(𝑟,𝑡) is the displacement of the third blade in the 𝑥-axis, 𝑣3(𝑟,𝑡) is the displacement of the third blade on the 𝑦-axis, and 𝐸 is Young’s modulus. By examining the internal and external damping effects of blades, virtual work 𝛿𝑊𝑒 in wind turbine blades can be derived as𝛿𝑊𝑒=3𝑘=1𝐿0𝐶0𝜕𝑢𝑘𝜕𝑡𝛿𝑢𝑘𝑑𝑟𝐿0𝐶1𝜕2𝐸𝐼𝑦𝑦𝜕𝑟2𝜕3𝑢𝑘𝜕𝑡𝜕𝑟2𝛿𝑢𝑘𝑑𝑟2𝐿0𝐶1𝜕𝐸𝐼𝑦𝑦𝜕𝜕𝑟4𝑢𝑘𝜕𝑡𝜕𝑟3𝛿𝑢𝑘𝑑𝑟𝐿0𝐶1𝐸𝐼𝑦𝑦𝜕5𝑢𝑘𝜕𝑡𝜕𝑟4𝛿𝑢𝑘𝑑𝑟𝐿0𝐶0𝜕𝑣𝑘𝜕𝑡𝛿𝑣𝑘𝑑𝑟𝐿0𝐶1𝜕2𝐸𝐼𝑥𝑥𝜕𝑟2𝜕3𝑣𝑘𝜕𝑡𝜕𝑟2𝛿𝑣𝑘𝑑𝑟2𝐿0𝐶1𝜕𝐸𝐼𝑥𝑥𝜕𝜕𝑟4𝑣𝑘𝜕𝑡𝜕𝑟3𝛿𝑣𝑘𝑑𝑟𝐿0𝐶1𝐸𝐼𝑥𝑥𝜕5𝑣𝑘𝜕𝑡𝜕𝑟4𝛿𝑣𝑘,𝑑𝑟(8) where 𝐶0 and 𝐶1 are external and internal damping coefficients of blades, respectively. Equations (6)–(8) are integrated into the Hamilton equation as follows:𝑡2𝑡1(𝛿𝑇𝑒𝛿𝑈𝑒+𝛿𝑊𝑒)𝑑𝑡=0.(9) This leads to the following equations for the motion of the blades:𝜕2𝐸𝐼𝑦𝑦𝜕𝑟2𝜕2𝑢𝑘𝜕𝑟2+2𝜕𝐸𝐼𝑦𝑦𝜕𝜕𝑟3𝑢𝑘𝜕𝑟3+𝐸𝐼𝑦𝑦𝜕4𝑢𝑘𝜕𝑟4+𝜕2𝐸𝐼𝑥𝑦𝜕𝑟2𝜕2𝑣𝑘𝜕𝑟2+2𝜕𝐸𝐼𝑥𝑦𝜕𝜕𝑟3𝑣𝑘𝜕𝑟3+𝐸𝐼𝑥𝑦𝜕4𝑣𝑘𝜕𝑟4𝜌𝐴Ω2𝑢𝑘sin2𝜃𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝐴𝑟+𝑟0𝑑𝑟𝜕𝑢𝑘𝜕𝑟𝜌Ω2𝐿𝑟𝐴𝑟+𝑟0𝜕𝑑𝑟2𝑢𝑘𝜕𝑟2+𝐶0𝜕𝑢𝑘𝜕𝑡+𝐶1𝜕2𝐸𝐼𝑦𝑦𝜕𝑟2𝜕3𝑢𝑘𝜕𝑟2𝜕𝑡+2𝐶1𝜕𝐸𝐼𝑦𝑦𝜕𝜕𝑟4𝑢𝑘𝜕𝑟3𝜕𝑡+𝐶1𝐸𝐼𝑦𝑦𝜕5𝑢𝑘𝜕𝑟4𝜕𝑡+𝐶1𝜕2𝐸𝐼𝑥𝑦𝜕𝑟2𝜕3𝑣𝑘𝜕𝑟2𝜕𝑡+2𝐶1𝜕𝐸𝐼𝑥𝑦𝜕𝜕𝑟4𝑣𝑘𝜕𝑟3𝜕𝑡+𝐶1𝐸𝐼𝑥𝑦𝜕5𝑣𝑘𝜕𝑟4𝜕𝜕𝑡+𝜌𝐴2𝑢𝑘𝜕𝑡2𝜕=0for𝑘=1,2,3,2𝐸𝐼𝑥𝑥𝜕𝑟2𝜕2𝑣𝑘𝜕𝑟2+2𝜕𝐸𝐼𝑥𝑥𝜕𝜕𝑟3𝑣𝑘𝜕𝑟3+𝐸𝐼𝑥𝑥𝜕4𝑣𝑘𝜕𝑟4+𝜕2𝐸𝐼𝑥𝑦𝜕𝑟2𝜕2𝑢𝑘𝜕𝑟2+2𝜕𝐸𝐼𝑥𝑦𝜕𝜕𝑟3𝑢𝑘𝜕𝑟3+𝐸𝐼𝑥𝑦𝜕4𝑢𝑘𝜕𝑟4𝜌𝐴Ω2𝑣𝑘cos2𝜃𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝐴𝑟+𝑟0𝑑𝑟𝜕𝑣𝑘𝜕𝑟𝜌Ω2𝐿𝑟𝐴𝑟+𝑟0𝜕𝑑𝑟2𝑣𝑘𝜕𝑟2+𝐶0𝜕𝑣𝑘𝜕𝑡+𝐶1𝜕2𝐸𝐼𝑥𝑥𝜕𝑟2𝜕3𝑣𝑘𝜕𝑟2𝜕𝑡+2𝐶1𝜕𝐸𝐼𝑥𝑥𝜕𝜕𝑟4𝑣𝑘𝜕𝑟3𝜕𝑡+𝐶1𝐸𝐼𝑥𝑥𝜕5𝑣𝑘𝜕𝑟4𝜕𝑡+𝐶1𝜕2𝐸𝐼𝑥𝑦𝜕𝑟2𝜕3𝑢𝑘𝜕𝑟2𝜕𝑡+2𝐶1𝜕𝐸𝐼𝑥𝑦𝜕𝜕𝑟4𝑢𝑘𝜕𝑟3𝜕𝑡+𝐶1𝐸𝐼𝑥𝑦𝜕5𝑢𝑘𝜕𝑟4𝜕𝜕𝑡+𝜌𝐴2𝑣𝑘𝜕𝑡2=0for𝑘=1,2,3.(10) The following equations are the corresponding boundary conditions:𝑢𝑘(0,𝑡)=0for𝑘=1,2,3,𝜕𝑢𝑘(0,𝑡)𝜕𝑟=0for𝑘=1,2,3,𝐸𝐼𝑦𝑦𝜕2𝑢𝑘(𝐿,𝑡)𝜕𝑟2𝜕=0for𝑘=1,2,3,𝜕𝑟𝐸𝐼𝑦𝑦𝜕2𝑢𝑘(𝐿,𝑡)𝜕𝑟2𝑣=0for𝑘=1,2,3,𝑘(0,𝑡)=0for𝑘=1,2,3,𝜕𝑣𝑘(0,𝑡)𝜕𝑟=0for𝑘=1,2,3,𝐸𝐼𝑥𝑥𝜕2𝑣𝑘(𝐿,𝑡)𝜕𝑟2𝜕=0for𝑘=1,2,3,𝜕𝑟𝐸𝐼𝑥𝑥𝜕2𝑣𝑘(𝐿,𝑡)𝜕𝑟2=0for𝑘=1,2,3.(11) A solution is thus assumed as𝑢𝑘=𝑈𝑘(𝑟)𝑒𝜆𝑡𝑣for𝑘=1,2,3,𝑘=𝑉𝑘(𝑟)𝑒𝜆𝑡for𝑘=1,2,3.(12) This yields the following equations for the motion of blades: 𝑑2𝐸𝐼𝑦𝑦𝑑𝑟2𝑑2𝑈𝑘𝑑𝑟2+2𝑑𝐸𝐼𝑦𝑦𝑑𝑑𝑟3𝑈𝑘𝑑𝑟3+𝐸𝐼𝑦𝑦𝑑4𝑈𝑘𝑑𝑟4+𝑑2𝐸𝐼𝑥𝑦𝑑𝑟2𝑑2𝑉𝑘𝑑𝑟2+2𝑑𝐸𝐼𝑥𝑦𝑑𝑑𝑟3𝑉𝑘𝑑𝑟3+𝐸𝐼𝑥𝑦𝑑4𝑉𝑘𝑑𝑟4𝜌𝐴Ω2𝑈𝑘sin2𝜃𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝐴𝑟+𝑟0𝑑𝑟𝑑𝑈𝑘𝑑𝑟𝜌Ω2𝐿𝑟𝐴𝑟+𝑟0𝑑𝑑𝑟2𝑈𝑘𝑑𝑟2+𝜆𝐶0𝑈𝑘+𝜆𝐶1𝑑2𝐸𝐼𝑦𝑦𝑑𝑟2𝑑2𝑈𝑘𝑑𝑟2+2𝜆𝐶1𝑑𝐸𝐼𝑦𝑦𝑑𝑑𝑟3𝑈𝑘𝑑𝑟3+𝜆𝐶1𝐸𝐼𝑦𝑦𝑑4𝑈𝑘𝑑𝑟4+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑑𝑟2𝑑2𝑉𝑘𝑑𝑟2+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑑𝑑𝑟3𝑉𝑘𝑑𝑟3+𝜆𝐶1𝐸𝐼𝑥𝑦𝑑4𝑉𝑘𝑑𝑟4+𝜆2𝜌𝐴𝑈𝑘𝑑=0for𝑘=1,2,3,2𝐸𝐼𝑥𝑥𝑑𝑟2𝑑2𝑉𝑘𝑑𝑟2+2𝑑𝐸𝐼𝑥𝑥𝑑𝑑𝑟3𝑉𝑘𝑑𝑟3+𝐸𝐼𝑥𝑥𝑑4𝑉𝑘𝑑𝑟4+𝑑2𝐸𝐼𝑥𝑦𝑑𝑟2𝑑2𝑈𝑘𝑑𝑟2+2𝑑𝐸𝐼𝑥𝑦𝑑𝑑𝑟3𝑈𝑘𝑑𝑟3+𝐸𝐼𝑥𝑦𝑑4𝑈𝑘𝑑𝑟4𝜌𝐴Ω2𝑉𝑘cos2𝜃𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝐴𝑟+𝑟0𝑑𝑟𝑑𝑉𝑘𝑑𝑟𝜌Ω2𝐿𝑟𝐴𝑟+𝑟0𝑑𝑑𝑟2𝑉𝑘𝑑𝑟2+𝜆𝐶0𝑑𝑉𝑘𝑑𝑡+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑥𝑑𝑟2𝑑2𝑉𝑘𝑑𝑟2+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑥𝑑𝑑𝑟3𝑉𝑘𝑑𝑟3+𝜆𝐶1𝐸𝐼𝑥𝑥𝑑4𝑉𝑘𝑑𝑟4+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑑𝑟2𝑑2𝑈𝑘𝑑𝑟2+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑑𝑑𝑟3𝑈𝑘𝑑𝑟3+𝜆𝐶1𝐸𝐼𝑥𝑦𝑑4𝑈𝑘𝑑𝑟4+𝜆2𝜌𝐴𝑉𝑘=0for𝑘=1,2,3.(13) The corresponding boundary conditions are as follows:𝑈𝑘(0)=0for𝑘=1,2,3,𝑑𝑈𝑘(0)𝑑𝑟=0for𝑘=1,2,3,𝐸𝐼𝑦𝑦𝑑2𝑈𝑘(𝐿)𝑑𝑟2𝑑=0for𝑘=1,2,3,𝑑𝑟𝐸𝐼𝑦𝑦𝑑2𝑈𝑘(𝐿)𝑑𝑟2𝑉=0for𝑘=1,2,3,𝑘(0)=0for𝑘=1,2,3,𝑑𝑉𝑘(0)𝑑𝑟=0for𝑘=1,2,3,𝐸𝐼𝑥𝑥𝑑2𝑉𝑘(𝐿)𝑑𝑟2𝑑=0for𝑘=1,2,3,𝑑𝑟𝐸𝐼𝑥𝑥𝑑2𝑉𝑘(𝐿)𝑑𝑟2=0for𝑘=1,2,3.(14) By employing the radial basis function approach, (2) and (3) can be substituted into (10). The equations of motion of wind turbine blades can be rearranged in the following matrix form:𝑟𝜌𝐴𝑖𝐵1𝑟𝑖𝑟𝜌𝐴𝑖𝐵2𝑟𝑖𝑟𝜌𝐴𝑖𝐵𝑁𝑟𝑖𝜕2𝑎𝑘1𝜕𝑡2𝜕2𝑎𝑘2𝜕𝑡2𝜕2𝑎𝑘𝑁𝜕𝑡2𝑇+𝐶0𝐵1𝑟𝑖𝐶0𝐵2𝑟𝑖𝐶0𝐵𝑁𝑟𝑖𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+2𝐶1𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+2𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+2𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑦𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐸𝐼𝑦𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+2𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌𝐴(𝑟)Ω2(sin𝜃)2𝐵1𝑟𝑖𝜌𝐴(𝑟)Ω2(sin𝜃)2𝐵2𝑟𝑖𝜌𝐴(𝑟)Ω2(sin𝜃)2𝐵𝑁𝑟𝑖𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵1𝑟𝑖𝜕𝑟𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵2𝑟𝑖𝜕𝑟𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵𝑁𝑟𝑖×𝑎𝜕𝑟𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵1𝑟𝑖𝜕𝑟2𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵2𝑟𝑖𝜕𝑟2𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵𝑁𝑟𝑖𝜕𝑟2×𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0],for𝑟𝑖=1,2,,𝑁,𝑘=1,2,3,𝜌𝐴𝑖𝐵1𝑟𝑖𝑟𝜌𝐴𝑖𝐵2𝑟𝑖𝑟𝜌𝐴𝑖𝐵𝑁𝑟𝑖𝜕2𝑎𝑘1𝜕𝑡2𝜕2𝑎𝑘2𝜕𝑡2𝜕2𝑎𝑘𝑁𝜕𝑡2𝑇+𝐶0𝐵1𝑟𝑖𝐶0𝐵2𝑟𝑖𝐶0𝐵𝑁𝑟𝑖𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+2𝐶1𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝜕𝑎𝑘1𝜕𝜕𝑡𝑎𝑘2𝜕𝜕𝑡𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝐶1𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+2𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝐶1𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝜕𝑎𝑘1𝜕𝑡𝜕𝑎𝑘2𝜕𝑡𝜕𝑎𝑘𝑁𝜕𝑡𝑇+𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+2𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑥𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑥𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵1𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵2𝑟𝑖𝜕𝑟2𝜕2𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝑟2𝜕2𝐵𝑁𝑟𝑖𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+2𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵1𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵2𝑟𝑖𝜕𝑟32𝜕𝐸𝐼𝑥𝑦𝑟𝑖𝜕𝜕𝑟3𝐵𝑁𝑟𝑖𝜕𝑟3𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵1𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵2𝑟𝑖𝜕𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝜕4𝐵𝑁𝑟𝑖𝜕𝑟4𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌𝐴(𝑟)Ω2(cos𝜃)2𝐵1𝑟𝑖𝜌𝐴(𝑟)Ω2(cos𝜃)2𝐵2𝑟𝑖𝜌𝐴(𝑟)Ω2(cos𝜃)2𝐵𝑁𝑟𝑖𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵1𝑟𝑖𝜕𝑟𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵2𝑟𝑖𝜕𝑟𝜌Ω2𝜕𝜕𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝜕𝐵𝑁𝑟𝑖×𝜕𝑟𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵1𝑟𝑖𝜕𝑟2𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵2𝑟𝑖𝜕𝑟2𝜌Ω2𝐿𝑟𝑖𝜕𝐴(𝑟)𝔄𝑑𝑟2𝐵𝑁𝑟𝑖𝜕𝑟2×𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑖=1,2,,𝑁,𝑘=1,2,3,(15)where 𝔄 denotes (𝑟+𝑟0).

Based on the radial basis function technique, (11) take the following discrete forms: 𝐵1𝑟1𝐵2𝑟1𝐵𝑁𝑟1𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝜕𝐵1𝑟1𝜕𝑟𝜕𝐵2𝑟1𝜕𝑟𝜕𝐵𝑁𝑟1𝑎𝜕𝑟𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝐸𝐼𝑦𝑦𝑟𝑁𝜕2𝐵1𝑟𝑁𝜕𝑟2𝐸𝐼𝑦𝑦𝑟𝑁𝜕2𝐵2𝑟𝑁𝜕𝑟2𝐸𝐼𝑦𝑦𝑟𝑁𝜕2𝐵𝑁𝑟𝑁𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝜕𝐸𝐼𝑦𝑦𝑟𝑁𝜕𝜕𝑟2𝐵1𝑟𝑁𝜕𝑟2𝜕𝐸𝐼𝑦𝑦𝑟𝑁𝜕𝜕𝑟2𝐵2𝑟𝑁𝜕𝑟2𝜕𝐸𝐼𝑦𝑦𝑟𝑁𝜕𝜕𝑟2𝐵𝑁𝑟𝑁𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑦𝑦𝑟𝑁𝜕3𝐵1𝑟𝑁𝜕𝑟3𝐸𝐼𝑦𝑦𝑟𝑁𝜕3𝐵2𝑟𝑁𝜕𝑟3𝐸𝐼𝑦𝑦𝑟𝑁𝜕3𝐵𝑁𝑟𝑁𝜕𝑟3×𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝐵𝑘=1,2,3,1𝑟1𝐵2𝑟1𝐵𝑁𝑟1𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝜕𝐵1𝑟1𝜕𝑟𝜕𝐵2𝑟1𝜕𝑟𝜕𝐵𝑁𝑟1𝜕𝑟𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝐸𝐼𝑥𝑥𝑟𝑁𝜕2𝐵1𝑟𝑁𝜕𝑟2𝐸𝐼𝑥𝑥𝑟𝑁𝜕2𝐵2𝑟𝑁𝜕𝑟2𝐸𝐼𝑥𝑥𝑟𝑁𝜕2𝐵𝑁𝑟𝑁𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0],for𝑘=1,2,3,𝜕𝐸𝐼𝑥𝑥𝑟𝑁𝜕𝜕𝑟2𝐵1𝑟𝑁𝜕𝑟2𝜕𝐸𝐼𝑥𝑥𝑟𝑁𝜕𝜕𝑟2𝐵2𝑟𝑁𝜕𝑟2𝜕𝐸𝐼𝑥𝑥𝑟𝑁𝜕𝜕𝑟2𝐵𝑁𝑟𝑁𝜕𝑟2𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇+𝐸𝐼𝑥𝑥𝑟𝑁𝜕3𝐵𝑁𝑟𝑁𝜕𝑟3𝐸𝐼𝑥𝑥𝑟𝑁𝜕3𝐵2𝑟𝑁𝜕𝑟3𝐸𝐼𝑥𝑥𝑟𝑁𝜕3𝐵𝑁𝑟𝑁𝜕𝑟3×𝑎𝑘1𝑎𝑘2𝑎𝑘𝑁𝑇=[0]for𝑘=1,2,3.(16)By applying the radial basis function approach, (4) and (5) are substituted into (13). The following equations are then yielded𝜆2𝑟𝜌𝐴𝑖𝐵1𝑟𝑖𝜆2𝑟𝜌𝐴𝑖𝐵2𝑟𝑖𝜆2𝑟𝜌𝐴𝑖𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶0𝐵1𝑟𝑖𝜆𝐶0𝐵2𝑟𝑖𝜆𝐶0𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝜆𝐶1𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑦𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝐸𝐼𝑦𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝑟𝜌𝐴𝑖Ω2(sin𝜃)2𝐵1𝑟𝑖𝑟𝜌𝐴𝑖Ω2(sin𝜃)2𝐵2𝑟𝑖𝑟𝜌𝐴𝑖Ω2(sin𝜃)2𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵1𝑟𝑖𝑑𝑟𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵2𝑟𝑖𝑑𝑟𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵𝑁𝑟𝑖×𝑏𝑑𝑟𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵1𝑟𝑖𝑑𝑟2𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵2𝑟𝑖𝑑𝑟2𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵𝑁𝑟𝑖𝑑𝑟2×𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇=[0]for𝜆𝑖=1,2,,𝑁,𝑘=1,2,3,2𝑟𝜌𝐴𝑖𝐵1𝑟𝑖𝜆2𝑟𝜌𝐴𝑖𝐵2𝑟𝑖𝜆2𝑟𝜌𝐴𝑖𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶0𝐵1𝑟𝑖𝜆𝐶0𝐵2𝑟𝑖𝜆𝐶0𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝜆𝐶1𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝜆𝐶1𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝜆𝐶1𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑥𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑥𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵1𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵2𝑟𝑖𝑑𝑟2𝑑2𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑟2𝑑2𝐵𝑁𝑟𝑖𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+2𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵1𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵2𝑟𝑖𝑑𝑟32𝑑𝐸𝐼𝑥𝑦𝑟𝑖𝑑𝑑𝑟3𝐵𝑁𝑟𝑖𝑑𝑟3𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵1𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵2𝑟𝑖𝑑𝑟4𝐸𝐼𝑥𝑦𝑟𝑖𝑑4𝐵𝑁𝑟𝑖𝑑𝑟4𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝑟𝜌𝐴𝑖Ω2(cos𝜃)2𝐵1𝑟𝑖𝑟𝜌𝐴𝑖Ω2(cos𝜃)2𝐵2𝑟𝑖𝑟𝜌𝐴𝑖Ω2(cos𝜃)2𝐵𝑁𝑟𝑖𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵1𝑟𝑖𝑑𝑟𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵2𝑟𝑖𝑑𝑟𝜌Ω2𝑑𝑑𝑟𝐿𝑟𝑖𝐴(𝑟)𝔄𝑑𝑟𝑑𝐵𝑁𝑟𝑖×𝑑𝑟𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵1𝑟𝑖𝑑𝑟2𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵2𝑟𝑖𝑑𝑟2𝜌Ω2𝐿𝑟𝑖𝑑𝐴(𝑟)𝔄𝑑𝑟2𝐵𝑁𝑟𝑖𝑑𝑟2×𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇=[0]for𝑖=1,2,,𝑁,𝑘=1,2,3.(17)where 𝔄 denotes (𝑟+𝑟0).

According to the radial basis function approach, the boundary conditions in (14) have the following discrete forms: 𝐵1𝑟1𝐵2𝑟1𝐵𝑁𝑟1𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝑑𝐵1𝑟1𝑑𝑟𝑑𝐵2𝑟1𝑑𝑟𝑑𝐵𝑁𝑟1𝑏𝑑𝑟𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝐸𝐼𝑦𝑦𝑟𝑁𝑑2𝐵1𝑟𝑁𝑑𝑟2𝐸𝐼𝑦𝑦𝑟𝑁𝑑2𝐵2𝑟𝑁𝑑𝑟2𝐸𝐼𝑦𝑦𝑟𝑁𝑑2𝐵𝑁𝑟𝑁𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇=[0]for𝑘=1,2,3,𝑑𝐸𝐼𝑦𝑦𝑟𝑁𝑑𝑑𝑟2𝐵1𝑟𝑁𝑑𝑟2𝑑𝐸𝐼𝑦𝑦𝑟𝑁𝑑𝑑𝑟2𝐵2𝑟𝑁𝑑𝑟2𝑑𝐸𝐼𝑦𝑦𝑟𝑁𝑑𝑑𝑟2𝐵𝑁𝑟𝑁𝑑𝑟2𝑏𝑘1𝑏𝑘2𝑏𝑘𝑁𝑇+𝐸𝐼𝑦𝑦𝑟𝑁𝑑3𝐵1𝑟𝑁𝑑𝑟3𝐸𝐼𝑦𝑦𝑟𝑁𝑑3𝐵2𝑟𝑁𝑑𝑟3𝐸𝐼𝑦𝑦𝑟𝑁𝑑3𝐵𝑁