Phương pháp phần tử hữu hạn cho bài toán đàn nhớt tuyến tính tựa tĩnh

Tạp chí KHOA HỌC ĐHSP TP. HCM Trịnh Anh Ngọc 
 36 
PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN CHO BÀI TOÁN 
ĐÀN NHỚT TUYẾN TÍNH TỰA TĨNH 
Trịnh Anh Ngọc1 
1. Giới thiệu 
Bài toán tựa tĩnh của lý thuyết đàn nhớt tương tự như bài toán trong trường 
hợp đàn hồi ngoại trừ quan hệ ứng suất – biến dạng được thay bởi 
  ij ijkl kl
ijkl
kl
t
C t
C t s
s
s ds 
 
z( ) ( ( )) ( ) ( ( ))0 0u u , (1) 
trong đó u x ( ( , ))u ti là trường chuyển dịch,   ( ( , ))ij tx là trường tenxơ ứng 
suất,   ( ( , ))ij tx là tenxơ biến dạng, xác định từ chuyển dịch nhờ hệ thức 
 ij i j j iu u 
1
2
( ), , , (2) 
C x ( ( , ))C tijkl là tenxơ chùng ứng suất thỏa các điều kiện đối xứng 
 C Cijkl jikl , C Cijkl ijlk , C Cijkl klij . (3) 
Với vật liệu đàn nhớt ứng xử tức thời là đàn hồi [3], nghĩa là tồn tại hằng số 
c0 0 sao cho 
 C cijkl ij kl ij ij( )0 0    . (4) 
Để đơn giản cách viết, như trong phương trình (1), thường ta không ghi rõ 
sự phụ thuộc của các đại lượng vào biến không gian x và cả biến thời gian t nếu 
không gây ngộ nhận. 
Trong khoảng thời gian I T 0, , xét vật thể đàn nhớt tuyến tính chiếm 
miền dR (d=1,2,3) là tập mở bị chặn với biên    D N chính quy, giả 
thiết meas D( ) 0 . Lực tác dụng lên vật gồm: lực thể tích f x ( ( , ))f ti x  , 
t T [ , ]0 ; lực mặt g x ( ( , ))g ti , x N , t T [ , ]0 . Trên phần biên D vật được giữ 
cố định. Bài toán tựa tĩnh của lý thuyết đàn nhớt tuyến tính được phát biểu như 
sau: Tìm hàm u u x ( , )t thỏa 
  ij j if, trong  I , (5) 
0 u trong D I , (6) 
1 TS. – Trường ĐH KHTN TP. HCM 
Tạp chí KHOA HỌC ĐHSP TP. HCM Số 14 năm 2008 
 37 
 ij j in g trong N I , (7) 
0(0) u u trong  , (8) 
trong đó ứng suất  liên hệ với chuyển dịch u thông qua (1) và (2). 
Có nhiều phương pháp giải số bài toán biên tựa tĩnh. Một trong các phương 
pháp thông dụng là phương pháp đặt cơ sở trên phép biến đổi Laplace và nguyên 
lý tương ứng của lý thuyết đàn nhớt tuyến tính [5]. Trong một số trường hợp đặc 
biệt bài toán nhận được bằng phương pháp này có dạng tương tự bài toán của lý 
thuyết đàn hồi cổ điển, điều này thu hút sự chú ý của nhiều nhà tính toán số. Tuy 
nhiên, như đã biết, bài toán biến đổi ngược Laplace là một bài toán không chỉnh 
do đó độ chính xác của kết quả phụ thuộc vào rất nhiều yếu tố. Một cách tiếp cận 
khác là áp dụng phép rời rạc hóa theo biến không gian dựa trên phát biểu biến 
phân “nửa yếu”, bằng cách này bài toán dẫn về một hệ phương trình tích phân 
Volterra loại hai. Sau đó áp dụng phương pháp chọn điểm (collocation method) 
giải hệ phương trình tích phân này. S. Shawetal., trong [6], đã áp dụng phương 
pháp phần tử hữu hạn để xấp xỉ bài toán theo biến không gian, với thời gian tác 
giả xấp xỉ số hạng tích phân bằng phép cầu phương. Vấn đề sai số được tác giả 
dẫn theo [1]. Gần đây hơn, trong [7], S. Shaw đưa cách tiếp cận rời rạc cả không 
gian lẫn thời gian bằng phương pháp phần tử hữu hạn trên cơ sở công thức biến 
phân đầy đủ. Cách làm này dẫn đến một công thức xấp xỉ khác (với quy tắc cầu 
phương cổ điển) số hạng tích phân. Trong bài này chúng tôi áp dụng cách tiếp 
cận thứ hai để giải gần đúng bài toán (5)-(7) kếp hợp với (1), (2). Cách rời rạc 
hóa theo biến không gian tương tự như [6,7], nhưng ở đây, chúng tôi đặc biệt 
quan tâm đến cách xấp xỉ theo biến thời gian. Phép cầu phương được thực hiện 
bằng các công thức khác nhau, có chú ý đến sai số của công thức và sự tiện lợi 
khi cài đặt trên máy tính. Phần còn lại của bài được tổ chức như sau: Mục 2 trình 
bày công thức biến phân nửa yếu và bài toán xấp xỉ phần tử hữu hạn theo biến 
không gian. Mục 3 giới thiệu các công thức tích phân số để xấp xỉ bài toán trong 
Mục 2 theo biến thời gian. Một thí số được cho trong mục 4 để minh họa cách áp 
dụng và đánh giá (theo quan điểm thực hành) phương pháp. 
2. Rời rạc hóa theo biến không gian 
2.1. Phát biểu biến phân nửa yếu 
Ký hiệu 1 1( ) ( ( ))dH H không gian Hilbert với tích trong 
( , ) ( , )
( )
w v w vi i H1  
và chuẩn tương ứng   ( , ) . 
Tạp chí KHOA HỌC ĐHSP TP. HCM Trịnh Anh Ngọc 
 38 
