Luận án Nghiên cứu tính toán ống trụ composite chịu tác dụng của áp suất di động

 
thông qua các đạo hàm theo các biến 1, 2 như sau: 
 NN  
 Ni   1  2 ii
 r  r  r 11  
  J  (2.23) 
 NNNi  1  2  i  i 
 s   s  s 22    
 12 
 rr JJ11 12
với: J - ma trận Jacobi 
 12  JJ21 22
 ss
 Ni
 r
  - đạo hàm riêng hàm dạng theo (r, s). 
 Ni
  s
 Từ đây ta nhận được: 
 36 
 N
 i NNii  
 ** 
 1 1 rr JJ11 12 
  J   (2.24) 
 NNN ** 
 i i JJ21 22 i 
  2 ss   
2.4.1.2. Chuyển vị của 1 điểm trong phần tử: 
 Vectơ chuyển vị của 1 điểm bất kỳ trong phần tử: 
 u u v wT (2.25) 
 Vectơ chuyển vị và góc xoay tại điểm bất kỳ trên mặt trung bình: 
 T
 u u v w   (2.26) 
  0 0 0 12
 Từ (2.1) ta có quan hệ 2 vectơ chuyển vị u và u: 
 1 0 0 0 
 u 0 1 0  0 u L u (2.27) 
    
 0 0 1 0 0
 với: L là ma trận cỡ (3x6) trong (2.27). 
 Khi chưa kể đến bậc tự do xoắn  , vectơ chuyển vị tại nút thứ i của 
 i
 phần tử là: 
 T
 q u v w   với i = 1÷4 (2.28) 
 i 0i 0i 0i 1i 2i 
 Vectơ chuyển vị của phần tử có kích thước 20x1: 
 e T
 q q1 q 2 q 3 q 4 (2.29) 
 Suy ra chuyển vị của một điểm bất kỳ trên mặt trung bình thuộc phần 
tử được biểu diễn theo hàm dạng và vectơ chuyển vị nút như sau: 
 37 
 e
 u  Nu  q , (2.30) 
trong đó: Nu  là các ma trận hàm dạng và có cấu trúc như sau: 
 NNINININI,u 1 5 2 5 3 5 4 5  (2.31) 
với: [I5] là ma trận đơn vị kích thước 5 5. 
2.4.1.3. Tọa độ của điểm bất kỳ trong phần tử: 
 Trong hệ tọa độ cục bộ, tọa độ (1, 2, ) tại điểm bất kỳ trong phần tử 
có thể biểu diễn thông qua hàm dạng Ni , tọa độ các nút của phần tử trong 
hệ tọa độ (1i, 2i, i) và chiều dày vỏ tại các nút i, với i = 1÷4 [1], [38]: 
 1   1i  l 3i
 4 1 
 2  N i  2i  t.h i m 3i  (2.32) 
 i1 2 
    i  n 3i
trong đó: hi là chiều dày phần tử tại nút i; l3i, m3i, n3i là cosin chỉ phương 
của vectơ pháp tuyến đơn vị k tại nút i trong hệ tọa độ cục bộ. 
2.4.1.4. Ma trận ma trận cosin chỉ phương của phần tử: 
 Ma trận cosin chỉ phương giữa hệ tọa độ tổng thể ( ) và hệ tọa 
độ cục bộ ( ) có dạng: 
 cos(e , i) cos(e , i) cos(e , i)
 l1 l 2 l 3 1 2 3
 L m1 m 2 m 3 cos(e 1 , j) cos(e 2 , j) cos(e 3 , j) (2.33) 
 e1 ,e 2 ,e 3
 n n n 
 1 2 3 cos(e1 ,k) cos(e 2 ,k) cos(e 3 ,k)
 i, j,k
 Các cosin chỉ phương tại điểm bất kỳ trên mặt trung bình của phần tử 
theo các biến trong hệ tọa độ tự nhiên (r,s) được xác định như sau: 
 38 
 Cosin chỉ phương của vectơ pháp tuyến e3 tại điểm bất kỳ có tọa độ 
