103學年度第一學期休假研究報告書：異向性彈性力學程式化之整合

(1)

國立成功大學

休假研究報告書

研究人員：胡潛濱教授系所：航空太空工程學系

(2)

2

目錄

page 目錄 2-3 休假研究報告(2014. 8. 1. ~ 2015. 1. 31) ……….. 4 附件一：德國錫根大學訪問研究期間所撰寫之期刊論文初稿

1. L. Xie, Chyanbin Hwu and C. Zhang, 2015, “Investigation of Green’s function and Its Derivatives for Three-Dimensional Anisotropic Elastic Solids,” to be submitted for publication.

2. Chyanbin Hwu, 2015, “Green’s Functions for Three-Dimensional Anisotropic Elastic Solids via Radon-Stroh Formalism.”

附件二： Content of “Anisotropic Elasticity with Matlab” 附件三：休假研究期間所撰寫之異向性彈性力學 Matlab 程式

1. Matlab Program of Boundary Element Analysis 附件四：德國錫根大學訪問研究報告

附件五：休假研究期間所發表及撰寫之期刊論文

1. Chyanbin Hwu and H.W. Chang, 2015, “Coupled Stretching-Bending Analysis of Laminated Plates with Corners via Boundary Elements,” Composite Structures, Vol. 120, pp.300-314. 2. Y.C. Shiah, C.L. Hsu and Chyanbin Hwu, 2015, “Direct

Volume-to-Surface Integral Transformation for 2D BEM Analysis of Anisotropic Thermoelasticity,” Computer Modeling in Engineering & Sciences, accepted for publication.

3. Chyanbin Hwu, Shao-Tzu Huang, and C. C. Li, 2015, “Stress Analysis of Multiple Holes and Cracks in Two-Dimensional Anisotropic Elastic Solids via a Special Boundary Finite Element.” (僅附大綱)

4. Chyanbin Hwu and H.W. Chang, 2015, “Singular integrals in boundary element analysis for unsymmetric laminated composites. (僅附大綱)

附件六：休假研究期間所撰寫將要發表之會議論文 (僅附摘要)

1. Chyanbin Hwu and H.W. Chang, 2015, “Singular integrals in boundary element analysis for unsymmetric laminated composites,” International Conference on Advances in Composite Materials and Structures (CACMS), Istanbul, Turkey. (accepted for presentation)

2. Chyanbin Hwu and Shao-Tzu Huang, 2015, “A Boundary Finite Element for Anisotropic/Piezoelectric Materials Containing Multiple Cracks,” International Conference on Mechanical Behavior of Materials (ICM12), Karlsruhe, Germany. (submitted for presentation, Nov. 26, 2014)

3. Chyanbin Hwu, 2015, “A unified near tip solution for multi-material anisotropic/piezoelectric wedges,” 9th International Conference on Advanced Computational Engineering and Experimenting (ACE-X), Munich, Germany. (to be submitted for presentation)

(3)

4. Chyanbin Hwu and Yu-Kuei Yeh, 2015, “Estimation of strength and fracture toughness for nanomaterials,” 20th International Conference on Composite Materials (ICCM20), Copenhagen, Denmark. (submitted for presentation, Nov. 26, 2014)

5. Chyanbin Hwu, C.L. Hsu and Y.C. Shiah, 2015, “Fundamental Solutions of Three-Dimensional Anisotropic Elastic Solids,” the 16th International Conference on Boundary Element Techniques, Valencia, Spain. (submitted for presentation, Dec. 25, 2014) 附件七：休假研究期間所進行及申請之研究計畫 (僅附摘要)

1. Chyanbin Hwu, 2015, "Boundary Element Design for Composite Laminated Plates (3/3)," National Science Council, R.O.C., NSC 101 - 2221 - E- 006 - 056 -MY3.

2. Chyanbin Hwu, 2015, "Estimation of Mechanical Properties for Nanomaterials (1/3)," Ministry of Science and Technology, R.O.C., MOST 103-2221-E-006 -161 -MY3.

3. Chyanbin Hwu, 2015, "Analytical and Numerical Analysis for Three-Dimensional Anisotropic Elastic Solids," Ministry of Science and Technology, R.O.C.

(4)

4

休假研究報告

(2014. 8. 1. ~ 2015. 1. 31) 胡潛濱, 2015. 2. 1 七年前我個人休假研究的主要目標之一為撰寫英文專業書籍，由於該段時間之全力付出，配合後續之接力，終於如願於2010 年由國際知名書商 Springer 出版了我個人

第一部英文專業書籍: Anisotropic Elastic Plates (673 pages)。該書出版至今紙本及電子

書之銷售及下載數量已超過 2000 本/次，以進階專業書籍而論，已如預期引起學界之重視。在此同時我也聽到許多讀者及好友反應，期望我能將我著作中所建立之理論及解析解程式化，就如同將高階之理論轉化為平易近人之電腦軟體。由於這類聲音也反應出常年來理論分析曲高和寡之心聲，因此我個人很希望藉由這次休假研究的機會創立以解析解為立基並方便使用者運算之力學分析軟體，以與現有之商用有限元素軟體區隔，提供學術界及業界具理論基礎之另類力學分析軟體。此一軟體除程式化我個人專業書籍中所提供超過 100 個問題之解析解外，還預計包括邊界元素之相關程式，基於此一需求，本次休假研究之初期(2014.9.1~11.30)配合國科會國外短期研究計畫至德國錫根大學(University of Siegen, Germany)訪問國際知名教授Prof. Chuanzeng Zhang (中文名: 張傳增)，期望藉由訪問研究之機會更進一步瞭解與邊界元素及有限元素相關之近期發展，以充實我個人預計發展之異向性彈性力學電腦軟體。由於張教授貼心之安排，訪問研究期間我住宿於該校學人宿舍(離研究室僅 5 分鐘之徒步距離)，並有自己單人之研究室及相關設施，同時亦不定期地(平均一星期至少二次)與張教授討論相關研究，經由密集的接觸與討論，這段其間我個人完成了 2 篇研究論文之初稿(如附件一)，同時也完成了近九成之書籍章節(如附件二)及異向性彈性力學運用在邊界元素分析之Matlab 程式(如附件三)。訪問期間之研究報告如附件四。訪問回國後，一如平日在校從事相關研究及指導學生，在這期間所發表及撰寫之期刊論文有 4 篇如附件五，所撰寫將要發表之會議論文有 5 篇如附件六，所進行及申請之研究計畫共3 件如附件七，所指導之學生共 13 位如附件八。為節省報告空間，部份附件僅列摘要、名稱及相關資訊。