Đưa vào không gian các hàm thử 
H v v { ( )H1 0  treân }D . 
Cố định t I , nhân phương trình (5) với v H bất kỳ, nhờ quan hệ ứng 
suất – biến dạng (1), sau một số biến đổi ta được bài toán phân nửa yếu: 
Tìm hàm u( ) ( , ; )t L T H 0 thỏa 
A t L t B t s s ds
t
( ( ), ) ( ; ) ( , ; ( ), )u v v u v z
0
 (9) 
với mọi v H , trong đó 
A C dijkl kl ij( , ) ( ) ( ) ( )w v w v z 0   

, (10) 
B t s
C t s
s
dijkl kl ij( , ; , )
( )
( ) ( )w v w v 
 
z    , (11) 
L t d d
N
( ; )v v f v g  z z 
 
. (12) 
Như một hệ quả của định lý 7.2, [2], tr. 189, ta có định lý sau. 
Định lý 1. Dưới các giả thiết: 
(i) Các thành phần tenxơ chùng ứng suất C tijkl ( ) là hàm đơn điệu giảm theo 
t ( t I ), có đạo hàm theo cấp một theo t thuộc L T L1 0( , , ( ))  ; hơn nữa, tồn tại 
hằng số c1 sao cho 
C t c Cijkl ij kl ijkl ij kl( ) ( )    1 0 (13) 
với mọi tenxơ đối xứng ( ) ij ; 
(ii) 1 2(0, ; ( ))L T L f và 1 2(0, ; ( ))NL T L g . 
Thì bài toán (10) tồn tại và duy nhất nghiệm. 
2.2. Bài toán xấp xỉ phần tử hữu hạn – Hệ phương trình tích phân 
Volterra loại hai 
Phân hoạch miền  thành ne phần tử hữu hạn tuyến tính (n-đơn hình) 
trong dR . Mỗi phần tử hữu hạn có 1d nút. Chuyển dịch nút 
 1{ , , }Tdq q q  ( 1, , 1d  ). (14) 
Vectơ chuyển dịch phần tử 
Tạp chí KHOA HỌC ĐHSP TP. HCM Số 14 năm 2008 
 39 
 { } { , , }q q q 1  T . (15) 
Các hàm dạng 
 N I  n , (16) 
trong đó  là tọa độ trọng tâm, dI là ma trận đơn vị cấp d . Ma trận hàm 
dạng của phần tử: 
 1 1[ ] [ , , ]d N N N . (17) 
Trong phần tử tương ứng, vectơ chuyển dịch u được xấp xỉ bởi 
 u N q [ ]{ } . (18) 
Vectơ biến dạng phần tử e [ , , , , , ]     11 22 33 23 13 12 T , 
 e u E q D [ ]{ } , (19) 
trong đó [ ] [ ]E N D với D là toán tử đạo hàm 
 D 



 
 
 