(r,s) được xác định bằng công thức: 
 11   
 rs 
   
 22  
 rs 
   
 l3 
 rs (r,s)  (r,s)
 em33  (2.34) 
 11   
 n3
 (r,s) rs 
   
 22  
 rs 
   
 rs (r,s)  (r,s)
 Viết dưới dạng định thức ma trận với các vectơ đơn vị trong hệ 
tọa độ cục bộ: 
 i j k
               
e 1 2 ( 2 2 )i ( 1 1 ) j ( 1 2 2 1 )k
 3 rrr  rsrs  rsrs  rsrs 
 12  
 s  s  s
 a3 i b 3 j c 3 k 2.35 
 Các thành phần của vectơ hay các cosin chỉ phương của nó trong hệ 
tọa độ cục bộ được xác định bằng công thức: 
 i, j,k
 a3 b 3 c 3
 l ; m ; n (2.36) 
 32 2 2 3 2 2 2 3 2 2 2
 a3 b 3 c 3 a 3 b 3 c 3 a 3 b 3 c 3
 39 
 Cosin chỉ phương của vectơ pháp tuyến e2 được xác định bằng tích có 
hướng của vec tơ e3 và vec tơ đạo hàm theo biến r tại tâm phần tử (r = s = 0): 
 
 1
 r
 ll23  
 1  
 e m  m  2  (2.37) 
 2 2 3 r
 1 
 nn23 
 (r,s) (r,s) 
 l r
 3 r
 2 (r,s)
 m  
 3 r
 n 
 3 (r,s) 
 r (0,0)
 Viết dưới dạng định thức ma trận với các vec tơ đơn vị trong hệ 
tọa độ cục bộ: 
 i j k
      
 e l mn (m n2 )i(n 1 l)j(l 2 m 1 )k
 2 3 3 3 3r 3  r 3  r 3  r 3  r 3  r 
 12  
 r  r  r
 a2 i b 2 j c 2 k (2.38)
 Các thành phần của vectơ hay các cosin chỉ phương của nó trong hệ 
tọa độ cục bộ được xác định bằng công thức: 
 a b c
 l 2 ; m 2 ; n 2 (2.39) 
 22 2 2 2 2 2 2 2 2 2 2
 a2 b 2 c 2 a 2 b 2 c 2 a 2 b 2 c 2
 i, j,k
 Cosin chỉ phương của vectơ pháp tuyến e1 được xác định bằng tích có 
hướng của vec tơ và vec tơ : 
 40 
 l1  l 2   l 3
 e1 m 1  m 2  m 3  (2.40) 
 n n n 
 1 (r,s) 2  (r,s)  3 (r,s)
 Viết dưới dạng định thức ma trận với các vec tơ đơn vị trong hệ 
tọa độ cục bộ: 
 i j k
 e12 l m 22 n (mn 2323 nm)i (nl 2323 ln)j (lm 23 ml)k 23
 (2.41) 
 l3 m 3 n 3
 a1 i b 1 j c 1 k
 Các thành phần của vectơ hay các cosin chỉ phương của nó trong hệ 
tọa độ cục bộ được xác định bằng công thức: 
 a b c
 l 1 ; m 1 ; n 1 (2.42) 
 12 2 2 1 2 2 2 1 2 2 2
 a1 b 1 c 1 a 1 b 1 c 1 a 1 b 1 c 1
 Các cosin chỉ phương (2.36), (2.39), (2.42) phụ thuộc biến r và s trong 
hệ tọa độ tự nhiên. Khi tính cosin chỉ phương tại nút i lấy giá trị ri và si tại nút 
i của phần tử. 
2.4.1.5. Phương trình vi phân dao động phi tuyến của phần tử vỏ: 
 Áp dụng nguyên lý Hamilton cho phần tử [26] ta có: 
 t1
  Te U e W e dt 0 i, j,k (2.43) 
 e1
 t0
 41 
 eeeee e
trong đó: T U W He q , q ,t là hàm tác dụng Hamilton, T 
và Ue là động năng và thế năng biến dạng của phần tử, We là công ngoại lực 
tác động lên phần tử, qe là vectơ vận tốc nút phần tử. 
 Trường hợp không kể đến lực cản, phương trình (2.43) dẫn đến: 
 d HH
 ee 0, (2.44) 
 dt qqee
   
 Động năng của phần tử được xác định theo biểu thức [1], [12]: 
 1 T
 Te u u dV (2.45) 
 2 
 Ve