(5)

Investigation of the Green’s functions for Three-Dimensional

Anisotropic Elastic Solids

Chyanbin Hwu, ………….. Institute of Aeronautics and Astronautics National Cheng Kung University, TAIWAN, R.O.C.

Abstract Keywords:

1. Introduction

2. Existing Green’s functions and their calculation methods

The governing differential equations, which the displacements uk, k=1,2,3, must

satisfy for a point force f applied at the origin x=0 of the infinite three-dimensional (3D) linear anisotropic elastic solids, can be written as

, + ( ) 0,

ijkl k lj i

C u  x f  (2.1) where ( )x is the Dirac delta function, repeated indices imply summation, and a comma stands for differentiation with respect to x . The solution to (2.1), which is _i the fundamental solution for the boundary element method, is usually called Green’s function.

Solutions in integral form

(i) Solution obtained from Fourier transform

By taking Fourier transforms on (2.1) and performing the inverse transform, it has been proved that the final solution of uk can be written in matrix form as (Ting and

Lee, 1997) 1 ( ) [ ] , 4 r  u x H x f (2.2a) where 1 * * * * 1 [ ] ( ) , ( ) . 2 s ds Qik C n nijkl j l  



 H x Q n n (2.2b) | |

r x and the integral is taken around the unit circle _|_n*_{| 1}_ _{on the normal plane to}

x. Note that [ ]H x is one of the three Barnett-Lothe tensors for the oblique plane whose normal is x, which depends not only on the elastic properties but also on the

(6)

2 direction of x.

(ii) Solution obtained from Radon transform

By taking a two-dimensional (2D) Radon transform on (2.1) without including the point force, it has been shown that the solution in 2D Radon space can be expressed by the Stroh formalism for 2D anisotropic elasticity (Wu, 1998). With the available 2D Green’s function (Ting, 1996; Hwu, 2010) together with the inverse Radon transform, the solution to (2.1) can be written as (Wu, 1998; Buroni and Denda, 2014; Hwu, et al., 2014)



1



3 2 ₀ sgn( ) ( ) Re ( ) ( ) ( ) , 2 T x z d          



  u x A A f (2.3a) where 1 2 3 ( ) cos sin . z_  x x   _x (2.3b) In (2.3a), sgn( ) 1x₃  for x₃  , 0 sgn( )x₃   for 1 x₃ ; Re stands for the real 0 part of a complex number; the angular bracket <> stands for a 3×3 diagonal matrix in which each component is varied according to its subscript  , e.g.,

1 2 3

diag[ , , ]

z z z z

  ; the superscript T stands for the transpose; _ and ( )A are, respectively, the Stroh’s eigenvalues and eigenvector matrix (see Appendix A) in the 2D Radon space spanned by the following two directions

(cos , sin , 0),   (0, 0, 1).

 

n m (2.4)

Calculation methods without integration

Till now most of the calculation methods proposed in the literature focus on the solution (2.2) obtained from the Fourier transform. As to the solution (2.3) obtained from the Radon transform, only the standard Gaussian quadrature rule is mentioned to perform the integral ranging from 0 to  . To calculate the integral involved in the matrix [ ]H x of (2.2), two main approaches have been proposed in the literature. (i) Calculate via Stroh’s eigenvectors (Ting, 1996; Nakamura and Tanuma, 1997)

[ ] 2_ _i T,

H x AA (2.5) where A is the Stroh’s eigenvector matrix associated with the oblique plane whose normal is the position vector x, and the oblique plane can be formed by any two mutually orthogonal unit vectors n and m. Since it has been proved that the result will be independent of the choice of n and m (Ting and Lee, 1997), we may choose

(7)

2 2 2 1 1 3 2 3 1 2 1 2 1 1 ( x x, ,0), ( x x , x x x, x ), if or x x 0, r          n m (2.6a) where 2 2 2 2 2 1 2, 1 2 3. x x r x x x       (2.6b) When x₁x₂ 0, 0x₃  , different but similar choice should be made such as

3 1

( x ,0, ) /x   

n , or n(0,x x₃, ) /₁ , and their associated m.

(ii) Calculate via Stroh’s eigenvalues (Ting and Lee, 1997)

4 ( ) 0 1 _ˆ [ ] , | | n n n q  



H x D T (2.7a) where 1 2 1 2 3 1 2 1 3 2 1 2 3 2 1 2 3 1 2 1 3 1 Re ... , for 0,1, 2, 2 ( )( ) 1 Re ... , for 2,3, 4. 2 ( )( ) n n n _n n n q n  _           _             _ _  _    _ _                  _ _ _      _ _ _ _      (2.7b)

In the above, the … represent two more terms obtained from the first term by a cyclic permutation of the subscripts; _k is the imaginary part of the Stroh’s eigenvalues

k

 , i.e., _k Im{ },_k k1, 2,3; _n₂  if n=2 and 1 _n₂  if 0 n2; T is defined in (A.1c); and _{D , n=0,1, 2,3,4, are the matrices related to the adjoint of D (see}ˆ( )n eqn.(A.4b)) by 4 ( ) 0 ˆ_{( )} nˆ n _, _{( ) ( ) | ( ) | ,}ˆ n       



 D D D D D I (2.8) where I is a 3 3 unit matrix.

Note that the oblique plane considered for the calculation of Stroh’s eigenvalues and eigenvectors can be formed by two unit vectors such as the ones shown in (2.6), which is different from the two vectors shown in (2.4) for 2D Radon transform.

The calculation methods stated in (2.5) and (2.7) only mentioned the Green’s function of displacements. To employ the Green’s function to the boundary element formulation, it is necessary to know the explicit expressions of its 1st and 2nd derivatives. To achieve this goal, till now most of the efforts emphasize on the expressions in terms of the Stroh’s eigenvalues such as Lee (…), Buroni… (…), and Xie and Zhang (2014). Since most of these expressions are complicated, for the

(8)

4

readers’ convenience only the most recent results presented by Xie and Zhang (2014) are listed in the Appendix B.

3. Calculation via Stroh’s eigenvectors