L
N
MMMMMMM
O
Q
PPPPPPP
1
2
3
3 2
3 1
2 1
0 0
0 0
0 0
0
0
0
. (20) 
Vectơ ứng suất phần tử s [ , , , , , ]     11 22 33 23 13 12 T , 
s D e D e
D E q D E q
 

 

z
z
[ ( )] ( ) [ ( )] ( )
[ ( )][ ]{ ( )} [ ( )][ ]{ ( )}
0
0
0
0
t t s
s
s ds
t t s
s
s ds
t
t (21) 
trong đó [ ]D là ma trận các hàm chùng ứng suất suy từ tenxơ C . 
Từ các công thức trên ta đưa vào ma trận độ cứng của phần tử e 
 [ ( )] [ ] [ ( )][ ]K E D Et t de eT e
e
 z 

. (22) 
Ma trận 
 [ ( )] [ ( )]G Kt s t s
s
e
e
 

 (23) 
Tạp chí KHOA HỌC ĐHSP TP. HCM Trịnh Anh Ngọc 
 40 
liên quan đến ảnh hưởng của lịch sử biến dạng, gọi là ma trận di truyền 
phần tử. 
Vectơ tải phần tử: 
 { } [ ] [ ]p N f N ge eT e eT ed d
e
N
e
 z z 
 
. (24) 
Bằng cách “lắp ghép” và từ (10) ta được bài toán xấp xỉ phần tử hữu hạn: 
Tìm hàm vectơ t t { ( )}q thỏa phương trình tích phân Volterra loại hai 
 [ ( )]{ ( )} { ( )} [ ( )]{ ( )}K q p G q0
0
t t t s s ds
t
 z , (25) 
trong đó [ ( )],[ ( )],{ },{ }K G q pt t lần lượt là các ma trận và các vectơ toàn cục. 
Định lý 2. Dưới các giả thiết của định lý 1, bài toán nửa rời rạc (25) tồn tại 
và duy nhất nghiệm. 
3. Rời rạc hóa theo biến thời gian 
3.1. Các công thức cầu phương 
Để giải phương trình (25) ta tính xấp xỉ tích phân ở vế phải phương trình 
 [ ( )]{ ( )}G qt s s ds
t
 z
0
, 
nghĩa là, rời rạc hóa phương trình này theo biến thời gian. Chia khoảng thời gian 
I thành nt khoảng con bằng lưới 
 { }0 1 2 1 1t t t t t Tn n nt  , t t ti i i 1 . (26) 
Có nhiều phương pháp xấp xỉ nhưng đơn giản hơn cả là quy tắc cầu phương 
[ ( )]{ ( )} (( )[ ( )]{ ( )} [ ( )]{ ( )})G q G q G qt s s ds t t t t t t t
t
t
j j j j j
j
j
 z 1 1 1 1   , (27) 
trong đó  [ , ]0 1 , các trường hợp đặc biệt: quy tắc Euler ( 0 ), quy tắc 
Euler lùi ( 1 ), quy tắc hình thang ( 1 2/ ). Theo [4], cấp của sai số trong 
quy tắc cầu phương với  0 và  1 là t , và bằng t 3 khi  1 2/ . 
3.2. Bài toán xấp xỉ không-thời gian 
Từ công thức (27), phương trình (25) có thể rời rạc thành bài toán: Tìm dãy 
các vectơ { ( )}q t j , j nt 1 2 1, , , sao cho, với mọi n nt { , , , }1 2 1 
Tạp chí KHOA HỌC ĐHSP TP. HCM Số 14 năm 2008 
 41 
[ ]{ ( )} { ( )} ( )[ ( )]{ ( )} [ ( )]{ ( )}( )K q p G q G qn n n j n j j n j j
j
n
t t t t t t t t t 
 1 1 1 1 1 1