trong đó: 
 u 
 u 0 2
 u v v  L u : vận tốc của chất điểm; 
   0 1  
 w 
  w0
 - khối lượng riêng vật liệu tại vị trí của chất điểm. 
 Vậy động năng của phần tử vỏ có thể biểu diễn như sau: 
 1 T T T
 e e LL e
 T q Nu  N u dV e q  (2.46) 
 2 
 Ve
 Thế năng biến dạng toàn phần của phần tử [12],[26]: 
 42 
 11TT
 Ue   dV    dV
 22 b b s s
 VVee
 (2.47) 
 11TT
  [Q ]  dV  [Q ]   dV
 22 b b b s s s
 VVee
 Thay {b},{s} từ ( 2.7 ) và (2.8) vào (2.47) ta có: 
 T
 e1 L 1 N  L 1 N 
 U [D]b D u  Q b [D] b D udV  
 2 2  2 
 V
 e (2.48) 
 1 T
 Ds u Q s  D s  u dV
 2 
 Ve
 Thay biểu thức (2.30) vào (2.48) dẫn đến: 
 T
 T 1 1 1
 e e L N  L N  e
 U q [Bb ] B  Q b [B b ] B  dV e q  
 2 2  2 
 Ve (2.49) 
 T
 qee  BT Q B dV q ,
  s s s e 
 Ve
trong đó: 
 BDNBBBB,LLLLLL (2.50) 
 b b u b1 b2 b3 b4
 L
 [B]bi  [B] i {B Ri } [B], i  (2.51) 
 Bi  D i I Ni , {B Ri } {D Ri } I Ni  ,
 N0i (2.52) 
 BDI,I,i  i Ni  Ni  
 0Ni
 BDNBBBB,s  s u  s1  s2  s3  s4  (2.53) 
 B 0ww B N I , B D N , (2.54) 
  si   i i s i N i
 43 
 BDNBBBB,NNNNNN (2.55) 
  u 1 2 3 4
 Công ngoại lực do tải trọng ngoài tác động [26]: 
 TTT
 We qe fdV e q e fdA e q e f, e
  b e  s e  c  (2.56) 
 VAee
 e
trong đó: Ae, Ve tương ứng là diện tích và thể tích phần tử, fb - vectơ lực 
 e f e
khối, fs - vectơ lực bề mặt phần tử, c - vectơ lực tập trung phần tử. 
 Thay (2.46), (2.49), (2.56) vào (2.44) dẫn đến phương trình vi phân mô 
tả dao động không cản của phần tử vỏ trong hệ tọa độ cục bộ: 
 e e e e e
 M q K q f , (2.57) 
 EE   E
 e e
trong đó: q  là vectơ gia tốc nút phần tử, M E là ma trận khối lượng 
 e e
của phần tử, K E là ma trận độ cứng của phần tử; f  E là vectơ lực nút 
của phần tử trong hệ tọa độ cục bộ do tải trọng động gây ra, chúng được 
xác định như sau: 
 - Ma trận khối lượng của phần tử: 
 T
 Me NT L N L dV ,
 E  u   u e (2.58) 
 Ve
 T
 Đặt: I L L d  
  0  
 h
 Me NT I N dA
 Khi đó: E  u   0u  
 Ae
 44 
 11
 e T
 Hay: M N I N J drds 
 E  u  0 u 
 11
trong đó: J là định thức của ma trận Jacobi được xác định theo (2.24). 
 - Ma trận độ cứng của phần tử: 
 e ue Q e u e
 K, KKK (2.59) 
 E LN E EE 
 - Vectơ lực nút của phần tử: 
 e e e e
 f = fb dV+ f s dA+f c  , (2.60) 
 E 
 VAee
với: 
 + Ma trận độ cứng uốn tuyến tính: 
 11
 ue L T L L T L
 KL B b QBdV b b B b  ABJdrds b (2.61) 
 E 
 Ve 1 1
 + Ma trận độ cứng cắt: 
 11
 Q e TT
 K  BHBSSSS   dV BHB   Jdrds (2.62) 
 E 
 Ve 1 1
 + Ma trận độ cứng uốn phi tuyến: 
 T
 LNNL 
 e [Bbb ] B B B B [B ]
 u 
 KN  dV
 E T
 BDBNN 
 Ve    
 (2.63) 
 T
 11 LNNL 
 [Bbb ] B B B  B [B ] 
 J drds
 T
 BDBNN 
 11    
 45 
 Khi kể đến bậc tự do xoắn  (i = 14), vectơ chuyển vị nút của phần 
 i