As shown in (2.2a) and (2.5), the Green’s function has been presented in terms of Stroh’s eigenvectors. However, no further discussions about the derivatives of Green’s function have been stated in the literature. Here, we like to go further to obtain the explicit expressions of the Stroh’s eigenvectors, and to obtain the derivatives of Green’s function based upon the Stroh’s eigenvectors. Differentiating (2.2a) with respect to x and further to _i x , we get _j





, 2 , 2 , 3 2 , , , 1 , 4 3 1 , 4 i i i i j ij ij i j j i ij x r r r x x x x r r r       _  _         __  _    _      u H H f u H H H H f (3.1a) where









, , , , , , , , , , 2 [ ] , 2 [ ] . T T T i i i T T T T T ij ij i j ij i j i i       H A A A A H A A A A A A A A (3.1b) From the results of (3.1), we see that that to calculate the derivatives of the Green’s function we need to know how to calculate the derivatives of the Stroh’s eigenvector matrix.

The Stroh’s eigenvector matrix for 2D anisotropic elasticity is a constant matrix related to the elastic compliances, and its explicit expression has been obtained through the Stroh-Lekhnitskii connection (see Appendix C). For 3D anisotropic elasticity, the Stroh’s eigenvector matrix depends not only on the elastic properties but also on the oblique plane whose normal is the position vector x. And hence, to get the explicit expression from the available solution of 2D anisotropic elasticity we need to consider a coordinate system based upon the position vector x and the oblique plane formed by any two mutually orthogonal unit vectors n and m. With this understanding, by following the concept of dual coordinate systems introduced in (Ting, 1996; Hwu, 2010) for the construction of the generalized eigen-relation of 2D anisotropic elasticity, the Stroh’s eigenvector matrix A of 3D anisotropic elasticity can be expressed as

*_,



(9)

where, if or x₁ x₂  0, 1 11 2 12 3 13 2 1 3 1 * 1 2 3 2 1 21 2 22 3 23 2 2 1 2 3 1 31 2 32 3 33 / / / / / / , . 0 / / c a c a c a x x x r x r x x x r x r c a c a c a x x r x r c a c a c a                _  _ _{ } _          Ω A (3.2b)

In (3.2b),  and r are defined in (2.6b), and a ,_ij i,j1,2,3 are the constants related to the reduced elastic compliances ˆS and the Stroh’s eigenvalues _ij _k, and

3 , 2 , 1 ,i

c_i are the normalization factors (see Appendix C). Note that both ˆS and _ij k

 are determined based upon the rotated coordinate system (n, m, x).

With the explicit expressions (3.2), the derivatives of the Stroh’s eigenvector matrix can be obtained as

* * , , , * * * * , , , , , , , , , i i i ij ij i j j i ij       A Ω A ΩA A Ω A Ω A Ω A ΩA (3.3a) where 1, 11 1 11, 2, 12 2 12, 3, 13 3 13, * , , 1, 21 1 21, 2, 22 2 22, 3, 23 3 23, 1, 31 1 31, 2, 32 2 32, 3, 33 3 33, + + + , + + + , + + + i i i i i i i i i i i i i i i i i i i i c a c a c a c a c a c a c a c a c a c a c a c a c a c a c a c a c a c a         _ _ _{ } _       _ _ Ω A (3.3b) and , 1 , ..., ...., k i k i c a   (待補) (3.3c) To calculate A and _,i A from (3.3), we need to know the derivatives of the Stroh’s _,ij eigenvalues and the reduced compliances, i.e., _{k i}_, , _{k ij}_, , Sˆ_{ij k}_, , and Sˆ_{ij kl}_, . By differentiating the eigen-relation (A.1) with respect to x and using the relations _i

T __ T

η N η where _η_( , )_{b a is the left eigenvector, we obtain}T

,i T ,i ,

  η N ξ (3.4a) where N is determined by its sub-matrices _,i N N_1,_i, _2,_i,N as _3,_i

1 1 1 1, ( ), , , 2, ( ) , , 3, , 1 1, ,, T T i    i   i i   i i   i  i i N T R T R N T N R N RN Q (3.4b) In (3.4b), 1 ,

(_T )_i_{can be obtained by differentiating}_TT1 __{I with respect to}

i x , which leads to 1 1 1 , , (T )_i  T T T _i  . (3.5)

(10)

6

Further differentiating (3.4a) with respect to x , we obtain _j

, , , , , ,,

T T T

ij i i ij i i

 η N ξ η N ξ η N ξ  (3.6) in which N can be determined by the way similar to (3.4b) and (3.5). _,ij

To get the derivatives of the reduced compliances Sˆ_{ij k}_, and Sˆ_{ij kl}_, , the transformation relation of 4th_{order tensor is differentiated as}

* *

, ( ), , , ( ), .

pqrs m pi qj rk sl m ijkl pqrs mn pi qj rk sl mn ijkl

S      S S      S (3.7) With the relation between the tensor notation and contracted notation (Hwu, 2010) and the relation (C.2) for the reduced elastic compliances, S and ˆ_{ij k}_, Sˆ_{ij kl}_, can be obtained.

With the results of (3.1), the Green’s function of stresses, tractions (on the surface with normal v) and their first derivatives can be obtained directly from the following relations

,, , , , ( , , , ).

ij C uijkl k l ti ij jv C u vijkl k l j ti j Cipkl u vk lj p u vk l p j

      (3.8)

4. Calculation via fundamental elasticity matrix in 2D Radon space

Although the solution obtained from 2D Radon transform has been obtained in (2.3), if no further reduction is made not only the numerical integration is necessary but also the calculation of Stroh’s eigenvalues and eigenvector matrices cannot be avoided. To have a further improvement on the numerical calculation, the following identity obtained for 2D anisotropic elasticity (Hwu, 2010)

1 2 2_r _z T sin ( ) _i(...),         A A N (4.1) is employed, and the results of (2.3) can be reduced to

3 2 2 ₀ sgn( ) sin ( ) ( , ) , 4 x d r       



u x N f (4.2) in which N₂( , )  is the Stroh’s generalized fundamental elasticity matrix in the 2D Radon space spanned by n and m set in (2.4).

By 2D Radon transform, Green’s function of the traction _{t on the surface with}n*

normal n* ( cos nsinm ) has also been obtained as (Hwu, et al., 2014)



2



3 2 ₀ sgn( ) ( ) Re ( ) (sin cos ) ( ) ( ) . 2 T x z d             



   * n t x B A f (4.3)

(11)

2 2 1 2_r (sin cos )_z T T( ) _i(...),             B A N (4.4a) where 2 2 1 1 1 1 1 1 1 2 3 ˆ