1
1
1  n s
(28) 
trong đó 
[ ] [ ( )] [ ( )]( )K K Gn nt 0 0 . (29) 
Định lý 3. Dưới các giả thiết của định lý 1, với t ti max đủ bé, hệ 
phương trình (28) có nghiệm duy nhất. 
Công thức cầu phương (27) với  1 2/ có sai số bé hơn hai công thức còn 
lại. Tuy nhiên, về phương diện tính toán thì công thức (27) với  0 lại có lợi 
hơn do [ ] [ ( )]( )K Kn 0 với mọi n . Khi đó việc tính nghịch đảo (giải phương trình) 
chỉ thực hiện một lần là đủ, trong khi với hai công thức còn lại ở mỗi bước tính 
(xác định { ( )}q t j ) ta phải tính ma trận nghịch đảo [ ]( )K n 1 . Có thể thoát khỏi khó 
khăn này nếu bước thời gian là đều. 
3.3. Áp dụng 
Phương pháp xấp xỉ trình bày trong các mục trên được áp dụng cho bài toán 
phẳng đối xứng trục. 
4. Bài toán - Nghiệm giải tích 
Cho ống trụ đàn nhớt đẳng hướng tiết diện ngang hình vành khăn bán kính 
a b, ( a b ). Hệ tọa độ trụ có trục z trùng với trục ống. Dưới áp suất p t( ) trên 
thành trong,  rr a t p t( , ) ( ) , ống chỉ biến dạng theo hướng kính 
u r t u r t ur ( , ) ( , ),  0 . Thành ngoài giữ cố định u b t( , ) 0 . 
Trạng thái biến dạng:  rr r t
u
r
( , ) 

,  ( , )r t
u
r
 . Giả thiết vật liệu là đồng 
nhất đẳng hướng, quan hệ ứng suất – biến dạng cho 
   
  
   
  

rr
t
t
r t u
r
u
r
t s
s
u s
r
t s t s
s
u s
r
ds
r t u
r
u
r
t s
s
u s
r
t s t s
s
u s
r
ds
( , ) ( ) ( ( ) ( ))
( ) ( ) ( ( ) ( )) ( ) ,
( , ) ( ) ( ( ) ( ))
( ) ( ) ( ( ) ( )) ( ) ,


 

 



L
NM
O
QP


 



 

L
NM
O
QP
z
z
0 0 2 0
2
0 0 2 0
2
0
0
trong đó 
Tạp chí KHOA HỌC ĐHSP TP. HCM Trịnh Anh Ngọc 
 42 


( ) expt K t FHG IKJ23 ,  ( ) expt G
t
 FHG IKJ0 
Phương trình chuyển động tựa tĩnh: 


   rr rr
r r
0 . 
Bằng phương pháp biến đổi Laplace ta có nghiệm giải tích của bài toán 
(trường chuyển dịch): 
u r t p b r
Ka
G c
K G c
K
K G c
t( , ) ( ) ( )
( )
exp
( )
L
NM
O
QP
RST
UVW
2 2
0
2
0
2
0
22
1 1 3
1 3 1 3 
, 
trong đó c b a . Tính toán trực tiếp ta thu được ứng suất vòng   ( , )b t tại 
r b . 


( , )
( )
( )
exp
( )
b t p G c
K G c
K
K G c
t
L
NM
O
QP
RST|
UVW|1
1
1 3 1 3
0
2
0
2
0
2 . 
Các kết quả này được dùng để so sánh với giá trị xấp xỉ tìm được bằng 
phương pháp trình bày ở đây. 
4.1. Áp dụng phương pháp xấp xỉ 
Đưa vào không gian các hàm thử: 
H v v H a b v b { ( , ), ( ) }1 0 
với tích vô hướng 
( , ) ( ) ( )w v w r v r rdr
a
b
 z , 
ta có bài toán biến phân: tìm hàm u( )t H thỏa 
A t L t B t s s ds
t
( ( ), ) ( ; ) ( , ; ( ), )u v v u v z
0
với mọi v H . Ở đây 
A u v
r
v u
r
r u
r
v
r
uv
r
dr
a
b
( , ) ( ) ( ( ) ( )) ,w v 



FHG IKJ 


 FHG IKJLNM
O
QPz  0 0 2 0 
Tạp chí KHOA HỌC ĐHSP TP. HCM Số 14 năm 2008 
 43 
B t s u s v
r
v u s
r
t s t s
s
r u s
r
v
r
u s v
r
dr
a
b
( , ; , ) ( ) ( )
( ( ) ( )) ( ) ( ) ,
w v 



FHG IKJLNM
 





 FHG IKJOQP
z
 2
L t av a p t( ; ) ( ) ( )v . 
Công thức phần tử hữu hạn tuyến tính. Phần tử [ , ]r r1 2 
[ ] [ , ]N 1 2 1l
r r r r , l r r 2 1 , [ ] [ ] ( ) ( )
E N 
 L