tử có dạng: 
 T
 e
 q u v w  u v w  
  1 1 111  21 4 4 4  14  24  1  2  3  4
 Lúc này, ma trận khối lượng, ma trận độ cứng và vectơ tải trọng nút 
của phần tử vỏ cong mỗi nút có 6 bậc tự do như sau: 
 e
 u Q u 
 e KLNN K k K 0
 K, (2.64) 
 0K r
 24 24E E
 eT
 e 
 e M0   e e
 M , f f 0 , (2.65) 
 E   
 00   
 24 24EE 1 20 1 4
trong đó: kN - hệ số phi tuyến, nhận giá trị 1 khi tiến hành giải bài toán phi 
tuyến và nhận giá trị 0 khi giải bài toán tuyến tính. Các thành phần độ cứng 
xoắn quanh trục  của ma trận [Kr] nhận giá trị rất nhỏ, khi xét trong tổng 
 u Q u
thể của ma trận KLNN K k K , chúng nhận giá trị vào 
 E
khoảng 1/1000 so với giá trị phần tử lớn nhất trong ma trận trên. Điều này 
 -3
có nghĩa là kru(i,j) = 10 max(k(m,n)), với k(m,n) là các phần tử của ma trận 
 u Q u
 KLNN K k K . 
 E
 Lúc này phương trình dao động của phần tử vỏ: 
 eeeee
 M q K q f , (2.66) 
 EE   E
 46 
 e
 Do ma trận độ cứng của phần tử K trong phương trình (2.66) chứa 
 E
 e e
vectơ chuyển vị nút K K q e , nên (2.66) là phương trình vi 
 E  
 E
phân mô tả dao động phi tuyến của phần tử. 
2.4.2. Phần tử vỏ trong nền đàn hồi chịu tác dụng của áp suất di động 
 Do ống trụ composite được bao bọc bởi nền đàn hồi, nên xét trên 
phương diện phần tử thì mỗi phần tử được xem là đặt trong nền đàn hồi. 
 Xét phần tử vỏ nằm trong nền đàn hồi với mô hình nền Winkler, chịu 
tác dụng của áp suất phân bố đều di chuyển trên phần tử với vận tốc đều v, 
theo chiều trục 1 (Hình 2.4). 
 Hình 2.4. Mô hình phần tử vỏ cong trong nền đàn hồi chịu áp suất di động 
 Xét phần tử bất kỳ tại thời điểm t. Trong trường hợp tổng quát, trên 
phần tử tồn tại hai vùng diện tích A1 và A2, với A1 là diện tích chịu áp suất 
phân bố đều theo cường độ p(t), còn A2 là diện tích mà trên đó áp suất bằng 
0 (do môi chất chưa di chuyển đến). Các diện tích A1, A2 xác định bởi: 
 47 
  1  le
 A R  d  d  ; A R  d  d  (2.67) 
 1 1 1 2 1 1
 00 0 1
Nếu phần tử hoàn toàn nằm trong vùng chịu áp suất thì A2 = 0. Còn nếu phần 
tử hoàn toàn nằm trong vùng không chịu áp suất thì A1 = 0. 
 Từ đó ta tính vectơ lực nút [12]: 
 e TT
 f  Nu p(t)dA N u p f dA
 E 
 AA1e
 NTT p(t)dA N k w t dA (2.68) 
 u u 0
 AA1e
 1 le
 R NTT p(t)dd   R N kNdd    q,e
 u 1 u 0 u 1
 0 0 0 0
trong đó: 
 pf = k0w(t) - phản lực nền; k0 - hệ số nền; w(t) – chuyển vị hướng 
kính;  - góc mở của phần tử; R – bán kính ống trụ composite. 
 Ta đặt: 
  1  le
 e T e T
 fm R N u p(t)d  d  1 , Kf R N u k 0 N u d  d  1 
 E E 
 00 00