( ) cos [sin ( )] sin [sin ( )] cos 2 [cos sin ],

ˆ ( ) ( ) ( ).                      N N N N I N N N N N  (4.4b) Note that in (4.4b) the prime  denotes the derivative with respect to  , and the derivative of the generalized fundamental matrices N_i( ) can be performed by (Hwu, 2010) 2 1 1 2 3 2 1 2 2 1 3 3 1 1 3 ( ) { ( ) ( ) ( )}, ( ) { ( ) ( ) ( ) ( )}, ( ) { ( ) ( ) ( ) ( )}. T T                            N I N N N N N N N N N N N N N (4.5)

Substituting (4.4) into (4.3), we get

3 1 2 ₀ 2 sgn( ) 1 ( ) ( , ) . 4 T x d r      



* n t x N f (4.6) To find the derivatives of the Green’s function, firstly we should clarify the 2D coordinate used in the Radon transform, i.e.,

1cos 2sin cos , sin ,3

x x r x r         (4.7a) where 2 2 2 2 3 1 2 3 1 3 1 3 1 2 ( cos sin ) , tan tan . cos sin r x x x x x x x x                  (4.7b)

With the coordinate relations (4.7), the derivatives of ( )u x and _tn*( )_x _{can then be} obtained by differentiating (4.2) and (4.6) with respect to x , and the results are _i

* * 3 , 2 ₀ 3 2 2 * * 3 , 2 ₀ 4 1 1 sgn( ) 1 _ˆ _ˆ [ ( , ) ( , )] , 4 sgn( ) 1 [2 ( , ) ( , )] , 4 i i i T T i i i x x d r x x d r                  _    _  



* n u N N f t N N f (4.8a) where * * * 1 2 3 3 * * * 1 3 2 3 3 2 2 cos , sin , , cos , sin , , ˆ ( , ) sin ( , ). x x x x x x                         N N (4.8b)

(12)

8 5. Discussions

The calculation processes of the methods stated above are summarized in Table 1. From this Table, we see that…….

……….

(i) In this paper, both of the solutions obtained via the Stroh’s eigenvectors and the Radon transform are written in matrix form. According to the corresponding results of two-dimensional anisotropic elastic solids, they can all be employed to the three-dimensional piezoelectric solids by suitable expansion of the elastic constants to include the piezoelectric coupling effects (Hwu, 2010).

(ii) It should be noted that whether the solutions are obtained from Fourier transform or Radon transform, a two-dimensional plane is formed for the final mathematical formulation. However, these two planes are different (see Figure 1). The oblique plane used in the solutions obtained from Fourier transform is normal to the position vector x, and can be formed by any two mutually orthogonal unit vectors n and m such as the one shown in (2.6). Whereas the plane in the Radon transform is spanned by the unit vectors n and m given in (2.4). If the position vector x is written in terms of the spherical coordinate, i.e.,

( cos cos , cos sin , sin ),r   r   r  

x (5.1) The plane whose normal is x can be spanned by

( sin , cos , 0),   ( sin cos , sin sin , cos ),    

    

n m (5.2)

which is a special case of (2.6). In the Radon space, the normal to the plane is ( sin , cos , 0)   instead of the position vector x, since during the transform the variable x keeps unchanged. Thus, in the Fourier transform the plane is ₃

(13)

tangent to the point x located on the sphere, whereas in the Radon transform the plane is a slice vertical to the horizontal plane and normal to the tangent of the cylinder (see Figure 1).

(iii) Since the solutions obtained from Radon transform are expressed in terms of the elasticity fundamental matrix N, it will not have the degenerate problem because we donot need to calculate the Stroh’s eigenvalues as well as the Stroh’s eigenvectors. The only drawback when compared with the one via Stroh’s eigenvectors is the necessity of the numerical integration, which can be performed by the standard Gaussian quadrature rule.

(iv) Although the solutions obtained via Stroh’s eigenvectors may encounter problems when their associated eigenvalues have repeated roots, however in numerical calculation it can be avoided by adding a small perturbation on the elastic constants.

(v) In two-dimensional problem the Barnett-Lothe tensors are defined by

1 0 2 3 0 0 1 (2 ) ( ) , 1 1 2 ( ) , 2 ( ) . T T T i d i d i d                     



S AB I N H AA N L BB N (5.3)

When we employ these relations to the 2D Radon space spanned by (cos ,sin ,0) 



n and m(0,0,1), it may be written as

2 0 1 ...., ( ) 2 ( ) ( )_ _i _ _ T  ( , )_{  }_d , ...., etc.   



H A A N (5.4)

However, it should be careful for the angle  since here it denotes the direction of the 2D Radon space, and is not the angle between the dual coordinates of Stroh formalism. In other words, the integral appeared in the solution (4.2) cannot be written in terms of H( ) , i.e.,

2 0 1 ( , )d ( ).      



N H (5.5)

(vi) Comparing with the one obtained via Stroh’s eigenvalues, ….. advantage:…. Disadvantage:…..

(i) Stroh’s eigenvalues & eigenvectors:

Problems: repeated roots (例子: 強調重根處理除須留意外並無額外困擾)

Solution methods: (1) small perturbation (numerical, simple in calculation) (2) …… (analytical, ….)

(14)

10 (ii) Fundamental elasticity matrix:

Problems: line integral (例子: 強調積分點的選取大致需多少點, 因 integrand 簡易計算效率佳) Solution methods: (1) Gauss quadrature rule (numerical, simple in calculation)

(iii) Fourier transform & Radon transform: Different cutting plane in 3D domain 例子: 須告知若不同切面的處理情形, 例1: (2.6a)須改為: when 1 2 1 2 1 2 or 0, (2.6) or 0, ... or 0, .... x x x x x x      

例2: (4.2) & (4.6) & (4.8): when x₃0, ...?? 6. Comparison through numerical calculation

Comparison on Green’s function of displacements, tractions, and their first derivatives

7. Conclusions

Acknowledgements

The authors would like to thank National Science Council, TAIWAN, R.O.C. for support through Grant NSC ……, and ……….., GERMANY.

Appendix C: Stroh’s eigen-relation

The standard form of Stroh’s eigen-relation is for the 6 6 fundamental elasticity matrix N, and can be written as

,   Nξ ξ (A.1a) where 1 2 T 3 1 1 1 1 1 2 2 3 , , , , , T T T        _ _ _{  }           N N a N ξ b N N N T R N T N N RT R Q (A.1b) and , , .

