关于3种计算张量函数方法的比较

第33卷第4期北　京　交　通　大　学　学　报Vol. 33No. 4

　　　　　　　　　　　　　

文章编号:167320291(2009) 0420106204

关于3种计算张量函数方法的比较

兑关锁1, 王志乔1, 马　莎2

(1. 北京交通大学土木建筑工程学院, 北京100044;2. 华北油田一中, 河北任丘, 062552)

摘　要:对计算各向同性对称张量函数的3种方法(谱分解法、) 进行了分析比较, 函数进行了实例计算. 结果表明, , 何选择计算张量函数的方法提供参考.

关键词:连续介质力学; 谱分解; ; 中图分类号:O331:A

on Three Methods for the Computation

of the Isotropic T ensor Function

DU I Guansuo 1, W A N G Zhiqiao 1, M A S ha 2

(1. School of Civil Engineering , Beijing Jiaotong University , Beijing 100044,China ; 2. First Middle School of Huabei Oilfield Company , Renqiu Hebei 062552,China )

Abstract :Three methods for the computation of the isotropic tensor function , spectral decomposition method , eigenprojection operator method and representation theorem method , are compared. Different approaches to calculate the eigenvalues and their corresponding eigenvectors are also discussed. As an example , the exponent tensor function is calculated. It is shown that the representation theorem method takes shortest time , followed by the spectral decomposition method. The results could provide guidance for choosing different methods in writing programs for calculation.

K ey w ords :continuum mechanics ; spectral decomposition ;eigenprojection operator ;representation the 2orem

　　各向同性张量函数的计算在连续介质力学和计算力学等领域中都有着广泛的应用. 例如在有限变形理论中的对数应变为伸长张量的函数, 应力为应变张量的函数等. 特别在非线性计算中, 张量函数的计算速度将直接影响到各种最终结果的计算效率. 通常计算张量函数的方法可归纳为3种:

1) 利用Hill [1]提出的主轴法给出. 该方法一般能将复杂的张量函数简洁地表出, 但由于所给出的表示是在主轴坐标下, 要获得在任意坐标下的表示, 首先要计算特征向量再进行坐标变换. 对一些复杂

问题, 如张量函数的导数等, 其计算工作量较大[2].

2) 根据特征投影算子的形式获得[3]. 近年来, 一些文献利用此表示方法对对数、指数及级数形式函数等的导数进行了数值分析和研究[4-6]. 虽然这种方法不需要直接求解特征向量, 但是需要计算3个特征投影算子.

3) 基于各向同性张量函数表示定理[7]. 此方法在计算张量函数时, 只需确定3个系数. 对于求张量函数导数等复杂问题, 利用这3个系数也可获得更简洁的表示公式[8-10].

　　收稿日期:2008-11-06