và thay (2.68) vào (2.66) ta có phương trình mô tả dao động phi tuyến của 
phần tử vỏ trong nền đàn hồi chịu tác dụng của áp suất di động: 
 eee e e e
 M q KK q f . (2.69) 
 EEE  f mE
 Các ma trận khối lượng, ma trận độ cứng và vectơ tải trọng phần tử 
trong (2.69) được chuyển về hệ tọa độ tổng thể 0xyz thông qua ma trận 
 48 
chuyển hệ tọa độ [Te], thông qua các công thức sau [1], [12], [14]: 
 e TTeee
 MMK TT,KTT,e e e e 
 G EE G 
 (2.70) 
 e e T e
 fmm T f .
 G E
 e e e
trong đó M , K , fm tương ứng là ma trận khối lượng, ma trận độ 
 G G G
cứng và vectơ tải trọng nút của phần tử trong hệ tọa độ tổng thể. 
Ma trận chuyển hệ tọa độ [Te] được xác định: 
 L . 0
 e L.
 T (2.71) 
 ....
 0 . L
     (6x6)
 với: L là ma trận cosin chỉ phương giữa các trục tọa độ của hệ trục tọa 
độ cục bộ của phần tử và hệ tọa độ tổng thể, được xác định bởi (2.33). 
Do đó, phương trình vi phân mô tả dao động của phần tử vỏ trong nền đàn 
hồi chịu áp suất di động (2.69) viết trong hệ tọa độ tổng thể như sau: 
 ee e e e e
 M q K Kf q fm (2.72) 
 G  G
 G
2.5. Thuật toán giải phƣơng trình vi phân dao động phi tuyến của 
ống trụ trong nền đàn hồi chịu áp suất di động 
2.5.1. Phương trình vi phân dao động phi tuyến của ống trụ trong nền 
đàn hồi chịu áp suất di động 
 Phương trình tổng thể mô tả dao động của hệ được xuất phát từ 
phương trình mô tả dao động của phần tử (2.72) bằng cách tập hợp các ma 
 49 
trận và vectơ tải trọng phần tử theo phương pháp độ cứng trực tiếp, trong 
đó áp dụng nguyên tắc ma trận chỉ số và sơ đồ Skyline, đây là phương pháp 
phổ biến trong lý thuyết PTHH. Theo phương pháp này, ma trận độ cứng 
và ma trận khối lượng của hệ được hình thành theo các bước sau: 
 - Tạo ma trận rỗng kích thước mn nd; 
 - Khởi tạo ma trận chỉ số, xác định vị trí nút của phần tử, tương ứng 
với chỉ số bậc tự do; 
 - Tại vị trí nút chung, theo tính chất bậc tự do của các phần tử chung 
nút, phần tử có chung tác dụng của cùng một loại bậc tự do được góp vào 
ma trận chung của hệ bằng phương pháp cộng trực tiếp, ngược lại chúng 
được cộng giá trị 0. 
 Tập hợp ma trận, vectơ tổng thể từ ma trận, vectơ phần tử: 
 ne
 e
 MM, []G
 e1 
 nnee
 ee
 K, [K]G [K]fG
  (2.73) 
 e 1 e 1
 nem
 f [f ]e ,
 m  m G
 e 1
trong đó: ne – số phần tử vỏ mô phỏng ống, nem – số phần tử vỏ chịu áp 
suất di động tác dụng. 
 Sơ đồ phương pháp mô tả các công thức (2.73) được thể hiện như sau: 
 50 
 Sau khi tập hợp các ma trận và vectơ tải trọng tổng thể, ta có phương 
trình vi phân dao động phi tuyến không cản của ống trụ composite trong 
nền đàn hồi chịu tác dụng của áp suất di động: 
 M q K q f
   m (2.74) 
 Khi xét đến cản, giả thiết lực cản của hệ tỷ lệ với vận tốc chuyển vị: 
 f Cq
 d  , thay vào (2.74) ta có phương trình vi phân dao động của 
hệ như sau: 
 M q C q K q f,
    m (2.75) 
với: C - ma trận cản của hệ. 
 Thông thường, đối với hệ nhiều bậc tự do, việc xác định ma trận cản 