ik ijkl j l ik ijkl j l ik ijkl j l

Q C n n R C n m T C m m (A.1c) In (A.1c), the two vectors n and m can be any two mutually orthogonal unit vectors, which form a plane for the Stroh’s eigen-relation. It has been proved that the result of eigenvalues and eigenvectors will be independent of the choice of vectors who form the plane (Ting and Lee, 1997). Also, the roots _k_,k1, 2,...,6 cannot be real if the

(15)

strain energy is positive (Eshelby, et al., 1953). Thus, there are three pairs of complex conjugates resulting from the vanishing of the determinant, |NI| 0 , and the roots are usually called the Stroh’s eigenvalue. Let Im_k 0, ,_k_₃ _k

3 , , 3 1, 2,3

k  k k  k k 

a a b b and assume that _k are distinct. The Stroh’s eigenvector matrices A and B are then the composition of the first three Stroh’s eigenvectors, i.e.,

1 2 3 1 2 3

[ , , ], [ , , ].

 

A a a a B b b b (A.2) Note that the eigen-relation (A.1) is constructed from the following relations

1

( T _ ) ( _ ) ,



    

b R T a Q R a (A.3) in which the second equality comes from

( )  , D a 0 (A.4a) where 2 ( )_ _{ }_( _ T)__ . D Q R R T (A.4b) In other words, if we have any difficulty to get solution from the standard eigen-relation (A.1), we may try the alternative approach to get  and a first from (A.4), and then b from (A.3).

Extended from the standard form (A.1), the generalized Stroh’s eigen-relation constructed in the dual coordinates can be expressed as (Hwu, 2010)

( )  ( ) , N ξ ξ (A.5a) where 1 2 T 3 1 ( ) ( ) _cos _sin ( ) , , ( ) , sin cos ( ) ( )   _ _ _              _ _ _{ }       N N a N ξ b N N (A.5b) and 1 1 1 2 2 1 3 ( ) ( ) ( ), ( ) ( ) ( ), ( ) ( ) ( ) ( ) ( ), T T T                     N T R N T N N R T R Q (A.5c) 2 2 2 2 2 2

( ) cos ( )sin cos sin ,

( ) cos ( )sin cos sin .

T T T                            Q Q R R T R R T Q R T T R R Q (A.5d)

Appendix B: Green’s function and its derivatives via Stroh’s eigenvalues ……….

(16)

12

Appendix C: Explicit expressions of 2D Stroh’s eigenvector matrix

1 11 2 12 3 13 1 1 2 2 3 3 3 1 21 2 22 3 23 1 2 3 3 1 31 2 32 3 33 1 1 2 2 3 , , c a c a c a c c c c a c a c a c c c c a c a c a c c c                   _ _ _ _           A B (C.1a) where 1 1 1 2 2 2 3 4 4 13 3 1 3 1 3 23 3 2 3 2 3 3 33 3 4 3 4 3 3 3 1 3 2 1 2 3 2 1 2 2 ( ) ( ) , [ ( ) ( )]/ , [ ( ) ( )]/ , 1, 2, ( ) ( ) , [ ( ) ( )]/ , [ ( ) ( )]/ , ( ) ( ) , , ( ) ( ) k k k k k k k k k k k k k k a p q a p q a p q k a p q a p q a p q l l l l                                                3 3 4 3 ( ) , ( ) l l    (C.1b) and 2 1 2 6 5 4 2 2 55 45 44 3 2 3 15 14 56 25 46 24 4 3 2 4 11 16 12 66 26 22 ˆ ˆ ˆ ˆ ˆ ( ) , ( ) , 1, 2, 4,5,6, ˆ ˆ ˆ ( ) 2 , ˆ ˆ ˆ ˆ ˆ ˆ ( ) ( ) ( ) , ˆ ˆ ˆ ˆ ˆ ˆ ( ) 2 (2 ) 2 , 1 j k k j j k j j k k j j k k k k k k k k k k k k p S S S q S S j l S S S l S S S S S S l S S S S S S k                                        , 2,3, (C.1c)

In the above, ˆS are the reduced elastic compliances which are related to the elastic _ij compliances S by _ij 3 3 33 ˆ i j _{ˆ .} ij ij ji S S S S S S    (C.2) k

 are the Stroh’s eigenvalues, and c_i,i1,2,3 are the normalization factors determined from the orthogonality relation, which requires that

2 2 3 2 1 3 23 3 3 13 3 33 1 1 , 1, 2, . 2( ) 2( ) k k k k k k c k c a  a  a a  a  a        (C.3) References

Buroni, F.C. and Denda, M., 2014, Radon-Stroh formalism for 3D theory of anisotropic elasticity, Beteq2014,….

Eshelby, J.D., Read, W.T. and Shockley, W., 1953, Anisotropic elasticity with applications to dislocation theory, Acta Metallurgica, 1: 251-259.

(17)

Hwu, C., 2010, Anisotropic Elastic Plates, Springer, New York.

Hwu, C., Xie, L.T. and Zhang, C.Z., 2014, Boundary element analysis for three-dimensional anisotropic elastic solids via Radon transform, under preparation.

Nakamura, G. and Tanuma, K., 1997, A formula for the fundamental solution of anisotropic elasticity, Q.Jl Mech. Appl. Math., 50 (2):179-194.

Ting, T.C.T., 1996, Anisotropic Elasticity: Theory and Applications, Oxford Science Publications, New York.

Ting, T.C.T. and V.G. Lee, 1997, The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids, Q.J. Mech. Appl. Math., 50 (3):407-426.

Wu, K.C., 1998, Generalization of the Stroh formalism to 3-dimensional anisotropic elasticity, Journal of Elasticity, 51: 213-225.

……….

List of Figures List of Tables

(18)

14

Table 1. Calculation processes for three different methods Given elastic constants C , and position _ijkl ( , , )x x x ₁ ₂ ₃

Via Stroh’s eigenvalues Via Stroh’s eigenvectors Via Radon transform (2.6a): n, m (2.6a): n, m (2.4): n, m (set  for integration) (A.1c): Q, R, T (A.1c): Q, R, T (A.1c): Q, R, T

(A.1b): N1, N2, N3, N (A.1b): N1, N2, N3, N (A.1b): N1, N2, N3, N

(A.1a):  (A.1a), (A.2): , A, B (4.7b): r,