NM
O
QP 
L
NM
O
QP
r
r l r r r r r r1
1 1 1
2 1
, 
[ ( )]
( ) ( ) ( )
( ) ( ) ( )
D t
t t t
t t t
L
NM
O
QP
  
  
2
2
, [ ( )] [ ] [ ( )][ ]K E D Et t rdrT
r
r
 z
1
2
, 
[ ( )] [ ] [ ( )][ ]G E D Et s
s
t s rdrT
r
r


 z
1
2
, 
{ ( )}
( )
p t
ap t
RST
UVW0 , nếu r a1 , và { ( )}p t 
RST
UVW
0
0
, nếu khác. 
4.2. Kết quả số – đánh giá 
Tính toán số được thực hiện với dữ liệu cụ thể: a b c b
a
 1 2 2, , ; 


( ) expt K t FHG IKJ23 ,  ( ) expt G
t
 FHG IKJ0 trong đó 
K G 2 3333 10769 10. , . , . 
Tính toán số theo quy tắc hình thang,  1 2/ , cho kết quả rất tốt. Với bước 
lưới không gian h 01. , bước thời gian t 01. sai số giữa nghiệm xấp xỉ và 
nghiệm chính xác chỉ cỡ 10 4 . Trên hình 1 là các đường cong chuyển dịch (c.x.) 
và nghiệm xấp xỉ (x.x.) của nó với t  0 10 20; ; . 
Kết quả tính toán khi dùng quy tắc Euler lùi,  1, tuy ổn định nhưng sai 
số ở những thời điểm lâu dài là khá lớn (cỡ 10 1 ) do sự tích tụ của sai số làm tròn 
(hình 2). Với quy tắc Euler,  0 , kết quả là không ổn định. 
Dùng công thức (21) có thể nhận được ứng suất xấp xỉ. Hình 3 cho kết quả 
xấp xỉ (theo quy tắc hình thang) ứng suất vòng tại r b ,   ( , )b t , so với ứng 
suất chính xác. Sai số cực đại cỡ 10 1 . 
Tạp chí KHOA HỌC ĐHSP TP. HCM Trịnh Anh Ngọc 
 44 
Hình 1: kết quả tính so với nghiệm chính xác theo quy tắc hình thang 
Hình 2: Kết quả tính so với nghiệm chính xác theo quy tắc Euler lùi 
Hình 3: Kết quả tính ứng suất   ( , )b t so với giá trị chính xác theo quy tắc hình thang 
Tạp chí KHOA HỌC ĐHSP TP. HCM Số 14 năm 2008 
 45 
TÀI LIỆU THAM KHẢO 
[1]. D.M. Bedivan and G.J. Fix (1997), Analysis of finite element 
approximation and quadrature of Volterra integral equations, Numer. Methods 
Partial Differential Equations, 13, pp. 663-672. 
[2]. G. Duvaut, J.L. Lions (1972), Les inéquations en mécanique et en physique, 
Dunod, Paris. 
[3]. Pipkin, A.C. (1972), Lectures on viscoelasticity theory, Springer – Verlag 
New York Inc. 
[4]. Mario G. Salvadori and Melvin L. Baron (1961), Numerical Methods in 
Engineering, Prentice-Hall International, London. 
[5]. Schapery, R.A. (1968), Stress analysis of viscoelastic composite materials, 
Composite Material Workshop, Technomic Publishing Co., Inc., pp.153-191. 
[6]. S. Shaw, M.K. Warby, J.R. Whiteman, C. Dawson, and M.F. Wheeler 
(1994), Numerical techniques for the treatment of quasistatic viscoelastic stress 
problems in linear isotropic solids, Comput. Methods Appl. Mech. Engrg., 118, 
pp. 211-237. 
[7]. S. Shaw and J.R. Whiteman (1999), Numerical solution of linear 
quasistatic hereditary viscoelasticity problems I: A priori estimates, BICOM: 
Tóm tắt 
Bài này trình bày một cách rời rạc hóa phần tử hữu hạn theo biến không 
gian và cầu phương theo biến thời gian cho bài toán đàn nhớt tuyến tính tựa tĩnh. 
Một thí dụ số được cho để minh họa cách áp dụng và thể hiện tính hiệu quả của 
phương pháp 
Abstract 
Finite element method for quasistatic linear viscoelasticity problems 
In this paper, we consider a finite-element-in-space, and quadrature-in-
time-discretization of a linear quasistatic viscoelasticity problem. A numerical 
application is presented in order to show validity of this discretization.