tổng thể của hệ theo ma trận cản phần tử là không thể, do khó xác định 
được tỉ số cản phần tử. Lúc này, phương pháp cản Rayleigh thường được 
sử dụng. Theo đó, ma trận cản C được tính thông qua ma trận khối 
lượng M và ma trận độ cứng K : 
 51 
 CMK  (2.76) 
 RR 
trong đó R, R là các hằng số cản Rayleigh được xác định thông qua tỷ số 
cản R và hai tần số dao động riêng đầu tiên 1,2 của hệ: 
 22RR
 R ;. R  1  2  R  1  2 (2.77) 
 1  2  1  2
 Trường hợp ống chịu áp suất di động, điểm khác biệt so với ống chịu 
tải trọng thông thường có điểm đặt lực không đổi là: 
 - Tại mỗi thời điểm chỉ có một số phần tử cùng với một phần của các 
phần tử chịu tác dụng của áp suất, các phần tử còn lại và một phần của các 
phần tử trước đó không chịu áp suất tác dụng. Hình 2.5 thể hiện ống có 
chiều dài L được chia ra thành n phần tử hữu hạn. 
 - Trong tính toán, tại mỗi thời điểm cần phải xác định phần tử chịu áp 
suất tác dụng và phần diện tích của phần tử chịu áp suất để xác định vectơ 
tải trọng vế phải của phương trình (2.75). Quá trình này được lặp theo bước 
thời gian cho đến khi toàn bộ ống đều chịu áp suất tác dụng, do giải lặp 
theo bước thời gian nên cho phép giải bài toán với các trường hợp áp suất 
p(t) là hằng số hoặc là hàm phụ thuộc thời gian. 
 52 
 Hình 2.5. Áp suất tác dụng vào ống tại thời điểm t 
2.5.2. Điều kiện biên và phương trình dao động của hệ sau khi khử biên 
 Với bài toán ống trụ composite đặt trong nền đàn hồi tựa hai đầu chịu 
áp suất di động như trong luận án, điều kiện biên được xác định: 
 w 0,t w L,t 0,
 (2.78) 
 M 0,t M L,t 0.
trong đó: w(x,t) - chuyển vị hướng kính tại điểm x và thời gian t, 
 M(x,t) - mô men uốn tại điểm x và thời gian t. 
 Với ống trụ composite đặt trên gối cứng (liên kết tựa): chuyển vị 
hướng kính tại các nút chứa liên kết tựa bằng 0. 
 Điều kiện biên được cụ thể hóa trong phương pháp PTHH theo 
phương pháp khử biên, nội dung cụ thể được thực hiện với tính chất của 
bậc tự do trên biên. Theo đó, tùy theo các loại liên kết sẽ biết được tính 
chất của các bậc tự do tại biên (ví dụ: tại ngàm, tất cả các bậc tự do q0i 
và theo đó hàng thứ i, cột thứ i trong hệ phương trình (2.75) bị xóa và do 
đó số phương trình và số ẩn số trong hệ phương trình (2.75) sau khi khử 
 53 
biên sẽ bé hơn khi chưa khử biên. Giả sử trước khi khử biên, hệ (2.75) có 
mn nd phương trình và mn nd ẩn số, hệ có biên ngàm chứa tất cả m nút, 
mỗi nút có n bậc tự do, dẫn đến có m×n bậc tự do bằng 0, vì vậy hệ (2.75) 
sau khi khử biên còn (mn nd - m×n) phương trình và (mn nd - m×n) ẩn số. 
 Sau khi khử biên như trên, ta có phương trình vi phân dao động của 
ống trụ composite trong nền đàn hồi chịu áp suất di động được biểu diễn 
bằng phương trình: 
 M q  C q  K q f.m (2.79) 
2.5.3. Thuật toán PTHH giải phương trình dao động của hệ 
 a) Bài toán dao động riêng tuyến tính, được mô tả bởi phương trình: 
 M q KL q {0}, (2.80) 
trong đó: KL - ma trận độ cứng tuyến tính của hệ ống - nền. 
 Từ phương trình (2.80) các tần số riêng của hệ được xác định bởi 
phương trình: 
 KML2   0 (2.81) 