(2.7b): q_n _{(2.5): H[x]} (A.5d): ( ), ( ), ( )Q R T (A.4b),(2.8): _Dˆ( )n _{(2.2a): u(x)} (A.5c):

1( ), 2( ), 3( )

N N N

(2.7a): H[x] (4.2): u(x) by Gauss integration

(2.2a): u(x)

Continue for the derivatives

……. (3.4b): N N_1,_i, _2,_i,N_3,_i _(4.8b): _N_{ˆ ( , )}₂ _{ } (4.5)2: Nˆ ( , )2   (3.4a): _,i (4.4b): N₁( , )  (4.5): N₁( , )  (3.7): S ˆ_{ij k}_, (4.6): t(x) by Gauss integration (3.3): A _,i (4.8a,b): u x t x_,_i( ), ( )_,_i by Gauss integration (3.6): _,ij (3.7): Sˆ_{ij kl}_, (3.3): A _,ij (3.1b): H H _,_i, _,_ij (3.1a): u x u x_,_i( ), _,_ij( ) (3.8): _ij( ), ( ), ( )x t x t x_,_i

Note: In Radon transform, for simplicity sometimes the dependence on variable  is not written explicitly. For example, N₁( ) means N₁( , )  , and N₁( ) means

1( , ) 

(19)

Boundary Element Analysis for Three-Dimensional Anisotropic

Elastic Solids via Radon Transform

Chyanbin Hwu, ……

1_{Institute of Aeronautics and Astronautics}

National Cheng Kung University, TAIWAN, R.O.C. E-mail: [email protected]

Abstract Keywords:

1. Introduction ……….

2. Radon-Stroh formalism for 3D anisotropic elasticity

In a fixed rectangular coordinate system x_i,i1,2,3, let ui, , and ij ij be, respectively, the displacement, stress and strain. The constitutive laws for the linear anisotropic materials, the strain-displacement relations for the small deformations, and the equilibrium equations for static loading conditions (if the body forces are neglected) can be written as (Hwu, 2010)

, , , 1 , ( ), 0, , , , 1, 2,3, 2 ij Cijkl kl ij ui j uj i ij j i j k l          (2.1)

where the repeated indices imply summation; a comma stands for differentiation and the elastic constants C_ijkl are assumed to be fully symmetric and positive definite. Substituting (2.1)2 into (2.1)1 with C_ijkl fully symmetric, we get

. ,l k ijkl ij C u  (2.2) The governing differential equations which the displacements must satisfy can then be obtained by substituting (2.2) into (2.1)3, i.e.,

. 0 ,lj  k ijklu C (2.3) Consider a Radon transform from the space of rectangular Cartesian coordinate

1 2 3

( , , )x x x 

x to the space of cylindrical coordinate ( , , )  x₃ by

3 ( , , ) ( ) ( ) , f  x  f   d  

_

  n x n x x  (2.4) where  is the Dirac delta function, and n( , ,0)n n₁ ₂ with n₁ cos , sin n₂  

(20)

2

is a unit vector normal to the line n x . Note that in this transformation the coordinate variable x keeps unchanged, and the argument written as ₃ ( , , )n  x₃ has the same meaning as ( , , )  x₃ since the unit vector n is a function of  .

To apply the Radon transform (2.4), the governing equations (2.3) are rewritten by separating the indices into ( , )x x and ₁ ₂ x as ₃

, ( 3 3 ) ,3 3 3 ,33 0, , =1,2,3, , 1, 2.

i k k i k i k k i k k

C_{ }u _  C_ C _ u _ C u  i k    (2.5) Knowing that R f{ }_,_ n f__,_ and R f{ }_,3  f_,3, application of the Radon transform (2.4) to (2.5) leads to

, ( 3 3 ) ,3 3 3 ,33 0.

i k k i k i k k i k k

C_{   }n n u _  C_ C _ n u_ _ C u  (2.6) Since the governing differential equations (2.6) in the domain of Radon transform are a set of homogeneous second-order differential equations which depend only on two variables  and x , a general solution for ₃ u_k can be written as

( ),

k k

u a f z or in matrix form u a f z( ), (2.7a) where

3.

z  x (2.7b) In (2.7a), a and  are constants to be determined and f is an arbitrary function of z. Substituting (2.7a)1 into (2.6), we get

2

{Q_ik (R_ikR_ki) T a f z_ik} _k( ) 0, (2.8a) or in matrix form

2

{_Q__(_{R R}_ T)__ _{T a 0 (2.8b)}} _ , where the superscript T stands for the transpose, and

, , ,

(cos , sin , 0), (0, 0, 1).

ik ijkl j l ik ijkl j l ik ijkl j l

Q C n n R C n m T C m m

 

  

 

n m (2.8c)

A nontrivial solution of a exists if

2

( T) 0,

 

   

Q R R T (2.9) which gives a sextic equation for . Note that  and a determined from (2.8), which are generally called material eigenvalue and material eigenvector (sometimes are called Stroh’s eigenvalue and Stroh’s eigenvector), now depend not only on the elastic constants C but also on the variable _ijkl  .

(21)

After obtaining  and a for the displacements in the transform domain, the _k stresses in the transform domain can be determined from (2.2) by employing (2.4) and (2.7). The result is

3

( ) ( ),

ij C nijk  Cijk a f zk

    (2.10) which then leads the tractions on the surfaces with normals n and m to

( ) ( ), ( ) ( ).

i ij j ik ik k i ij j ki ik k

tn  n  Q R a f z tm  m  R T a f z _(2.11) Introducing a new vector b as

1

( T _ ) ( _ ) ,



    

b R T a Q R a (2.12) in which the second equality of (2.12) is obtained from (2.8b). With (2.12), equation (2.11) can now be written as

( ), ( ).

i i i i

tn  b f z tm b f z _(2.13)

Introduce the stress functions _i, i1,2,3, such that

,3, , . i i i i tn   tm _ _(2.14) we get ( ) i b f zi    , or in matrix form  bf z( ). (2.15) In two-dimensional anisotropic elasticity, it has been proved that the six roots of material eigenvalues  cannot be real and are three pairs of complex conjugates (Ting, 1996). Thus, by letting Im_k 0, ,_k_₃ _k a_k_₃ a_k, b_k_₃b_k,k=1,2,3, and assuming that _k are distinct, the general solution obtained by superposing six solutions of (2.7a)2 and (2.15) can be written as