基金项目:国家自然科学基金资助项目(10772021)

) , 男, 河南清丰人, 博士, 博士生导师. em ail :gsdui @bjtu.edu. cn. 作者简介:兑关锁(1963—

第4期　　　　　　　　　　　　　兑关锁等:关于3种计算张量函数方法的比较

107

　　本文目的是对众所周知的张量函数3种计算方法进行定量性分析比较, 从而为在复杂问题及数值计算中, 如何选用张量函数计算方法提供依据. 由于实际问题多数为二维和三维问题, 而二维问题过于简单, 故后面的讨论限于三维空间. 首先从计算二阶对称张量特征根与特征向量出发, 给出了特征根的两种计算公式和计算特征向量的显式表示公式. 最后通过利用各向同性对称张量函数的3种计算方法(谱分解法、特征投影算子法和表示定理法) 对指数张量函数进行了实例计算.

式(7) 中的±与(2I 31-9I 1I 2+27I 3) 正负相同.

若令C 的分量为C ij , N i =

(N i 1, N i 2, N i 3) T 和

A k =C -Λk I , 则特征向量可通过下面方程得到

A 11A 21A 31

k k k

A 12A 22A 32

k k

k

A 13A 23A k k

k

N k 1N k 2=N k 00(8)

(k =1, 2, 3)

根据定义知, A k 的3. 再由式(8) 可知, . 因(8) , 而通.

, 且令k k k k k b 1=A 12A 23-A 13A 22,

b 2=A 13A 21-A 11A 23, b 3=A 11A 22-A 12A 21,

k

k

k

k

k

k

k

k

k

k

1　二阶对称张量特征值与特征向量

设C 为二阶对称张量, i , 向为N i . I 1=tr C 22

I 2=[(tr C ]=Λ1Λ2+Λ2Λ3+Λ1Λ3

2

I 3=det C =Λ1Λ2Λ3

(

1)

|b k |=

k 2k 2k 2(b 1) +(b 2) +(b 3) 1

则3个特征向量为

N k =|b |

k

-1

k k

(b 1　b 2　b k 3)

T

可求得特征根为

Λ1=I 1+2(I 21-3I 2)

3

(9)

2

×

(2)

当C 有两个相同特征根时, N 1可以通过相同的方法获得.

3arccos 　32(I 2-3I ) 2

1

2

2　计算各向同性张量函数的3种方法

若张量函数T (C ) 对任意正交张量Q 满足

T (QCQ T ) =Q T (C ) Q T

(T , C ) 共轴, 即TC =CT .

(10)

利用特征方程

Λ3i -I 1Λ2i +I 2Λi -I 3=0可得恒等式

Λ2j +Λj Λi +Λ2i -I 1(Λj +Λi ) +I 2=0

(i ≠j )

(3)

则T (C ) 称之为各向同性张量函数. 由此知共轭对

当C 给定时, 通常可以通过3种方法计算各向同性对称张量函数T (C ) . 211　谱分解法

(4)

当Λ1、I 1和I 2给定, 则另外两个特征根Λ2和Λ3可由恒等式(4) 解出[11]

Λ2

2(I 1-Λ1±I 1Λ1-3Λ2=1-4I 2+I 1) 2Λ(5)

当C 具有3个不同特征根时, 由谱分解定理知

C =

ΛN ∑

i

i

i

N i N i

(11)

则各向同性张量函数T (C ) 可以由主轴N i 表示为

T =

t N ∑

i

i

或直接解特征方程(3) 得Λ22

() I +2I -3I =1122×3Λ3arccos 332(I 2-3I ) 2

1

2

i

(12)

式中标量函数t i (Λ1, Λ2, Λ3) 具有对称性[1], 即

t 1(Λi , Λj , Λk ) =t 1(Λi , Λk , Λj ) =t i (Λ1, Λ2, Λ3)

(6)

(13)

式中(i , j , k ) 为(1, 2, 3) 的偶排列.

Λ2=Λ3=当C 有两个相同特征根时, 设Λ1≠Λ0, 则有

C =Λ2I +(Λ1-Λ2) N 1 N 1

(14) (15)

当C 有相同特征根时, 设Λ2=Λ3=Λ0, 则特征根

的计算公式退化为

Λ1=I 1±2(I 21-3I 2) 2,

3

Λ2=Λ3=(I 1-Λ1) (7)

2

和

T =t 2I +(t 1-t 2) N 1 N 1

式中I 为二阶单位张量.

108

北　京　交　通　大　学　学　报　　　　　　　　　　　　　　　第33卷

当C 有三重特征根时, 设Λ1=Λ2=Λ3=Λ和

t 1=t 2=t 3=t , 则

C =ΛI T =t I

(16) (17)

φ 0=

ΛΛΛ2-Λ1,

(29)

φ 1=ΛΛ

2-1

212　特征投影算子法

若定义特征投影算子P i 为

P i =N i N i

(18)

3　计算实例与分析

为了比较不同计算方法的计算效率, 具体对给定的张量指数函数进行计算. 设T (C ) =e C , 当C 给定, 且其分量为

C ij =

由特征向量N i 的正交性, 得特征投影算子的特性为

P i N j =δij N j , P i P j =δij P j ,

-5, 6

∑

i

P i =I ,

5

) 可求得不变量为

I 1=0, 　I 2=-84, 　I 3=295.

(P i N i P i =

() ()

(Λi -Λj ) (Λi -Λk ) T =

t P ∑

i

i

当C 具有其他量可按不同的方法求得(本文所有结果均使用

Celeron (R ) CPU 210GHz , 512MB 内存计算机计算获得. ) . 311　特征根计算

(20) (21)

因此

i

特征根可通过3种不同的方法获得:①雅可比法; ②由式(2) 和式(6) 给出; ③由式(2) 和式(5) 给出. 所得结果为

Λ1=-51578, 　Λ2=-5, 　Λ3=101578. 　　将运算循环100万次, 则所花费的时间分别为:

当C 有两个相同特征根时, 特征投影算子为[3]

(Λ) P 1=,

(Λ1-Λ2)

P 2

(Λ) =

(Λ2-Λ1)

(22) (23)

方法①,414s ; 方法②,316s ; 方法③,217s. 计算结果表明, 采用方法②, 即3次方程求根公式, 所用时间仅为数值求解的雅可比法的80%.而利用方法③, 即不变量求根公式(5) , 特征根计算时间仅为方法①的60%.312　特征向量计算

和

T =t 2I +(t 1-t 2) P 1

213　表示定理法

对于对称各向同性张量函数T , 由表示定理可知[7]

2

T =φ0I +φ1C +φ2C

(24) (25)

式中φ0, φ1和φ2为C 的不变量函数, 由式(24) 有

2t i =φ0+φ1Λi +φ2Λi

特征根可通过2种不同的方法获得:①雅可比

法; ②由式(9) 给出. 所得结果为

T

N 1=(01781, -01592, -01201) ,

当C 有3个不同特征根时, 方程组(25) 有唯一解

φ2=

f ∑

i

N 2=(-1/3, -2/3, 2/3) T ,

T

N 3=(01529, 01453, 01718) .

i

φ1=I 1φ2+φ0=I 3

f Λ∑

i

i

i

(26)

利用数值求解的雅可比法及计算特征根的第③种方法和计算特征向量的第②种方法进行计算, 若循环100万次, 所花的时间分别为:方法①,517s ; 方法

∑

i

f i /Λi

式中:f i =t i (Λk -Λj ) /Δ, 其中

Δ=(Λ1-Λ2) (Λ2-λ3) (Λ3-Λ1) .

当C 有两个相同特征根时, 则T 的表示退化为

T =φ 0I +φ 1C

(27) (28)

②,313s. 计算结果表明, 理论公式计算用时不到数值计算用时的60%.313　张量函数计算

张量函数可通过3种不同的方法获得:①谱分

φ 0+φ 1Λi =t i 　　(i =1, 2)

解法; ②特征投影算子法; ③表示定理法.

由式(26) 有

φ0=413479×103,

因此

求解方程组(28) 有

第4期　　　　　　　　　　　　　兑关锁等:关于3种计算张量函数方法的比较

109

φ1=116498×103, φ2=1551967.

所得结果为

1019662

(e C )

ij

tion of Fourth -Order Isotropic Tensor Functions [J ].Comput. Struct. , 1998,66:37-43.

[3]Luehr C P , Rubin M B. The Significance of Projection Op 2

erators in the Spectral Represention of Symmetrical 2nd -Op 2

3

[***********]9

[1**********]719×10. 201=[1**********]96

erators[J].C omput. Method Appl. M. ,1990,84:243-246. [4]Ortiz M , Radovitzky R A , Repetto E A. The Computation

of the Exponential and Logarithmic Mappings and Their First and Second Linearizations[J].Int. J. Numer. Meth 2ods Engrg. ,2001,52:1431-1441.

[5]Itskov On the Theory -Order Tensors and

Their ].C om p. Mech. :419-438. ]of Isotropic Tensor Functions[J].

