Phương pháp phân tích động phi tuyến kết cấu theo lịch sử thời gian không có điều kiện ổn định

Nội dung text: Phương pháp phân tích động phi tuyến kết cấu theo lịch sử thời gian không có điều kiện ổn định

1. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG PHƯƠNG PHÁP PHÂN TÍCH ĐỘNG PHI TUYẾN KẾT CẤU THEO LỊCH SỬ THỜI GIAN KHÔNG CÓ ĐIỀU KIỆN ỔN ĐỊNH GS. TS. SHUENN-YIH CHANG Trường Đại học Công nghệ Quốc lập Quốc gia Đài Bắc ThS. TRẦN NGỌC CƯỜNG Viện KHCN xây dựng Tóm tắt: Trong các phương pháp phân tích động hiển thức cho các bài toán cần quan tâm đến các phi tuyến theo lịch sử thời gian hiện nay, đã có một số dạng dao động bậc cao, nếu giá trị của các bước thời phương pháp không có điều kiện ổn định và có khả gian tính toán được lựa chọn để thỏa mãn điều kiện năng kiểm soát hệ số tiêu tán của hệ kết cấu. Tuy ổn định thì độ chính xác của kết quả cũng sẽ tự động nhiên, các phương pháp này đều là các phương pháp được thỏa mãn. nội ẩn thức, do đó quy trình tính toán đều yêu cầu tính lặp trong mỗi bước. Trong bài báo này, tác giả đề xuất Đã có một số các phương pháp tính tích phân phụ một họ phương pháp phân tích động phi tuyến mới. thuộc kết cấu được đề xuất bởi Chang [4, 5, 6, 7] Họ phương pháp này, tuy là ngoại hiển thức nhưng lại nhằm tích hợp đồng thời các ưu điểm của hai phương không có điều kiện ổn định. Phương pháp này còn có pháp nội ẩn thức và ngoại hiển thức, các phương hệ số tiêu tán thích hợp và có thể kiểm soát được, có pháp này đều có đặc điểm không có điều kiện ổn định thể điều chỉnh để hệ số cản nhớt số bằng không. Ưu (giống ưu điểm như phương pháp nội ẩn thức) và điểm lớn nhất của phương pháp này là không cần tính không cần tính lặp phi tuyến (giống ưu điểm của lặp trong mỗi bước, do vậy tiết kiệm được rất nhiều phương pháp ngoại hiển thức). Các phương pháp công sức tính toán so với các phương pháp nội ẩn này rất phù hợp để giải các bài toán thuộc dạng thứ thức hiện có. nhất. Tuy nhiên, khi giải các bài toán thuộc dạng thứ 1. Đặt vấn đề hai thì các phương pháp đó không có hệ số tiêu tán Về cơ bản, các bài toán về dao động của hệ kết (numerical dissipation) phù hợp để loại bỏ nhiễu do cấu được chia làm hai dạng chính: Dạng thứ nhất là ảnh hưởng bởi các dạng dao động bậc cao. Trong các hệ kết cấu bị chi phối bởi lực quán tính, chiếm đa các phương pháp được biết đến rộng rãi hiện nay, có số trong các bài toán động lực học công trình, dao một số họ phương pháp tính toán đã được phát triển động của dạng kết cấu này chủ yếu ảnh hưởng bởi có hệ số tiêu tán thích hợp, như các phương pháp các dạng dao động có tần số thấp; Dạng thứ hai là trong tài liệu [1, 10, 15, 16, 17], nhưng các phương các hệ kết cấu bị chi phối bởi cả các dạng dao động pháp này đều thuộc nhóm phương pháp nội ẩn thức, có tần số thấp và tần số cao, ví dụ như khi hệ kết cấu do vậy đều cần tính lặp phi tuyến khi giải các bài toán công trình bị tác động bởi lực gây ra do nổ và va kết cấu phi tuyến. chạm. Phương pháp để giải các bài toán thuộc dạng Vì những lý do trên, phương pháp được đề xuất thứ nhất thường được đề xuất là phương pháp nội ẩn trong bài báo này là một phương pháp ngoại hiển thức (implicit) như các phương pháp [1, 2, 3, 9, 11, thức không có điều kiện ổn định và không cần tính lặp 12, 13, 16]. Trong khi đó, ngoại hiển thức (explicit) là phi tuyến, giống các phương pháp đã được đề xuất phương pháp thích hợp để giải các bài toán dạng thứ trước đó [4, 5, 6, 7]. Điểm khác là phương pháp này hai, ví dụ như phương pháp của Newmark [13]. Điều còn có một hệ số tiêu tán phù hợp, có thể điều chỉnh, này là do các phương pháp nội ẩn thức không có điều dùng để giải các bài toán thuộc dạng thứ hai. kiện ổn định do vậy có thể sử dụng các bước thời 2. Phương pháp đề xuất gian (time step) lớn hơn, ngoài ra còn do các dạng dao động bậc cao không ảnh hưởng nhiều đến kết Phương trình dao động của hệ một bậc tự do quả bài toán. Ngược lại, ưu điểm chính của các được viết như sau: phương pháp ngoại hiển thức là các giá trị của bước mu cu  ku f (1) sau được tính trực tiếp từ bước trước mà không cần sử dụng các thuật toán tính lặp, thường khá phức tạp trong đó: trong các bài toán phi tuyến, do vậy khối lượng tính m, c, k, f lần lượt là khối lượng, độ cản nhớt vật lý, toán trong một bước sẽ ít hơn nhiều so với phương độ cứng và ngoại lực tác dụng lên hệ kết cấu; pháp nội ẩn thức. Khi áp dụng phương pháp ngoại Tạp chí KHCN Xây dựng - số 4/2015 3
2. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 2 u , u và u tương ứng là các giá trị chuyển vị, vận Để tiện cho việc tính toán, các hệ số ξΩ0 và Ω0 2 tốc, gia tốc. được viết lại ξΩ0 = ξω0(Δt) = c0(Δt)/(2m) và Ω0 = (Δt)2(k /m), theo đó, các hệ số β đến β và D trở Phương pháp đề xuất để giải phương trình dao 0 0 3 thành: động (1) được viết như sau: 3 2p p 1 2 p p 1 1p 1 2 2 ma c v k d k d f f  t k i 1 0 i 1 i 1 i 1 i i i 1 i 0 0 p 1 p 1 p 1 p 1 D 8 p 1 2 (2) d  d  d  t v  t a 3 i 1 0 i 1 1 i 2 i 3 i 1 p 1 2 2  1 t k 3p 1 p 3 1D8 p 1 0 v v t a t a i 1 i i i 1 2 p 1 2 p 1 1 p 3 (5)  m t c trong đó: 2 0 D 2 p 1 ai+1, vi+1, di+1, fi+1 là gia tốc, vận tốc, chuyển vị và 2  1 1 1 2 p 3 lực tác dụng lên hệ kết cấu ở bước thứ (i+1);  m t c  3D2 4 p 1 p 1 0  ai, vi, di, fi là gia tốc, vận tốc, chuyển vị và ngoại 3 lực tác dụng ở bước (i); p 3 p 2 2 D m t c t k 2 p 1 0 4 p 1 0 di-1 là chuyển vị nút ở bước tính toán thứ (i-1); ki, ki+1 là độ cứng của hệ ở bước (i) và (i+1); Có thể thấy rằng, các hệ số β0 đến β3 chỉ phụ thuộc vào độ cản nhớt c0 và độ cứng ban đầu k0 của c là độ cản nhớt vật lý của hệ kết cấu (giả thiết c 0 0 hệ kết cấu, nó không thay đổi trong suốt quá trình tính không thay đổi trong suốt quá trình tính toán, điều này toán, do vậy với mỗi bài toán chỉ cần tính duy nhất thường đúng với phần lớn các hệ kết cấu); một lần. Δt là bước thời gian tính toán. Trong dòng thứ hai của công thức (1), giá trị di+1 Các hệ số β0 đến β3 được tính như sau: được tính phụ thuộc vào các đặc điểm của kết cấu tính toán như khối lượng, độ cứng, độ cản nhớt, do đó 3 1p 1 2 2   phương pháp đề xuất ở đây gọi là phương pháp phụ 0 0 D 8 p 1 thuộc kết cấu. Nó khác so với phương pháp không phụ thuộc kết cấu [13], vì trong phương pháp Newmark, 3 1p 1 2 2  1  các hệ số β0 đến β3 đều là các hằng số cố định. Ngoài 1D8 p 1 0 (3) ra, giá trị di+1 khi tính toán không chỉ phụ thuộc vào các giá trị ở bước (i) trước đó mà còn phụ thuộc vào giá trị 1 p 3  1   ở bước (i-1), do vậy, khi áp dụng phương pháp này để 2D p 1 0 tính toán cần có một chút lưu ý khi tính bước đầu tiên, 2  1 1 1 2 p 3 điều này sẽ được nói rõ ở mục 4.     3D2 2 p 1 p 1 0  So sánh lời giải của phương pháp này so với lời giải của Zhou & Tamma [15, 16], mặc dù lời giải của trong đó: Zhou & Tamma cũng có độ chính xác tương tự ξ là hệ số cản nhớt; phương pháp này nhưng đó lại là phương pháp nội ẩn thức, khác với phương pháp đề xuất ở đây là Ω0 =ω0(Δt) với  k/ m là tần số dao động tự 0 0 phương pháp ngoại hiển thức. nhiên của hệ kết cấu tương ứng với độ cứng ban đầu k0; 3. Khảo sát tính chất của phương pháp đề xuất p là hệ số dùng để kiểm soát hệ số tiêu tán của phương pháp này (việc lựa chọn giá trị p sẽ được Khi khảo sát tính chất của một phương pháp phân trình bày rõ hơn ở mục 3). tích động phi tuyến áp dụng cho kết cấu công trình, các đặc tính thường quan tâm đến là khoảng ổn định Hệ số D được tính như sau: (stability), độ chính xác (accuracy) và tính hiệu quả 3 p 3 p 2 2 (efficiency) của phương pháp đó. Các tính chất này D 1    (4) p 10 4 p 1 0 được biểu hiện thông qua sai số tương đối của chu kỳ (relative period error), hệ số cản nhớt số (numerical 4 Tạp chí KHCN Xây dựng - số 4/2015
3. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG damping ratio) và sự biến thiên đột biến (overshooting). Sẽ rất khó để trình bày một cách chi tiết các bước để xây dựng lời giải này cũng như mô tả chi tiết việc khảo sát các tính chất của phương pháp tính chỉ trong khuôn khổ một bài báo, do vậy, ở đây chỉ nêu các đặc điểm chính, nếu bạn đọc quan tâm chi tiết hơn có thể tham khảo thêm các tài liệu [4, 5, 6, 7, 8]. Trong phần này, tác giả sử dụng một số thuật ngữ sau: Hình 1. Quan hệ giữa bán kính phổ và đại lượng Δt/T tương ứng với từng trường hợp p - NEM (Newmark Explicit Method): phương pháp ngoại hiển thức Newmark [13]; 3.2 Sai số tương đối của chu kỳ - AAM (Average Acceleration Method): phương Sai số tương đối của chu kỳ (relative period error) pháp nội ẩn thức gia tốc trung bình [13]; là đại lượng được định nghĩa bằng TTT / , trong đó T là chu kỳ dao động chính xác của hệ, T là chu - PFM1 (Proposed Family Method): phương pháp kỳ dao động tính toán của hệ. Đại lượng này đặc đề xuất với trường hợp p=1; trưng cho độ chính xác của phương pháp phân tích. - PFM2 (Proposed Family Method): phương pháp Đại lượng này càng nhỏ thì phương pháp phân tích đề xuất với trường hợp p= 0,5. càng chính xác. 3.1 Lựa chọn khoảng giá trị cho tham số p Hình 2 biểu diễn mối quan hệ giữa sai số tương Việc lựa chọn khoảng giá trị p cho phương pháp đối của chu kỳ của phương pháp đề xuất ở đây với các trường hợp p = 1; 0,75; 0,5. Sai số tương đối của này được căn cứ vào khoảng giá trị mà ở đó kết quả chu kỳ của phương pháp AAM cũng được thể hiện tính của phương pháp này là hội tụ. Muốn khảo sát trong hình vẽ để so sánh. các tính chất này, đầu tiên ta đi lập một ma trận đặc trưng rồi xem xét đến tính hội tụ (convergence), bán Hình 2 cũng cho thấy sai số tương đối của chu kỳ kính phổ (spectral radius). Một quy trình chung cho tỷ lệ nghịch với giá trị của p, p càng giảm thì sai số cách khảo sát này được trình bày trong tài liệu [1, 10, tương đối càng tăng, đồng nghĩa với độ chính xác của kết quả của phương pháp giảm xuống. Biểu đồ cũng 11] hoặc trình bày chi tiết hơn cho phương pháp này cho thấy rằng với các giá trị Δt/T nhỏ thì sai số cũng trong tài liệu [8]. nhỏ. Với giá trị Δt/T ≤ 0,1 thì sai số của phương pháp Kết quả khảo sát cho thấy với phương pháp này này là dưới 5%. Với trường hợp p = 1, đường cong khoảng giá trị thích hợp nhất của tham số p là 0,5 ≤ p sai số trùng với đường cong sai số của phương pháp ≤1, trong đó với p = 1 cho giá trị bán kính phổ luôn AAM, nói cách khác, độ chính của phương pháp bằng 1, chứng tỏ hệ số cản nhớt số (xem mục 3.3) sẽ PFM1 tương đương với độ chính xác của phương bằng 0. Qua hình 1 ta thấy bán kính phổ cũng bằng 1 pháp AAM. với các giá trị p khác khi Δt/T nhỏ, nhưng nó sẽ giảm xuống (tương ứng với hệ số cản nhớt số tăng lên) khi Δt/T lớn hơn khoảng 0,1. Giá trị p = 0,33 cho giá trị bán kính phổ không phù hợp nên không được xét đến trong phương pháp này. Hình 2. Biểu đồ quan hệ giữa sai số tương đối của chu kỳ và Δt/T với các giá trị p khác nhau Tạp chí KHCN Xây dựng - số 4/2015 5
4. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 3.3 Hệ số cản nhớt số thái ban đầu d0 = 1,0 và v0 = 0. Chọn bước thời gian Δt = 10T = 62,8 s. Kết quả tính toán với 100 bước đầu Như nói ở phần đặt vấn đề, phần lớn các bài toán tiên được thể hiện trong hình 4. phân tích phi tuyến trong xây dựng thuộc dạng thứ nhất, chỉ cần quan tâm đến những dạng dao động bậc Hình 4 cho thấy đường cong biểu diễn kết quả thấp mà không quan tâm đến ảnh hưởng của những tính theo phương pháp PFM1 và AAM hoàn toàn dạng dao động bậc cao, do ảnh hưởng của dạng dao trùng lên nhau, thêm vào đó chuyển vị cũng không bị động bậc cao đến tổng thể kết cấu là không lớn, hơn suy giảm. Trong khi đó, giá trị chuyển vị tính theo nữa, kết quả tính cho những dạng dao động này phương pháp PFM2 bị tắt rất nhanh chỉ sau 10 bước thường kém chính xác nên để đơn giản người ta sẽ đầu tiên, đó là do phương pháp PFM2 có hệ số cản loại bỏ nó đi. Hệ số cản nhớt số (numerical damping nhớt số rất lớn. Với PFM2, giá trị chuyển vị bị vượt ratio) của mỗi phương pháp tính là đại lượng đặc quá không đáng kể trong vài bước đầu tiên, sau đó tắt trưng cho khả năng loại bỏ sự ảnh hưởng của các rất nhanh, trong khi giá trị v0/ω0 bị vượt quá dạng dao động bậc cao mà không làm ảnh hưởng (overshoot) đáng kể trong những bước tính toán đầu đến độ chính xác trong kết quả tính toán của các tiên, tuy nhiên sau đó tắt rất nhanh nên xét về lâu dài dạng dao động bậc thấp. thì kết quả tính vẫn không bị ảnh hưởng. Điều này phù hợp với các kết quả khảo sát như đã trình bày ở mục 3.3. Hình 3. Biểu đồ quan hệ giữa hệ số cản nhớt số của chu kỳ và Δt/T với các giá trị p khác nhau Hình 4. So sánh ảnh hưởng của dao động bậc cao Trong hình 3, đường cong biểu diễn mối quan hệ giữa hệ số cản nhớt số với đại lượng Δt/T được thể 4. Áp dụng với hệ nhiều bậc tự do hiện với các trường hợp tương ứng với p = 1,0; 0,75, Phần lớn các bài toán động lực học công trình là 0,5 cũng như với phương pháp AAM để tham khảo. các bài toán hệ có nhiều bậc tự do, do đó, phần này Biểu đồ cho thấy, với trường hợp p = 0,75 và p = 0,5 sẽ đưa ra cách áp dụng phương pháp đề xuất ở đây các dạng dao động bậc cao (Δt/T) sẽ dễ dàng bị loại để giải các bài toán dạng này. bỏ do có hệ số cản nhớt số lớn, trong khi với các dạng dao động bậc thấp sẽ không bị ảnh hưởng. Với Với hệ nhiều bậc tự do, công thức (2) sẽ được viết như sau: trường hợp p = 1, phương pháp đề xuất ở đây cũng giống phương pháp AAM đều không có hệ số cản 2p p 1 2 p p 1 Ma C v R R f f nhớt số, được dùng để tính toán với các bài toán có i 1 0 i 1p 1 i 1 p 1 i p 1 i 1 p 1 i xét đến đến ảnh hưởng của cả các dạng dao động 2 (6) d B d B d B t v B t a bậc thấp và bậc cao. i 1 0 i 1 1 i 2 i 3 i 3p 1 p 3 3.4 Ảnh hưởng của dao động bậc cao v v t a t a i 1 i2 p 1 i 2 p 1 i 1 Để làm rõ hơn đặc điểm về hệ số cản nhớt số của trong đó: phương pháp đề xuất, ta khảo sát thêm tính chất nữa. Tính chất này thường được biết đến trong thuật ngữ M, C0 là các ma trận khối lượng, ma trận độ cản tiếng Anh với tên gọi là overshooting. Xét một hệ một nhớt vật lý của kết cấu; bậc tự do có khối lượng m = 1kg và độ cứng k = a, v, d, f tương ứng là các vec-tơ gia tốc, vận tốc, 1N/m, chu kỳ dao động tự do của hệ sẽ là: chuyển vị và ngoại lực tác dụng; T 2 m k 6,28 s . Cho hệ dao động tự do từ trạng 6 Tạp chí KHCN Xây dựng - số 4/2015
5. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Các hậu tố (0), (i-1), (i), (i+1) chỉ thứ tự các bước 5. Một số ví dụ tính toán tính toán; Để làm rõ hơn các đặc điểm của phương pháp Ri, Ri+1 là vec-tơ nội lực của hệ kết cấu ở bước này, một số ví dụ sẽ được trình bày ở đây. Các ví dụ tính toán thứ (i) và (i+1). này sẽ so sánh các đặc điểm của phương pháp đề xuất PFM với hai phương pháp phân tích phi tuyến Ma trận B0 đến B3 được tính như sau: được biết đến rất rộng rãi là NEM [13], đại diện cho 3 1 p 1 2 2 phương pháp ngoại hiển thức và AAM [13], tiêu biểu BDK t 08 p 1 0 cho phương pháp nội ẩn thức. 3 Nhìn chung, với tất cả các phương pháp phân tích 1 p 1 2 2 BIDK t động phi tuyến đều cần sử dụng máy tính do khối 18 p 1 0 lượng tính toán lớn. Với những bài toán đơn giản, (7) 1 p 3 BDMC t người dùng có thể tự lập trình bằng phần mềm như 2 0 2 p 1 Excel hoặc viết các chương trình dựa trên ngôn ngữ 2  lập trình Fortran, Matlab, C++ Với các bài toán 1 2 p 3 BDMC1 1 t 2 4  phức tạp hơn thì nên sử dụng những phần mềm 3p 1 p 1 0  chuyên dụng như Sap, Etabs hay OpenSees. Các ví 3 dụ tính toán ở đây được tác giả tự lập trình bằng p 3 p 2 2 DMCK t t Matlab. 2 p 1 0 4 p 1 0 5.1 Ví dụ 1: Hệ một bậc tự do đàn dẻo tuyệt đối với I là ma trận đơn vị đường chéo chính (ma trận vuông với các giá trị ở đường chéo chính bằng 1), K0 là ma trận độ cứng của hệ kết cấu ở thời điểm ban di đầu (initial stiffness). m =104 kg Cần nói thêm rằng, các ma trận từ B0 đến B3 được xác định chỉ dựa vào các đặc điểm ban đầu của kết cấu (điều kiện trước khi biến dạng) là M, C , K và 0 0 k =106 N/m giá trị bước thời gian Δt, do đó chỉ cần tính một lần trong suốt quá trình tính toán, điều này làm tiết kiệm được nhiều công sức. Dòng thứ hai của công thức (6) cho thấy lời giải của phương trình vi phân này là lời ag giải gồm hai bước, việc tính chuyển vị ở bước (i+1) được tính từ các giá trị ở bước (i) và bước (i-1) trước Hình 5. Mô hình thí nghiệm trên bàn rung với hệ một bậc tự do đó, do vậy cần có một bước đệm khi tính toán với Giả thiết có một hệ một bậc tự do được thí bước đầu tiên. Để ý rằng, tham số B0di-1 sẽ triệt tiêu nếu p = 1, do vậy, ta gán giá trị p = 1 cho bước đầu nghiệm trên bàn rung như hình 5. Tải trọng tập trung 4 tiên. Với các bước tiếp theo, giá trị hệ số p được lựa ở đầu cột bằng m = 10 kg, cột giả thiết như một thanh 6 chọn tùy theo yêu cầu tính toán. đàn dẻo tuyệt đối, độ cứng k = 10 N/m, do đó chu kỳ dao động tự nhiên ban đầu của hệ ω0 = 10 rad/s. Quy trình tính toán với hệ nhiều bậc tự do được Cường độ chịu kéo và chịu nén của vật liệu thanh giả thực hiện như sau: thiết bằng Rt = Rc = 50kN. Bỏ qua hệ số cản nhớt vật Bước 1: Tính giá trị vec-tơ chuyển vị di+1 từ dòng lý của hệ. Gia tốc nền tác dụng lên hệ được điều thứ hai của công thức (6); khiển thông qua kích thủy lực của bàn rung, được chọn theo gia tốc nền ghi nhận được từ trận động đất Bước 2: Thế giá trị vec-tơ chuyển vị d vừa tìm i+1 Chi-Chi xảy ra ở Đài Loan vào năm 1999 (tên phổ ghi được và giá trị vec-tơ vận tốc v ở dòng thứ ba vào i+1 gia tốc này theo ký hiệu quốc tế thường dùng là CHY dòng thứ nhất cùng của công thức (6), giải và tìm vec- 028), đỉnh gia tốc nền được lấy bằng 0,5g. Dùng tơ gia tốc a ; i+1 phương pháp PFM1 để dự đoán chuyển vị của hệ Bước 3: Giá trị vec-tơ vận tốc được tính bằng đồng thời so sánh kết quả với hai phương pháp NEM dòng thứ ba của công thức (6) sau khi đã có giá trị và AAM. vec-tơ vận tốc vi+1. Tạp chí KHCN Xây dựng - số 4/2015 7
6. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 6 thể hiện kết quả tính toán chuyển vị theo làm chuẩn để so sánh với các phương pháp khác, thời gian của hệ và biểu đồ quan hệ giữa chuyển vị và trong khi bước thời gian tính toán theo AAM và PFM1 nội lực thanh. Giá trị bước thời gian tính toán theo là Δt =0,02 s. phương pháp NEM là Δt =0,005 s, đủ nhỏ để coi như kết quả tính từ phương pháp này là chính xác, Hình 6. Chuyển vị của hệ dưới tác dụng của tải động đất và quan hệ chuyển vị - nội lực của thanh Kết quả tính toán thể hiện trong hình 6a cho thấy giá trị chuyển vị thu được từ phương pháp PFM1 rất sát với kết quả thu được từ phương pháp AAM, chỉ sai lệch rất nhỏ so với kết quả thu được từ phương pháp NEM. Kết quả từ hình 6b cho thấy cột đã hoàn toàn ứng xử phi tuyến. Như vậy, có thể thấy rằng phương pháp đề xuất ở đây hoàn toàn có thể sử dụng để phân tích phi tuyến với độ chính xác cao. 5.2 Ví dụ 2: Dao động tự do hệ nhiều bậc tự do Hình 7. Mô hình 5 tầng – 5 bậc tự do Hình 7 thể hiện mô hình một hệ kết cấu 5 tầng, k0-i là độ cứng ban đầu (trước khi biến dạng) của mỗi tầng được mô hình hóa như một tấm cứng tuyệt tầng thứ i; đối, toàn bộ hệ kết cấu 5 tầng được coi như có 5 bậc |ui – ui-1| là chuyển vị lệch tầng của tầng thứ i so tự do. Các giá trị tải trọng và độ cứng mỗi tầng được với tầng thứ (i-1) bên dưới; thể hiện như trong hình vẽ. Chu kỳ dao động tự nhiên hệ, vec-tơ dạng dao động thứ 1, 4, 5 cũng được thể Hệ số pi đặc trưng cho tính phi tuyến của mỗi hiện trong hình vẽ. Độ cứng của mỗi tầng được thể tầng, pi = 0 tương ứng với thanh đàn hồi, pi 0 1/5 tương ứng với hệ cứng hóa sau biến dạng k k 1 p u u , i 1~ 5 (8) j i0 i i i i 1 (hardening). Hình 8 minh họa mối quan hệ này với đường nét liền thể hiện cho mối quan hệ tuyến tính trong đó: giữa chuyển vị lệch tầng và nội lực trong cột, đường kj-i là độ cứng tức thời của tầng thứ i tại thời điểm nét đứt thể hiện mối quan hệ phi tuyến, ở đây là quan cuối cùng của bước thứ j; hệ mềm hóa. 8 Tạp chí KHCN Xây dựng - số 4/2015
7. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 8. Quan hệ giữa chuyển vị lệch tầng và nội lực trong cột Hình 9. Dao động tầng 5 của hệ kết cấu dưới tác dụng Trong ví dụ này, mô hình công trình chịu tác dụng của tải trọng động đất của lực động đất, chuyển động với gia tốc nền CHY 5.3 Ví dụ 3: So sánh về hiệu quả tính toán 028, đỉnh gia tốc bằng 0,5g như đã dùng ở ví dụ 1. Độ cứng của mỗi tầng được tính như công thức (8) với giá trị p1 = -1,0, p2~p5 = -0,5. Để đánh giá ảnh hưởng của các dạng dao động bậc cao, chuyển vị ban đầu của hệ được gán bằng d(0) = Ф5/10 (m), trong đó giá trị của vec-tơ Ф5 thể hiện trong hình 7. Chuyển vị của tầng thứ 5 được tính bằng phương pháp AAM, PFM1 và PFM2 được thể hiện trong hình 9. Ngoài ra, kết quả của một trường hợp tính bằng phương pháp AAM trong đó chuyển vị ban đầu của hệ kết cấu được Hình 10. Mô hình hệ nhiều bậc tự do tính từ trạng thái nghỉ d(0) = 0 được in trong hình 9a Xem xét mô hình một hệ phi tuyến nhiều bậc tự để so sánh. Tất cả các trường hợp này đều tính với do như trong hình 10. Hệ gồm n vật nặng, mỗi vật có bước thời gian Δt = 0,01s. khối lượng 1kg, nối với nhau bằng lò xo có độ cứng ki, Kết quả từ hình 9 cho thấy chuyển vị tại tầng 5 trong đó giá trị ki thay đổi phụ thuộc vào biến dạng của hệ kết cấu bị ảnh hưởng rất nhiều của các dạng của lò xo, có công thức tính được ghi trong hình 10. dao động bậc cao khi tính bằng phương pháp AAM và Toàn bộ hệ chịu tác dụng của chuyển vị nền, dao PFM1. Các giá trị dao động trở nên rất “nhiễu” (không động theo dạng hình sin như trong hình vẽ. Xem xét chính xác) và gần như không thể sử dụng được để hai trường hợp, hệ có 500 bậc tự do (n = 500) và hệ phân tích kết cấu. So sánh hình 9a và 9b cho thấy, độ có 1000 bậc tự do (n = 1000). Về nguyên tắc thì mô “nhiễu” theo tính toán bằng phương pháp PFM1 thậm hình hệ lò xo này cũng giống mô hình hệ 5 tầng 5 bậc chí còn lớn hơn độ “nhiễu” theo tính toán bằng AAM. tự do trong ví dụ 2 (khác số bậc tự do). Chu kỳ dao Tuy nhiên, độ “nhiễu” này đã hoàn toàn bị loại bỏ động tự nhiên nhỏ nhất được tính toán dựa vào độ trong kết quả của phương pháp PFM2. Kết quả cứng ban đầu trước khi biến dạng với hệ 500 bậc tự chuyển vị thu được từ hình 9c cho thấy đường cong do có kết quả bằng ω0 = 31,38 rad/s, trong đó với hệ biểu diễn chuyển vị rất trơn và trùng với kết quả thu 1000 bậc tự do bằng ω0 = 15,70 rad/s, trong khi chu được từ phương pháp AAM trong trường hợp hệ dao kỳ dao động tự nhiên lớn nhất với cả hai trường hợp động từ trạng thái nghỉ. đều bằng ω0 = 20000 rad/s. Hệ được tính toán bằng 3 phương pháp là NEM, AAM và PFM2. Qua ví dụ này cho thấy, chỉ bằng cách điều chỉnh hệ số p của phương pháp đề xuất ở đây, ta hoàn toàn Bước thời gian tính toán Δt đối với phương pháp có thể điều chỉnh được hệ số cản nhớt số để loại bỏ NEM được lựa chọn căn cứ theo điều kiện ổn định ảnh hưởng của các dạng dao động bậc cao mà không Ω0 = ω0(Δt)≤2, theo đó với chu kỳ dao động lớn nhất làm ảnh hưởng đến độ chính xác của kết quả tính ω0 = 20000 rad/s thì bước thời gian được lựa chọn toán. Với giá trị p = 1, phương pháp này không có hệ Δt ≤ 2/20000 = 0,0001s. Giá trị tính toán của bước số cản nhớt số, giá trị p giảm cho hệ số cản nhớt số thời gian này cũng đủ nhỏ để có thể coi kết quả tính càng tăng, và hệ số cản nhớt số lớn nhất khi p = 0,5. toán thu được là chính xác, làm cơ sở để so sánh với Tạp chí KHCN Xây dựng - số 4/2015 9
8. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG các phương pháp khác. Trong khi đó, do hai phương trận hai lần, NEM chỉ cần giải một lần). Tuy nhiên, do pháp AAM và PFM2 đều là những phương pháp NEM bị giới hạn về điều kiện ổn định nên số bước không có điều kiện ổn định nên bước thời gian Δt tính toán bị tăng lên rất nhiều, ở đây số bước tính được lựa chọn cho hai phương pháp AAM và PFM2 toán của phương pháp NEM lớn hơn phương pháp căn cứ vào yêu cầu chính xác của kết quả tính toán PFM2 tới 50 lần. Do vậy, xét về tổng thể thì thời gian với Δt càng nhỏ thì độ chính xác càng cao, nhưng thời gian tính toán bằng NEM vẫn nhiều hơn PFM2. ngược lại khối lượng tính toán càng lớn. Ở đây lựa Mặt khác, khi so sánh với phương pháp AAM, tuy chọn Δt = 0,005s giây cho cả hai phương pháp AAM PFM2 và AAM có số bước tính toán bằng nhau và PFM2. nhưng trong mỗi bước AAM phải tính lặp rất nhiều lần để tìm ra đáp số chính xác (ví dụ: sử dụng phương Hình 11 thể hiện kết quả kết quả tính toán chuyển pháp Newton Raphson thường mất trung bình khoảng vị của lò xo thứ 500 và thứ 1000 tương ứng với hai hệ 10-25 lần lặp trong một bước). Việc tính lặp này rất đã nêu ở trên. Biểu đồ cho thấy đường cong biểu diễn mất thời gian, công sức và kết quả không phải lúc nào chuyển vị gần như trùng nhau, cho thấy kết quả tính cũng hội tụ hoặc hội tụ đến kết quả không mong toán từ phương pháp AAM và PFM2 là chính xác. muốn. Do vậy, trong một bước tính toán thì tính bằng Tiếp theo so sánh về thời gian tính toán, cả ba AAM mất thời gian và công sức hơn PFM2 khá nhiều, phương pháp này đều được tính toán bằng chương tùy thuộc vào tốc độ hội tụ và độ chính xác yêu cầu. trình viết bằng Matlab, chạy trên máy tính cá nhân Trong phân tích ứng xử phi tuyến cho công trình dùng chip Intel®CoreTMi5 CPU M460 @2.53GHz, RAM xây dựng, việc giảm thời gian tính toán có ý nghĩa rất 4.00 GB. Thời gian tính toán (tính bằng giây) cho từng lớn. Để giảm thời gian tính toán có thể sử dụng hệ 500 và 1000 bậc tự do tương ứng với mỗi phương những hệ máy tính mạnh hoặc chạy song song nhiều pháp được trình bày trong bảng 1. máy tính, tuy nhiên cách này thường gây tốn kém chi Kết quả so sánh từ bảng 1 cho thấy, thời gian tính phí, do đó sử dụng một phương pháp có hiệu quả tính toán của phương pháp PFM2 chỉ bằng khoảng 5% so toán cao như phương pháp đề xuất ở đây cũng là một với phương pháp NEM và bằng dưới 2% so với cách tốt. Trong ví dụ này chỉ tính toán với hệ 1000 phương pháp AAM. bậc tự do, với hệ có số bậc tự do càng lớn, chẳng Điều này được lý giải là do trong một bước tính hạn như mô hình một công trình xây dựng có thể lên toán thì tính bằng phương pháp PFM2 mất nhiều đến hàng triệu bậc tự do, thì việc tiết kiệm này sẽ công sức hơn (PFM2 phải giải hệ phương trình ma càng có ý nghĩa. Hình 11. Kết quả tính toán chuyển vị 10 Tạp chí KHCN Xây dựng - số 4/2015
9. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Bảng 1. Bảng so sánh về hiệu quả tính toán N-DOF NEM (1) AAM(2) PFM2 (3) (3)/(1) (3)/(2) 500 640,69 s 1733,47 s 32,94 s 0,051 0,019 1000 2955,91 s 14224,70 s 140,81 s 0,048 0,0099 4. Kết luận [7]. Chang, S.Y. (2010). “A new family of explicit method Trong bài báo này, một họ phương pháp phân for linear structural dynamics”, Computers & tích phi tuyến động theo lịch sử thời gian được đề Structures, Vol. 88, No.11-12:755-772. xuất, trong đó các đặc trưng số học của phương pháp [8]. Chang, S.Y., N.C. Tran. (2014). “A two-step được điều chỉnh chỉ thông qua một hệ số p duy nhất. unconditionally stable explicit method with Hệ số này được lựa chọn trong khoảng 0,5 ≤ p ≤ 1, controllable numerical dissipation for structural với p = 1 cho phép tính toán không có hệ số cản nhớt dynamics”, Advances in Civil Engineering and số, p giảm cho hệ số cản nhớt số tăng. Phương pháp Building Materials IV. 379-383. thuộc họ ngoại hiển thức không có điều kiện ổn định, [9]. Dobbs, M.W. (1974). “Comments on ‘stability and do vậy tiết kiệm được rất nhiều công sức tính toán so accuracy analysis of direct integration methods’ với hai phương pháp truyền thống được so sánh là by Bathe & Wilson”, Earthquake Engineering and NEM và AAM, trong khi độ chính xác trong kết quả Structural Dynamics, Vol. 2: 295-299. thu được là tương đương. Thông qua các ví dụ đã chỉ [10]. Hilber, H.M., Hughes, T.J.R. and Taylor, R.L. ra rằng phương pháp này là phù hợp để giải các bài (1977). “Improved numerical dissipation for time toán phi tuyến trong xây dựng. integration algorithms in structural dynamics”, Lời cảm ơn Earthquake Engineering and Structural Nhóm tác giả xin cảm ơn Hội đồng Khoa học Quốc Dynamics, Vol. 5: 283-292. gia Đài Loan đã tài trợ kinh phí cho nghiên cứu này [11]. Hughes, T.J.R. (1987). “The Finite element method”, thông qua hợp đồng số NSC-99-2221-E-027-029. Bài Prentice-Hall, Inc., Englewood Cliffs, N.J., U.S.A. báo được gửi đăng với sự đồng ý của ©2015 Taylor & [12]. Krieg, R.D. (1973). “Unconditional stability in Francis Group, London, UK. numerical time integration methods”, Journal of Applied Mechanics, Vol. 40: 417-421. TÀI LIỆU THAM KHẢO [13]. Newmark, N.M. (1959). “A method of computation [1]. Bathe, K.J. and Wilson, E.L. (1973). “Stability and for structural dynamics”, Journal of Engineering accuracy analysis of direct integration methods”, Mechanics Division, ASCE, Vol. 85: 67-94. Earthquake Engineering and Structural Dynamics. [14]. Zhou X, Tamma KK. (2004). “Design, analysis Vol. 1: 283-291. and synthesis of generalized single step single [2]. Belytschko, T. and Schoeberle, D.F. (1975). “On solve and optimal algorithms for structural the unconditional stability of an implicit algorithm dynamics”, International Journal for Numerical for nonlinear structural dynamics”, Journal of Methods in Engineering; Vol. 59: 597-668. Applied Mechanics, Vol. 17: 865-869. [15]. Zhou X, Tamma KK. (2006). “Algorithms by design [3]. Belytschko, T. and Hughes, T.J.R. (1983). with illustrations to solid and structural mechanics/ “Computational methods for transient analysis”, dynamics”, International Journal for Numerical Elsevier Science Publishers B.V., North-Holland. Methods in Engineering; Vol. 66: 1841-1870. [4]. Chang, S.Y. (2002). “Explicit pseudodynamic algorithm [16]. Zienkiewicz, O.C. (1977). “The Finite Element with unconditional stability”. Journal of Engineering Method”, McGraw-Hill Book Co (UK) Ltd. Third Mechanics, ASCE, Vol. 128, No.9: 935-947. edition. [5]. Chang, S.Y. (2007). “Improved explicit method for [17]. W.L. Wood, M. Bossak, and O.C. Zienkiewicz. structural dynamics”, Journal of Engineering “An alpha modification of Newmark’s method”, Mechanics, ASCE, Vol. 133, No. 7: 748-760. International Journal for Numerical Methods in [6]. Chang, S.Y. (2009). “An explicit method with Engineering, Vol. 15: 1562-1566, 1981. improved stability property”, International Journal Ngày nhận bài: 03/8/2015. for Numerical Method in Engineering, Vol. 77, No. Ngày nhận bài sửa lần cuối: 06/11/2015. 8:1100-1120. Tạp chí KHCN Xây dựng - số 4/2015 11
10. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG ỨNG DỤNG PHẦN MỀM MÃ NGUỒN MỞ OPENSEES TRONG LẬP TRÌNH MÔ PHỎNG CẦU CHỊU ĐỘNG ĐẤT KS. TRẦN TIẾN ĐẠT, KS. NGUYỄN ĐỨC PHÚC, TS. TRẦN ANH BÌNH Trường Đại học Xây dựng Tóm tắt: OpenSees (Open System of Earthquake Navigator, GiD, TclBuilder hoặc toolbox trong Matlab Engineering Simulation) là một phần mềm mã nguồn như OpenSees Pre- and Post- Processing [2] để xuất mở với thư viện mã code được viết chủ yếu bằng và nhập kết quả. Tuy nhiên, việc sử dụng các phần C++, một ngôn ngữ lập trình hướng đối tượng. Nó cho mềm hỗ trợ này vẫn chưa đầy đủ, hoàn chỉnh như phép người dùng tạo ra các mô hình phần tử hữu hạn các phần mềm thương mại hiện nay. Bởi vì để mô phỏng phản ứng của hệ kết cấu và đất nền OpenSees hỗ trợ nhiều phương pháp giải khác nhau, dưới tác dụng của động đất. Bài báo này sẽ giới thiệu do đó đối với người dùng thông thường sẽ rất khó về OpenSees từ đó xây dựng thuật toán ứng dụng khăn trong việc lựa chọn, thậm chí dẫn đến tính toán OpenSees vào trong lập trình mô phỏng một ví dụ kết đưa ra kết quả sai nếu người dùng không hiểu hoặc cấu cầu chịu tác dụng của động đất. dùng sai các phương pháp tính. Về cơ bản, ở thời điểm hiện tại OpenSees dành cho người dùng trong Từ khóa: OpenSees, động đất, mô phỏng. công tác nghiên cứu và phát triển chứ chưa thực sự 1. Mở đầu hướng đến những người dùng phổ thông. OpenSees (Open System of Earthquake Sử dụng OpenSees trong phân tích kết cấu xây Engineering Simulation) là một phần mềm mã nguồn dựng nói chung và đặc biệt trong kết cấu cầu nói mở được phát triển bởi nhóm nghiên cứu ở PEER riêng được nhiều nhóm nghiên cứu quan tâm, Jinchi (Pacific Earthquake Engineering Research) của Lu, Kevin Mackie, and Ahmed Elgamal [3] đã xây trường đại học University of California, Berkeley. Thư dựng chương trình phân tích Pushover cho kết cấu viện mã code của OpenSees được viết chủ yếu bằng nhịp giản đơn; Pallavi Gavali, Mahesh S. Shah, Gouri C++ là ngôn ngữ lập trình hướng đối tượng. Kadam và Kranti Meher [4] đã xây dựng chương trình OpenSees là một phần mềm mã nguồn mở miễn phân tích 3D kết cấu cầu và đất nền chịu động đất. phí. Phần mềm có khả năng tính toán mạnh với một OpenSees hiện tại và trong tương lai gần vẫn là thư viện hỗ trợ nhiều loại vật liệu, người dùng có thể một công cụ mạnh trong phát triển và mô phỏng kết tự phát triển các loại vật liệu khác tùy thuộc vào nhu cấu chịu động đất. Bài báo này sẽ giới thiệu về cầu của người sử dụng. OpenSees hỗ trợ nhiều OpenSees từ đó xây dựng sơ đồ thuật toán ứng dụng phương pháp tính khác nhau, do đó tùy vào mục đích vào trong lập trình mô phỏng một ví dụ kết cấu cầu cụ và đặc điểm của bài toán kết cấu mà người dùng có thể chịu tác dụng của động đất. thể linh hoạt lựa chọn phương pháp thích hợp để 2. Tổng quan về OpenSees phân tích. Mã nguồn của OpenSees thường xuyên được cập nhật và đóng góp của nhiều cơ quan 2.1 Những tính năng của OpenSees nghiên cứu lớn trên thế giới, các chuyên gia từ các Như đã đề cập trong phần mở đầu, OpenSees là trường đại học, viện nghiên cứu, đây là phần mềm có một phần mềm mã nguồn mở phục vụ cho phân tích lượng người dùng lớn do đó việc trao đổi, thảo luận là tính toán kết cấu theo phương pháp phần tử hữu hạn. tương đối dễ dàng [1]. Tính năng mạnh nhất của OpenSees là mô phỏng Tuy nhiên, OpenSees có nhược điểm là một phần phản ứng của kết cấu chịu tác động của động đất, các mềm không có giao diện đồ họa, vì vậy gây khó khăn bài toán phân tích phi tuyến và phân tích kết cấu cho những người dùng thông thường, hiện nay có tương tác với đất nền. Sơ đồ hình 1 thể hiện các ưu một số phần mềm đã được viết như OpenSees điểm mà OpenSees mang lại: 12 Tạp chí KHCN Xây dựng - số 4/2015
11. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 1. Các ưu điểm của OpenSees Mô hình kết cấu trong phần mềm mã nguồn mở OpenSees người sử dụng có thể tùy chọn hoàn toàn thông qua API (các hàm, thủ tục được lập sẵn) về mô hình, hiển thị, phương pháp giải, kết quả đầu ra. So với phần mềm mã nguồn đóng truyền thống ví dụ như Midas/Civil chưa hỗ trợ bộ API [5], Sap2000 chỉ cho phép đưa dữ liệu đầu vào, xuất kết quả qua các hàm API [6]. Khả năng tùy chỉnh của các phần mềm mã nguồn đóng phụ thuộc vào bộ API do phần mềm cung cấp do đó rất hạn chế (hình 2). Hình 2. So sánh phần mềm mã nguồn đóng và OpenSees [7] Ngoài ra, bộ mã nguồn mở được cung cấp miễn phí [8], người sử dụng có thể xây dựng thêm cơ sở dữ liệu cho bài toán của mình như tạo phần tử, vật liệu, phương pháp giải, lặp, thuật giải mới và sau đó xây dựng thành trình dịch của riêng mình. Hình 3. Lớp trừu tượng chính (main abstractions) 2.2 Mô hình trong OpenSees trong OpenSees [7] OpenSees bao gồm một tập hợp các mô-đun để Các loại đối tượng chính này thực hiện các nhiệm thực hiện công việc tạo ra các mô hình phần tử hữu vụ khác nhau, cụ thể như: hạn. Trong một mô hình phân tích phần tử hữu hạn ModelBuilder: Là một đối tượng trong chương của kết cấu trong OpenSees, các lệnh được sử dụng trình chịu trách nhiệm tạo ra mô hình kết cấu như nút, với mục đích để tạo ra 4 loại đối tượng chính (hình 3). phần tử, vật liệu, định nghĩa các loại tải trọng và định nghĩa điều kiện biên. Tạp chí KHCN Xây dựng - số 4/2015 13
12. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Domain: Có trách nhiệm lưu trữ các đối tượng được tạo ra từ ModelBuilder và cho phép các lớp chính Analysis và Recoder truy nhập tới các đối tượng này (hình 4). Hình 4. Domain lưu trữ các đối tượng tạo ra từ ModelBuilder [7] Recorder: Theo dõi các tham số mà người dùng đã định nghĩa trong suốt quá trình phân tích, ví dụ: theo dõi chuyển vị tại một nút theo lịch sử thời gian. Lớp trừu tượng này cũng điều khiển ghi file và xuất kết quả. Analysis: Chịu trách nhiệm thực hiện phân tích. Trong OpenSees lớp trừu tượng chính Analysis bao gồm các đối tượng điều khiển kiểu phân tích cho mô hình (hình 5). Hình 5. Các đối tượng trong lớp Analysis [7] Để thực hiện xây dựng và phân tích một mô hình phần tử hữu hạn trong OpenSees người dùng sử dụng các hàm đã được lập sẵn trong OpenSees và sử dụng ngôn ngữ TCL (một ngôn ngữ thông dịch mạnh và dễ sử dụng) để tạo ra mô hình hình học, tải trọng (file nguồn). 2.3 Thuật toán ứng dụng mã nguồn mở OpenSees mô phỏng cầu chịu động đất Để minh họa cho phương pháp mô hình kết cấu trong OpenSees. Bài báo trình bày sơ đồ thuật toán (hình 6) được xây dựng để áp dụng trực tiếp vào ví dụ ở mục 3 như sau: Hình 6. Sơ đồ thuật toán 14 Tạp chí KHCN Xây dựng - số 4/2015
13. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 3. Ví dụ số và phân tích kết quả 3.1 Số liệu hình học và gia tốc nền 95000 l=30000 l=35000 l=30000 10000 Hình 7. Mặt đứng kết cấu cầu Trong ví dụ sử dụng một kết cấu dầm bản rỗng 3 nhịp (30m + 35m + 30m) trụ dạng cột tiết diện tròn vật liệu bê tông cốt thép, trụ và dầm được liên kết cứng với nhau, hai gối di động được đặt ở hai mố đầu cầu. Sơ đồ kết cấu đối xứng qua tim cầu theo cả phương dọc và ngang (hình 7). 24-D32 D13@75 Hình 8. Mô hình kết cấu đưa vào OpenSees Một mô hình 3D của kết cấu được lập trình trong OpenSees như sau (hình 8): kết cấu dầm chủ được mô hình bằng phần tử dầm – cột đàn hồi (elasticBeamColumn). Để mô hình liên kết cứng giữa trụ và dầm sử dụng một liên kết tuyệt đối cứng trong OpenSees (rigidLink). Trụ cầu được mô hình xét đến phi tuyến vật liệu sử dụng phần tử dầm – cột phi tuyến (nonlinearBeamColumn) với mô hình tiết diện bê tông cốt thép (Fiber Section). Khối lượng kết cấu được gán vào các nút trên các phần tử dầm và trụ. Trọng lượng bản thân kết cấu và hiệu ứng P – Delta được đưa vào cùng phân tích động đất. Mô hình liên kết ngàm (SP-Constraint-fix) ở đáy đài cọc. Dữ liệu trận động đất được sử dụng cho ví dụ được lấy từ PEER Ground Motion Database [9]. Cụ thể ở đây là 2 trận động đất EI Centro 12 và EI Centro 01 của trạm đo Imperial Valley 10/15/79 (bảng 1). Bảng 1. Chi tiết trận động đất sử dụng trong ví dụ Trận động Cấp động Bước thời Hướng Trạm đo PGA (g) đất đất (Mw) gian (s) Phương dọc Cầu (Trục X Imperial valley 10/15/79 El Centro12 6.5 0.005 0.116 hệ tọa độ tổng thể) EL CENTRO ARRAY #12 Theo phương ngang Cầu Imperial valley 10/15/79 El Centro01 6.5 0.005 0.139 (Trục Y hệ tọa độ tổng thể) EL CENTRO ARRAY #1 Tạp chí KHCN Xây dựng - số 4/2015 15
14. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Hình 9. Biểu đồ gia tốc nền Bảng 2. Bảng tổng hợp về mô hình kết cấu trong OpenSees Phần kết cấu Mô hình phân tích Chi tiết Diện tích tiết diện: 6.9 m2 Dầm chủ Mô hình 3D phần tử dầm – cột đàn hồi Mô men quán tính: 1.31 m4 3Mô hình 3D phần tử dầm – cột phi tuyến, phân Diện tích tiết diện: 1.7 m2 Trụ cầu tích ảnh hưởng hiệu ứng P – Delta. Mô men quán tính: 0.25 m4 3.2 Thông số vật liệu Mô hình bê tông: Mô hình bê tông được sử dụng trong ví dụ này được giới thiệu bởi Kent và Park và sau đó được mở rộng bởi Scott [10]. Sử dụng hàm Concrete01 đã được thiết lập sẵn trong OpenSees [7]. 2 ' εc ε c σc = Kf c 2 - εc ε o (1) ε ε co co ' σc =Kf c 1-Z(ε c -ε co ) εco ε c ε u (2) ' σc = 0.2Kf c εcu ε c (3) ρs fyh K = 1+ ' εco = 0.002K (4) fc 0.5 Z = ' ' 3 + 0.0284fc h (5) Hình 10. Mô hình bê tông [10] + 0.75ρ - 0.002K ' s 14.21fc -1000 sh Ở đây  co là biến dạng trong bê tông khi ứng suất đạt lớn nhất; K là hệ số tăng cường độ do bê tông bị kiềm ' Z f 2 f chế; là độ dốc đường biến dạng; c là cường độ chịu nén của bê tông kg/cm ; yh là cường độ chảy của cốt thép đai; ρs là tỉ số thể tích của thép ngang kiềm chế trên thể tích của lõi kiềm chế; sh là khoảng cách từ tim đến ' tim đai; h là chiều rộng của lõi kiềm chế tính đến mép ngoài thép đai. Bảng 3. Các thông số bê tông sử dụng trong mô hình H.lượng E γ Cốt thép Bê tông c f’ (Mpa) ε ε K Z thép dọc (Mpa) c (kg/m3) 0 u đai ρ (%) Thường 33994 40.000 24.5 0.0020 0.006 1.00 64.5 1.112 D13@75 Bị kiềm chế 33994 43.446 24.5 0.0217 0.02 1.086 469.15 1.112 D13@75 16 Tạp chí KHCN Xây dựng - số 4/2015
15. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Mô hình cốt thép: Quan hệ ứng suất – biến dạng cản tỉ số cản được giả thiết là ξ = 5% cho 2 tần số của cốt thép được mô hình đơn giản là các đoạn dao động đầu tiên [11]. Để tiết kiệm bộ nhớ máy tính thẳng tuyến tính (hình 11). Lý do sử dụng mô hình khi phân tích sử dụng phương pháp chuyển ma trận xấp xỉ này là tiện lợi cho mô hình tính toán. Mô hình dải của ma trận độ cứng về một ma trận dạng chữ này đã được OpenSees đưa vào hàm Steel01 [7]. Số nhật bằng lệnh system BandGeneral. Số bước lặp liệu cốt thép đưa vào ví dụ như sau: Cường độ chảy lớn nhất 100 bước và độ hội tụ 10-8 được đưa vào để của cốt thép: fy = 420 MPa; mô đun đàn hồi của thép: dừng phân tích khi vượt quá số bước lặp mà kết quả E0 = 200000 MPa; tỷ số biến cứng: b = 0.01; trọng chưa hội tụ, sử dụng hàm test NormUnbalance $tol lượng riêng thép: γ = 7500 kg/m3. $iter trong đó $tol là độ hội tụ và $iter số bước lặp lớn nhất. Dữ liệu động đất sẽ được phóng lên 3 lần để cho thấy ảnh hưởng lớn của các hiệu ứng phi tuyến trong kết cấu. 3.4 Mô hình trong phần mềm Midas/ Civil Để so sánh ưu nhược điểm của OpenSees với phần mềm thương mại hiện nay trong lĩnh vực phân Hình 11. Mô hình vật liệu thép [7] tích kết cấu, đặc biệt là kết cấu cầu. Bài báo sử dụng phần mềm Midas/Civil để phân tích và so sánh kết 3.3 Cài đặt các thông số phân tích cho mô hình quả. Midas/Civil là phần mềm thương mại của hãng trong OpenSees Midas IT (Hàn Quốc). Đây là phần mềm phân tích kết Để thực hiện phân tích động phi tuyến, phương cấu có độ chính xác cao cùng với các phần mềm như pháp tích phân trực tiếp Newmark với Δt = 0.01s trên Sap2000, ADINA và ANSYS. Midas/Civil hiện nay cơ sở thuật toán gia tốc trung bình được sử dụng, được nhiều kỹ sư và nhà khoa học sử dụng để phân phương pháp này đã được lập sẵn trong hàm tích kết cấu đặc biệt là trong kết cấu cầu. Các thông integrator Newmark $gamma $beta. Thuật giải lặp số đầu vào đưa vào phần mềm Midas/Civil tương tự Newton-Raphson được khai báo trong OpenSees như các thông số được sử dụng lập trình trong bằng hàm algorithm Newton. Để thiết lập ma trận OpenSees. Hình 12. Mô hình kết cấu trong Midas/Civil Tạp chí KHCN Xây dựng - số 4/2015 17
16. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 3.5 Kết quả phân tích a. Tần số dao động riêng Trong bài toán động lực học công trình thì tần số dao động riêng của kết cấu là một giá trị đặc biệt quan trọng, do đó kết quả của phân tích tần số dao động riêng của hai chương trình được đưa ra để so sánh (hình 13). Kết quả phân tích 5 tần số dao động riêng từ OpenSees và Midas/Civil cho giá trị gần giống nhau, cụ thể là ở dạng dao động đầu tiên sai số giữa 2 chương trình là 0.87%, sai số nhỏ nhất là ở dạng dao động thứ 3 với sai số là 0.042%, sai số tăng dần cho các dạng dao động tiếp theo với 1.85% và 2.23% cho dạng dao động thứ 4 và thứ 5. Hình 13. Đồ thị so sánh kết quả tần số dao động riêng b. Chuyển vị theo lịch sử thời gian  Chuyển vị dọc cầu tại vị trí đỉnh trụ kết quả phân tích từ OpenSees khi có xét và không xét đến hiệu ứng P-Delta. Hình 14. So sánh chuyển vị dọc khi xét hiệu ứng P-Delta trong OpenSees Kết quả so sánh chuyển vị theo lịch sử thời gian phân tích từ OpenSees trong 2 trường hợp: có xét và không xét đến hiệu ứng P-Delta của trụ cho thấy có sự khác nhau về giá trị. Tuy nhiên, giá trị khác nhau này là nhỏ kết quả ở đây là dưới 1% (hình 14) bởi vì trụ cầu có độ mảnh nhỏ.  So sánh kết quả phân tích chuyển vị dọc cầu tại vị trí đỉnh trụ trong Midas/Civil và OpenSees (cả hai chương trình đều đã xét đến hiệu ứng P-Delta). Hình 15. So sánh kết quả giữa OpenSees và Midas/Civil 18 Tạp chí KHCN Xây dựng - số 4/2015
17. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Bảng 4. Bảng giá trị kết quả chuyển vị dọc cầu Đại lượng so sánh Kết quả từ OpenSees Kết quả từ Midas/Civil Sai số Chuyển vị lớn nhất (+) 0.101 m 0.108 m 6.5 % Chuyển vị nhỏ nhất (-) - 0.125 m - 0.117 m 6.8 % Tại vị trí đỉnh trụ trái chuyển vị theo phương dọc cầu lớn nhất khi phân tích theo OpenSees là 0.101 m và theo Midas/Civil là 0.108 m sai số là 6.5 % so với kết quả từ Midas/Civil, cùng với vị trí đó khi xét chuyển vị nhỏ nhất tại vị trí giữa nhịp kết quả phân tích theo OpenSees và Midas/Civil lần lượt là 0.125 m, 0.117 m sai số là 6.8% so với kết quả từ Midas/Civil. Từ một ví dụ bằng số mà bài báo đã phân tích cho thấy kết quả phân tích từ OpenSees có độ chính xác cao khi so với phần mềm thương mại Midas/Civil, một phần mềm được dùng chủ yếu trong phân tích kết cấu cầu. Qua quá trình phân tích bằng 2 phần mềm với cùng một ví dụ cụ thể có thể rút ra được một số ưu điểm và nhược điểm của OpenSees so với Midas/Civil như sau: Bảng 5. Bảng so sánh một số ưu nhược điểm của OpenSees với Midas/Civil Đánh giá Hạng mục Midas/Civil OpenSees Giá thành Phần mềm thương mại phải trả phí Phần mềm mã nguồn mở không mất phần mềm bản quyền. phí Đối với ví dụ trên khi sử dụng chức Trong khi đó OpenSees cung cấp năng mô hình tiết diện bê tông cốt thép hơn 200 loại mô hình vật liệu khác nhau Midas/Civil chỉ cung cấp 7 loại mô hình và cho phép người dùng tạo ra mô hình bê tông và 4 loại mô hình thép. vật liệu mới. Tính linh Trong quá trình chia thớ của tiết diện OpenSees cho phép người sử dụng Một số ưu động điểm của trụ Midas/Civil chỉ cho phép chia 1000 chia không giới hạn. OpenSees phần tử. so với Chỉ cho phép liên hợp 6 loại vật liệu Cho phép mô hình liên hợp không giới Midas/Civil trong cùng một tiết diện. hạn số loại vật liệu trong cùng 1 tiết diện. Midas/Civil không cho phép người OpenSees cho phép người dùng lựa dùng lựa chọn phương pháp giải phù chọn phương pháp giải và độ hội tụ yêu hợp với từng loại bài toán. cầu với từng bài toán cụ thể của người Phân tích sử dụng. Thời gian phân tích trong Midas/Civil Trong ví dụ trên thời gian phân tích trong cùng một ví dụ trên là 301s. trong OpenSees là 61s. Có giao diện, không mất nhiều thời Không có giao diện, tuy nhiên người Một số Sử dụng gian để sử dụng. dùng có thể lập trình tạo giao diện nhược người dùng điểm của Người dùng không cần biết quá nhiều OpenSees yêu cầu người dùng phải OpenSees về kiến thức phần tử hữu hạn cũng có lựa chọn các phương pháp giải cụ thể Tính đơn so với thể giải chính xác được các bài toán do đó người dùng cần phải nắm vững giản Midas/Civil phức tạp. được về lý thuyết và các phương pháp giải để áp dụng. Kết luận Phù hợp cho kỹ sư thiết kế Phù hợp với nghiên cứu. 4. Kết luận thương mại. Đây là phần mềm phù hợp với mục đích nghiên cứu. Tuy nhiên, OpenSees không có giao diện Sử dụng OpenSees, việc mô hình, mô phỏng kết người dùng, hiện nay đã có một số chương trình trợ cấu rất linh động đặc biệt là các bài toán kết cấu chịu giúp xử lý nhập và xuất dữ liệu nhưng vẫn còn nhiều động đất, phân tích phi tuyến và mô hình tương tác hạn chế. kết cấu đất nền. Người dùng có thể tùy chỉnh, can thiệp vào hầu hết các thông số từ phần tử, vật liệu TÀI LIỆU THAM KHẢO đến phương pháp phân tích. Đặc biệt OpenSees là [1]. Website The OpenSees Community Forum: phần mềm mã nguồn mở miễn phí khi phân tích cho kết quả chính xác không thua kém các phần mềm Tạp chí KHCN Xây dựng - số 4/2015 19
18. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG [2]. WaiChing Sun (2004): Earthquake Engineering Research Center, University of California at Berkeley. ange/6113-opensees-pre-and-post-processing. [8]. Website: [3]. Jinchi Lu, Kevin Mackie, and Ahmed Elgamal. (2011). “OpenSees 3D Pushover and Earthquake er/download.php. Analysis of Single-Column 2-span Bridges”, UC, [9]. Website: Berkeley ( [4]. P. Gavali, M S. Shah, G. Kadam và K. Meher. _motion_db.html. (2013). ”Seismic response and simulations of [10]. B.D.Scott, R.Park, and M.J.Priestley. (1982) reinforced concrete bridge using OpenSees on “Stress-strain Behavior of concrete confined by high performance computing”, CSI Transactions overlapping hoop at low and high strain rates”, ACI on ICT, Volume 1, Issue 3, pp 215-220. Journal, January-February 1982, title no. 79-2. [5]. Website: [11]. Finley A. Charney. (April 2008). "Unintended Consequences of Modeling Damping in hnjs.htm. Structures", J. Struct. Engrg. Volume 134, Issue [6]. Website: 4, pp. 581-592. programming-interface. Ngày nhận bài: 05/8/2015. [7]. McKenna F, Fenves GL. (2001). The OpenSees Ngày nhận bài sửa lần cuối: 03/9/2015. command language manual, Version 1.2, Pacific 20 Tạp chí KHCN Xây dựng - số 4/2015