2 Re ( ) ,z 2 Re ( ) ,z   u Af 　 Bf (2.16a) where

 

1 1 21 32 2 3 3 1 2 3 [ ], [ ], [ ( ), ( ), ( )] ,T z f z f z f z    A a a a B b b b f     (2.16b)

and Re stands for the real part of a complex number. As stated in (Buroni and Denda, 2014), the key advantage of Radon-Stroh formalism when compared with the other methods is that the transformed space corresponds to physically meaningful 2D problem. With this understanding, to state the traction prescribed boundary conditions in the transformed domain, we need to know the relation between the stress function

(22)

4 vector  and the traction  _n*

t 

on the surface with normal n*_{(see Figure 1). The} normal n* of the surface edge of the transformed plane can be expressed as

* __cos_ __sin_ _,

n n m (2.17) where n and m are given in (2.8c). With (2.17) and the relations (2.11) and (2.14), by Cauchy’s formula we obtain

* i_{, or in matrix form} _, i ij j t n s s          * n _t       _(2.18)

where s is a parameter denoting the tangential direction which has to be chosen such that when one faces the direction of increasing s the solid body lies on the right side for the coordinate system shown in Figure 1. Otherwise, a negative sign should be put on the right hand side of equation (2.18). Note that during the derivation of (2.18), the chain rule as well as the relations cos  x₃/ and sins     have been / s used.

Figure 1.

In the above, all the formulae have been purposely organized into the mathematical form of Stroh formalism for 2D anisotropic elasticity. Thus, some important relations and identities developed in Stroh formalism (Ting, 1996; Hwu, 2010) can all be employed in the domain of Radon transform for 3D elastostatic analysis. For example, the well-known material eigenrelation Nξ ξ, and the Barnett-Lothe tensors S, H and L, can be stated as follows.

1 2 T 3 1 1 1 1 1 2 2 3 , , , , , . T T T          _ _ _{  }           Nξ ξ N N a N ξ b N N N T R N T N N RT R Q (2.19)

(23)

1 0 2 3 0 0 1 (2 ) ( ) , 1 1 2 ( ) , 2 ( ) . T T T i d i d i d                     



S AB I N H AA N L BB N (2.20)

As stated in the note below equation (2.9), for simplicity the dependence on variable  for some quantities of the Radon-Stroh formalism is not written explicitly unless stated otherwise. For example,

( ), ( ), ( ), ( ), ( ), ( ), ( ), ( ), ( ), ( ),...,                       Q Q R R T T A A B B N N S S H H L L (2.21)

and hence the explicit form of (2.20) is

1 0 1 ( )_ _i{2 ( ) ( )_ _ T }  ( , )_{  }_d , ...., etc.    



S A B I N (2.22)

In view of the expression of (2.22), one should not be confused about the roles of angular parameters  and  . The former is a variable in the transform domain, whereas the latter is the angle between the dual coordinates ( , ) x₃ and * *

3

( , ) x considered in the Stroh formalism for 2D problems.

Inverse of Radon transform

After finding the solutions in the transform domain, the 3D solutions for the original problem should be recovered by applying the inverse Radon transform. It’s known that (….) 0 1 ( ) ( , ) , 2 f  f  d  



 x  n x (2.23a) where , ( , ) 1 ( , ) f rr . f dr r          



  (2.23b) The symbol  is a notation for Hilbert transform, whose integral can be derived analytically by considering the contour integral of f z( , ) /( z) around a very large semi-circle involving a small semi-circular detour around the pole  in the upper (or lower) half-plane. If the Green’s function is considered for three-dimensional elastostatic analysis whose unit point force is applied on

1 2 3

ˆ ˆ ˆ

(24)

6

hence the Hilbert transform (2.23b) can be obtained as

3 ,

( , ) sgn( ) ( ),

f    i x f_  (2.24) where sgn(x₃) 1 for   , x₃ 0 sgn(x₃) 0 for   , x₃ 0 sgn(x₃)  for 1

3 0

x

  , and  x₃ x₃ . xˆ₃

3. 2D Boundary element in Radon transform domain

For 3D analysis of anisotropic elastic solids, several different boundary element methods have been proposed in the literature such as (….). Due to the complexity of the Green’s function, most of them consider how to expedite the calculation of Green’s function and its derivatives (…..). By using the Radon-Stroh formalism stated in the previous section, it is not difficult to get these solutions from their corresponding 2D problems. For example, we can get the Green’s function for 3D anisotropic solids from the corresponding 2D solution (2.16) in Radon domain by applying the inversion stated in (2.23) and (2.24). It’s known that the Green’s function for 2D problem can be expressed as (2.16) with (Ting, 1996; Hwu, 2010)

ˆ ( )z  ( ) ,z f F p (3.1a) where 3 1 _ˆ _ˆ _ˆ _ˆ ( ) ln( ) , . 2 T z z z z x i             F A (3.1b)

In (3.1), i 1, and pˆ ( ,p p pˆ ˆ ˆ₁ ₂, ₃) is the line force applied on the point ( , )ˆ ˆx₃ ; the angular bracket <> stands for a 3×3 diagonal matrix in which each component is varied according to its subscript  , e.g., z diag[ , , ]z z z₁ ₂ ₃ . Substituting (2.16) with (3.1) into (2.23a) with (2.24), we get









0 0 1 1 ( ) Re ( ) ( , )  d, ( ) Re ( ) ( , )  d,   

_



_

u x A f x 　 x B f x (3.2a) where 3 3 1 2 3 sgn( ) 1 _ˆ ( , ) ( ) , ˆ 2 ( ) ( ) cos sin . T x z z z x x x x                           f x A p n x  (3.2b)

Use of the identities obtained for 2D problems, i.e., (Hwu, 2010)

1 2 1 1 2 sin ( ) (...), 2 cos sin ( ) (...), T T T r z i r z i                    A A N B A I N (3.3)

(25)

the results of (3.2) can be further reduced to 3 2 2 ₀ 3 1 2 ₀ sgn( ) sin _ˆ ( ) ( , ) , 4 sgn( ) 1 ˆ ( ) [cos sin ( , )] . 4 T x d r x d r    _{ } _            



u x N p x I N p  (3.4)

With (3.2a)2, the traction

* n

t on the surface with normal n* can then be obtained from (2.18) via the inverse of Radon transform, and the result is

3 2 ₀ 2 sgn( ) sin cos ˆ ( ) Re ( ) ( ) , ˆ 2 ( ( ) ) T x d z z  _                _   _   