Numer. Methods Engrg. ,1993,9:889-896. [7]黄筑平. 连续介质力学基础[M ].北京:高等教育出版

　　采用最佳计算特征根和特征向量的方法, 将运算循环100万次, 每种方法计算张量函数所花费的时间分别为:方法①,418s ; 方法②,716s ; 方法③, 313s. 结果表明, 表示定理法计算张量函数所用的

时间仅占特征投影算子法所用时间的43%, 比任何一种方法计算特征根的时间要少.

4　结论

, 对, 而没有必要进行数值迭代求解, 从而大大提高了计算速度. 通过利用各向同性对称张量函数的3种计算方法(谱分解法、特征投影算子法和表示定理法) 对指数张量函数进行的实例计算表明, 表示定理法所用的时间最短, 谱分解法其次. 另外还可以发现, 在整个计算过程中, 计算特征根所用时间占的比重较大. 所以, 采用哪种方法计算特征根是非常重要的. 参考文献:

[1]Hill R. On the Constitutive Inequalities for Simple Materi 2

als -I[J].J. Mech. Phy. S olids. , 1968,16:229-242. [2]Miehe C. Comparison of Two Algorithms for the Computa 2

社, 2003.

HUANG Zuping. Fundamentals of C ontinuum Mechanics[M].Beijing :Higher Education Publisher ,2003. (in Chinese ) [8]兑关锁, 王正道, 金明. 各向同性张量函数的导数[J].

中国科学G ,2006,36(1) :89-102.

DU I Guansuo , WAN G Zhengdao , J IN Ming. Derivatives on the Isotropic Tensor Functions [J].Science in China :Series G , 2006,49(3) :321-334. (in Chinese )

[9]Dui G , Ren Q , Shen Z. Time Rates of Hill ’s Strain Ten 2

sors[J].J. Elast. ,1999,54:129-140.

[10]Dui G , Chen Y C. Basis 2Free Representations for the

Stress Rate of Isotopic Materials [J ].Int. J. S olids Struct. ,2004,41:4845-4860.

[11]Dui G. Determination of the Rotation Tensor in the Polar

Decomposition[J].J. Elast. ,1998,50:197-208.