với:  là tần số riêng của hệ. 
 Các dao động riêng {qi} của hệ được xác định theo các vectơ riêng 
{Qi} bởi phương trình: 
 L2
 KM ii  Q0  (2.82) 
 b) Bài toán dao động cưỡng bức phi tuyến: 
 54 
 Nhiệm vụ ở đây là xác định đáp ứng phi tuyến động lực học của hệ, 
bằng cách giải hệ phương trình phi tuyến (2.7). Trong trường hợp này 
phương pháp tích phân trực tiếp Newmark kết hợp với phương pháp lặp 
Newton-Raphson được tác giả áp dụng. Theo đó, nghiệm của phương trình 
(2.79) tại bước lặp thứ i, ở thời điểm tính t + t được xác định bởi phương 
trình [1]: 
 (i) (i) (i)
 M q C(i) q K (i) q f (i) (2.83) 
  tt   tt  tt   tt  tt  mtt 
 Điều kiện ban đầu cho mỗi cấp tải trọng được xác định như sau: 
 v(0) = v0, p(0) = p0. 
 q 0 q;f0 =f;K 0 K. (2.84) 
 t t  t m t+Δt t t t  t 
 Ma trận độ cứng tiếp tuyến hiệu quả được xác định: 
 (i) ii 
 K* K a M a C (2.85) 
 t t t t 0  1 t t  
 Ma trận vectơ tải trọng hiệu quả được tính theo công thức sau: 
 * (i) i-1 i 1 i 1 i 1 
 fm =f m +Maq  0tt  aq 2tt  aq 3tt  
 t+Δt t+Δt 
 i 1 i 1 i 1 
 C a q a q a q  2.86 
 tt1tt 4tt 5tt 
 (i)
 trong đó: K - ma trận độ cứng tiếp tuyến, 
 tt 
 (i)
 fm - Vectơ lực nút quy đổi.
 t+Δt
 Chuyển vị, vận tốc và gia tốc nút tại thời điểm t + t của bước lặp thứ i là: 
 (i) (i)
 (i) *
 Kt t q t t f m , (2.87) 
 t+Δt
 55 
 (i) (i1) (i1) (i)
 qt t q t t a 6 q t t a 7 q t t  (2.88) 
 (i) (i) (i1) (i1) (i1) 
 qtt  a 0tt q  q tt  a 2tt q  a 3tt q  
với: 
 1 1 1
 a ; a ; a ; a 1;
 02 1 t 2 t 3 2 
 t (2.89) 
  t 
 a 1;a ( 2);a  t(1 );a  .t
 4 52 6 7
các tham số  và được chọn theo [1]. 
 t
với  0,5; 0,25, điều kiện ổn định nghiệm có dạng: . Như vậy 
 Tn
điều kiện ổn định nghiệm luôn thỏa mãn với mọi giá trị của t. Tuy nhiên 
để đảm bảo độ chính xác, giá trị t phải đủ nhỏ, thông thường 
 11
 tT  1 , với T1 là chu kỳ dao động đầu tiên của hệ. 
 10 8
 Tiêu chuẩn dừng của phép lặp là sự hội tụ của chuyển vị nút, được xác 
định bởi điều kiện: qc  , 
 q(i)
 Hay:   , (2.90) 
 qc(i)
 qqt t t 
 -4 -3
với: c - độ chính xác theo chuyển vị, thường được chọn 10  10 . 
 Tóm tắt thứ tự các bước của thuật toán như sau: 
 1. Khai báo các thông số ban đầu của bài toán. Tính các ma trận [M], 
[K], giải bài toán dao động riêng, xác định ma trận cản [C]. 
 56 
 2. Cho tổng thời gian tính t, chọn bước thời gian t , các tham số , ; 
tính các hệ số a0 , a1 , a2 , a3 , a4 , a5 , a6 , a7 . Gán giá trị đầu tại thời điểm 
t = 0 cho chuyển vị và vận tốc. 
 3. Thực hiện vòng lặp tích phân theo thời gian: j = 1, 2, , n tương 
ứng với các thời điểm t , 2t , nt với n t / t . Tại mỗi thời điểm t + 
 t xác định vị trí của áp suất (các phần tử chịu áp suất và phần diện tích 
các phần tử còn lại chịu áp suất tác dụng), tí