* n t x B A p 　 (3.5)

in which the angle  is related to the direction n* by (2.17). Use of the identities for 2D problems (Hwu, 2010), equation (3.5) can also be reduced (Hwu, et al., 2014).

Using the Green’s functions obtained in (3.2)-(3.6), by following the standard procedure for boundary element formulation (Brebbia, et al., 1984) a 3D boundary element can be established for elastostatic analysis of linear anisotropic elastic solids. In this paper, to avoid the complexity of 3D formulation and implementation, a 2D boundary element in Radon transform domain is now established as follows.

Since in the Radon transform domain, only two variables,  and x , are ₃ involved, the boundary integral equations can be set as those of 2D problems, i.e., (Brebbia, et al., 1984) * * ( ) ( ) ( , ) ( ) ( ) ( , ) ( ) ( ), , 1, 2,3, ij j ij j ij j c u t u d u t d i j   

_

 

_

  ξ  ξ  ξ χ  χ χ  ξ χ  χ χ (3.6)

where χ( , ) x₃ and ξ( , )ˆ ˆx₃ ,  denotes the boundary of the elastic solid in the Radon domain; u_j( )χ and ( )t_j χ are the Radon transform of the displacements and surface tractions along the boundaries of the Radon domain; c_ij(ξ) is a coefficient dependent on the location of ξ , which equals to _ij/2 for a smooth boundary and c_ij _ij for an internal point. The symbol _ij is the Kronecker delta, i.e., 1_ij  when i and j _ij 0 when i . In practical applications, j c_ij(ξ) can be computed by considering rigid body motion. *_{( , )}

ij u ξ χ and *_{( , )} ij t ξ χ  are the fundamental solutions of displacements and tractions for 2D elastic problems, which can be obtained from (2.16) with (3.1), and written in matrix form as

* *

,

[ ] 2 Re{[ ( )] }, [ ] 2 Re{[T ( )] },T

ij ij s

u  AF z t  BF z (3.7) where the symbol F_,_s   F/ s denotes the differentiation with respect to the

(26)

8

tangential direction s of the body in Radon domain, and hence

, sin cos 1 ( ) . ˆ 2 T s z _i _z _z            F A (3.8) After getting the fundamental solutions in (3.7), the unknowns remained in the boundary integral equations (3.6) are u_j and t_j over the boundary  . In boundary element formulation, the boundary  is approximated by a series of elements, and the points χ , displacements u_j and tractions t_j on the boundary are approximated by the nodal points, nodal displacement and nodal traction through suitable interpolation functions. In this study, the linear variation within each element is assumed for the boundary points, and the displacements and tractions in Radon transform domain. Thus, by following the standard procedure for boundary element formulation, the Radon transform of the displacements and tractions of the boundary nodes can be obtained by solving a system of linear algebraic equations. Once all the Radon transform of the displacements and tractions on the boundary are determined, the Radon transform of the displacement at the internal point can be calculated by using the boundary integral equations (3.6) again, where c_ij( )ξ _ij, i.e.,

* * ( ) ( , ) ( ) ( ) ( , ) ( ) ( ), , 1, 2,3. j ij j ij j u u t d t u d i j   

_

 

_

  ξ ξ χ  χ χ  ξ χ χ χ    _(3.9)

The internal strains at point ξ can be found by differentiating (3.9) with respect to i

 and using the strain-displacement relation _ij (u_i_,_ju_j_,_i)/2. To find the internal stresses at point ξ , the stress-strain law _ij C_ijks_ks should be employed. Therefore, for the calculation of the internal strains and stresses, we need to know the derivatives of *

ij

u and *

ij

t . From the fundamental solutions (3.7), we get

* * ˆ ˆ ˆ ˆ , , , , [ ] 2 Re{[ ( )] }, [T ] 2 Re{[ ( )] },T ij k k ij k sk u  AF z t  BF z (3.10) where F_,_kˆ   F/ k, k 1, 2,  

(₁ ˆ, ₂ xˆ₃ ), denotes the differentiation with respect to the source point ξ of the Radon domain. From (3.1b) and (3.8), we obtain

ˆ ˆ ,1 ,2 ˆ 2 , 1 ˆ 2 , 2 1 1 1 ( ) , ( ) , ˆ ˆ 2 2 sin cos 1 ( ) , ˆ 2 ( ) (sin cos ) 1 ( ) . ˆ 2 ( ) T T T s T s z z i z z i z z z i z z z i z z                                             F A F A F A F A     (3.11)

(27)

After getting the displacements u_i and stresses _ij in the Radon transform domain through the 2D boundary element, their associated solutions in the 3D physical domain of can then be determined by numerical inversion of Radon transform.

4. Numerical inversion of the Radon transform ……… 5. Numerical Examples ……… 6. Conclusions ……… Acknowledgements

The authors would like to thank National Science Council, TAIWAN, R.O.C. for support through Grant NSC ……, and ………..

References

Brebbia, C.A., Telles, J.C.F. and Wrobel, L.C., 1984, Boundary Element Techniques: Theory and Applications in Engineering, Springer-Verlag, Berlin.

Buroni, F.C. and Denda, M., 2014, Radon-Stroh formalism for 3D theory of anisotropic elasticity, Beteq2014,….

Hwu, C., 2010, Anisotropic Elastic Plates, Springer, New York.

Kogl, M. and Gaul, L., 2000, A 3-D boundary element method for dynamic analysis of anisotropic elastic solids, CMES, 1 (4): 27-43.

Ting, T.C.T., 1996, Anisotropic Elasticity: Theory and Applications, Oxford Science Publications, New York.

Ting, T.C.T. and V.G. Lee, 1997, The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids, Q.J. Mech. Appl. Math., 50 (3):407-426.

Wu, K.C., 1998, Generalization of the Stroh formalism to 3-dimensional anisotropic elasticity, Journal of Elasticity, 51: 213-225.

………. List of Figures List of Tables

(28)

(29)

103學年度第一學期休假研究報告書：異向性彈性力學程式化之整合

國 立 成 功 大 學

休假研究報告書

目錄

休假研究報告

Investigation of the Green’s functions for Three-Dimensional

Anisotropic Elastic Solids

















































Boundary Element Analysis for Three-Dimensional Anisotropic

Elastic Solids via Radon Transform











 







































Anisotropic Elasticity with Matlab

Contents

國立成功大學

_

_

_

_

_

_

_