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

ầu tính lặp trong mỗi bước. Trong bài báo này, tác giả đề xuất một họ phương pháp phân tích động phi tuyến mới. Họ phương pháp này, tuy là ngoại hiển thức nhưng lại không có điều kiện ổn định. Phương pháp này còn có hệ số tiêu tán thích hợp và có thể kiểm soát được, có thể điều chỉnh để hệ số cản nhớt số bằng không. Ưu điểm lớn nhất của phương pháp này là không cần tính lặp trong mỗi bước, do vậy tiết kiệm được rất nhiều công sức tính toán so với các phương pháp nội ẩn thức hiện có. 1. Đặt vấn đề Về cơ bản, các bài toán về dao động của hệ kết cấu được chia làm hai dạng chính: Dạng thứ nhất là các hệ kết cấu bị chi phối bởi lực quán tính, chiếm đa số trong các bài toán động lực học công trình, dao động của dạng kết cấu này chủ yếu ảnh hưởng bởi các dạng dao động có tần số thấp; Dạng thứ hai là các hệ kết cấu bị chi phối bởi cả các dạng dao động có tần số thấp và tần số cao, ví dụ như khi hệ kết cấu công trình bị tác động bởi lực gây ra do nổ và va chạm. Phương pháp để giải các bài toán thuộc dạng thứ nhất thường được đề xuất là phương pháp nội ẩn thức (implicit) như các phương pháp [1, 2, 3, 9, 11, 12, 13, 16]. Trong khi đó, ngoại hiển thức (explicit) là phương pháp thích hợp để giải các bài toán dạng thứ hai, ví dụ như phương pháp của Newmark [13]. Điều này là do các phương pháp nội ẩn thức không có điều kiện ổn định do vậy có thể sử dụng các bước thời 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 quả bài toán. Ngược lại, ưu điểm chính của các phương pháp ngoại hiển thức là các giá trị của bước 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 các bài toán phi tuyến, do vậy khối lượng tính toán trong một bước sẽ ít hơn nhiều so với phương pháp nội ẩn thức. Khi áp dụng phương pháp ngoại hiển thức cho các bài toán cần quan tâm đến các dạng dao động bậc cao, nếu giá trị của các bước thời gian tính toán được lựa chọn để thỏa mãn điều kiện ổn định thì độ chính xác của kết quả cũng sẽ tự động được thỏa mãn. Đã có một số các phương pháp tính tích phân phụ thuộc kết cấu được đề xuất bởi Chang [4, 5, 6, 7] nhằm tích hợp đồng thời các ưu điểm của hai phương pháp nội ẩn thức và ngoại hiển thức, các phương pháp này đều có đặc điểm không có điều kiện ổn định (giống ưu điểm như phương pháp nội ẩn thức) và không cần tính lặp phi tuyến (giống ưu điểm của phương pháp ngoại hiển thức). Các phương pháp này rất phù hợp để giải các bài toán thuộc dạng thứ nhất. Tuy nhiên, khi giải các bài toán thuộc dạng thứ hai thì các phương pháp đó không có hệ số tiêu tán (numerical dissipation) phù hợp để loại bỏ nhiễu do ảnh hưởng bởi các dạng dao động bậc cao. Trong các phương pháp được biết đến rộng rãi hiện nay, có một số họ phương pháp tính toán đã được phát triển có hệ số tiêu tán thích hợp, như các phương pháp trong tài liệu [1, 10, 15, 16, 17], nhưng các phương pháp này đều thuộc nhóm phương pháp nội ẩn thức, do vậy đều cần tính lặp phi tuyến khi giải các bài toán kết cấu phi tuyến. Vì những lý do trên, phương pháp được đề xuất trong bài báo này là một phương pháp ngoại hiển thức không có điều kiện ổn định và không cần tính lặp phi tuyến, giống các phương pháp đã được đề xuất trước đó [4, 5, 6, 7]. Điểm khác là phương pháp này còn có một hệ số tiêu tán phù hợp, có thể điều chỉnh, dùng để giải các bài toán thuộc dạng thứ hai. 2. Phương pháp đề xuất Phương trình dao động của hệ một bậc tự do được viết như sau:    mu cu ku f (1) trong đó: m, c, k, f lần lượt là khối lượng, độ cản nhớt vật lý, độ cứng và ngoại lực tác dụng lên hệ kết cấu; KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 4 Tạp chí KHCN Xây dựng - số 4/2015 u , u và u tương ứng là các giá trị chuyển vị, vận tốc, gia tốc. Phương pháp đề xuất để giải phương trình dao động (1) được viết như sau:             2 1 2 1 1 0 1 1 1 11 1 1 1 2 1 0 1 1 2 3 3 1 3 1 12 1 2 1 p p p p ma c v k d k d f fi i i i i i i ip p p p d d d t v t ai i i i i p p v v t a t ai i i ip p                                   (2) trong đó: ai+1, vi+1, di+1, fi+1 là gia tốc, vận tốc, chuyển vị và lực tác dụng lên hệ kết cấu ở bước thứ (i+1); ai, vi, di, fi là gia tốc, vận tốc, chuyển vị và ngoại lực tác dụng ở bước (i); 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); c0 là độ cản nhớt vật lý của hệ kết cấu (giả thiết c0 không thay đổi trong suốt quá trình tính toán, điều này thường đúng với phần lớn các hệ kết cấu); Δt là bước thời gian tính toán. Các hệ số β0 đến β3 được tính như sau: 3 1 1 2 2 0 08 1 3 1 1 2 211 08 1 1 3 12 01 2 1 1 1 2 3 3 02 2 1 1 p D p p D p p D p p D p p                                                                   (3) trong đó: ξ là hệ số cản nhớt; Ω0 =ω0(Δt) với 0 0 /k m  là tần số dao động tự nhiên của hệ kết cấu tương ứng với độ cứng ban đầu k0; 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 trình bày rõ hơn ở mục 3). Hệ số D được tính như sau: 3 3 2 21 0 01 4 1 p p D p p                (4) Để tiện cho việc tính toán, các hệ số ξΩ0 và Ω02 được viết lại ξΩ0 = ξω0(Δt) = c0(Δt)/(2m) và Ω02 = (Δt)2(k0/m), theo đó, các hệ số β0 đến β3 và D trở thành:                 3 21 1 2 0 08 1 3 21 1 2 11 08 1 1 3 2 02 1 2 1 1 1 2 3 3 02 4 1 1 3 23 2 0 02 1 4 1 p t k D p p t k D p p m t c D p p m t c D p p p p D m t c t k p p                                                                               (5) 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 hệ kết cấu, nó không thay đổi trong suốt quá trình tính toán, do vậy với mỗi bài toán chỉ cần tính duy nhất một lần. Trong dòng thứ hai của công thức (1), giá trị di+1 đượ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 đó phương pháp đề xuất ở đây gọi là phương pháp phụ 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, các hệ số β0 đến β3 đều là các hằng số cố định. Ngoài 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ị ở bước (i-1), do vậy, khi áp dụng phương pháp này để tính toán cần có một chút lưu ý khi tính bước đầu tiên, điều này sẽ được nói rõ ở mục 4. 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 Zhou & Tamma cũng có độ chính xác tương 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à phương pháp ngoại hiển thức. 3. Khảo sát tính chất của phương pháp đề xuất Khi khảo sát tính chất của một phương pháp phân 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 (stability), độ chính xác (accuracy) và tính hiệu quả (efficiency) của phương pháp đó. Các tính chất này đượ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 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 5 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: - NEM (Newmark Explicit Method): phương pháp ngoại hiển thức Newmark [13]; - AAM (Average Acceleration Method): phương pháp nội ẩn thức gia tốc trung bình [13]; - PFM1 (Proposed Family Method): phương pháp đề xuất với trường hợp p=1; - PFM2 (Proposed Family Method): phương pháp đề xuất với trường hợp p= 0,5. 3.1 Lựa chọn khoảng giá trị cho tham số p Việc lựa chọn khoảng giá trị p cho phương pháp này được căn cứ vào khoảng giá trị mà ở đó kết quả tính của phương pháp này là hội tụ. Muốn khảo sát 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 kính phổ (spectral radius). Một quy trình chung cho cách khảo sát này được trình bày trong tài liệu [1, 10, 11] hoặc trình bày chi tiết hơn cho phương pháp này trong tài liệu [8]. Kết quả khảo sát cho thấy với phương pháp này khoảng giá trị thích hợp nhất của tham số p là 0,5 ≤ p ≤1, trong đó với p = 1 cho giá trị bán kính phổ luôn bằng 1, chứng tỏ hệ số cản nhớt số (xem mục 3.3) sẽ bằng 0. Qua hình 1 ta thấy bán kính phổ cũng bằng 1 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 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 3.2 Sai số tương đối của chu kỳ Sai số tương đối của chu kỳ (relative period error) là đại lượng được định nghĩa bằng   /T T T , trong đó T là chu kỳ dao động chính xác của hệ, T là chu kỳ dao động tính toán của hệ. Đại lượng này đặc trưng cho độ chính xác của phương pháp phân tích. Đại lượng này càng nhỏ thì phương pháp phân tích càng chính xác. Hình 2 biểu diễn mối quan hệ giữa sai số tương đố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 chu kỳ của phương pháp AAM cũng được thể hiện trong hình vẽ để so sánh. Hình 2 cũng cho thấy sai số tương đối của chu kỳ tỷ lệ nghịch với giá trị của p, p càng giảm thì sai số 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 cho thấy rằng với các giá trị Δt/T nhỏ thì sai số cũng nhỏ. Với giá trị Δt/T ≤ 0,1 thì sai số của phương pháp này là dưới 5%. Với trường hợp p = 1, đường cong sai số trùng với đường cong sai số của phương pháp AAM, nói cách khác, độ chính của phương pháp PFM1 tương đương với độ chính xác của phương pháp AAM. 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 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 6 Tạp chí KHCN Xây dựng - số 4/2015 3.3 Hệ số cản nhớt số Như nói ở phần đặt vấn đề, phần lớn các bài toán 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 thấp mà không quan tâm đến ảnh hưởng của những dạng dao động bậc cao, do ảnh hưởng của dạng dao động bậc cao đến tổng thể kết cấu là không lớn, hơn nữa, kết quả tính cho những dạng dao động này thường kém chính xác nên để đơn giản người ta sẽ loại bỏ nó đi. Hệ số cản nhớt số (numerical damping ratio) của mỗi phương pháp tính là đại lượng đặc trưng cho khả năng loại bỏ sự ảnh hưởng của các dạng dao động bậc cao mà không làm ảnh hưởng đến độ chính xác trong kết quả tính toán của các dạng dao động bậc thấp. 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 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ể hiện với các trường hợp tương ứng với p = 1,0; 0,75, 0,5 cũng như với phương pháp AAM để tham khảo. Biểu đồ cho thấy, với trường hợp p = 0,75 và p = 0,5 các dạng dao động bậc cao (Δt/T) sẽ dễ dàng bị loại 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 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 nhớt số, được dùng để tính toán với các bài toán có xét đến đến ảnh hưởng của cả các dạng dao động bậc thấp và bậc cao. 3.4 Ảnh hưởng của dao động bậc cao Để làm rõ hơn đặc điểm về hệ số cản nhớt số của 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ữ tiếng Anh với tên gọi là overshooting. Xét một hệ một bậc tự do có khối lượng m = 1kg và độ cứng k = 1N/m, chu kỳ dao động tự do của hệ sẽ là: 2 6,28 sT m k  . Cho hệ dao động tự do từ trạng 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 tiên được thể hiện trong hình 4. Hình 4 cho thấy đường cong biểu diễn kết quả tính theo phương pháp PFM1 và AAM hoàn toàn trùng lên nhau, thêm vào đó chuyển vị cũng không bị suy giảm. Trong khi đó, giá trị chuyển vị tính theo phương pháp PFM2 bị tắt rất nhanh chỉ sau 10 bước đầu tiên, đó là do phương pháp PFM2 có hệ số cản nhớt số rất lớn. Với PFM2, giá trị chuyển vị bị vượt quá không đáng kể trong vài bước đầu tiên, sau đó tắt rất nhanh, trong khi giá trị v0/ω0 bị vượt quá (overshoot) đáng kể trong những bước tính toán đầu tiên, tuy nhiên sau đó tắt rất nhanh nên xét về lâu dài 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 4. So sánh ảnh hưởng của dao động bậc cao 4. Áp dụng với hệ nhiều bậc tự do Phần lớn các bài toán động lực học công trình là các bài toán hệ có nhiều bậc tự do, do đó, phần này sẽ đưa ra cách áp dụng phương pháp đề xuất ở đây để giải các bài toán dạng này. Với hệ nhiều bậc tự do, công thức (2) sẽ được viết như sau:             2 1 2 1 1 0 1 1 11 1 1 1 2 1 0 1 1 2 3 3 1 3 1 12 1 2 1 p p p p i i i i i ip p p p t ti i i i i p p t ti i i ip p                              a C v R R f f d B d B d B v B a v v a a M (6) trong đó: M, C0 là các ma trận khối lượng, ma trận độ cản nhớt vật lý của kết cấu; a, v, d, f tương ứng là các vec-tơ gia tốc, vận tốc, chuyển vị và ngoại lực tác dụng; KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 7 Các hậu tố (0), (i-1), (i), (i+1) chỉ thứ tự các bước tính toán; Ri, Ri+1 là vec-tơ nội lực của hệ kết cấu ở bước tính toán thứ (i) và (i+1). Ma trận B0 đến B3 được tính như sau:                 3 21 21 0 08 1 3 21 21 1 08 1 31 2 02 1 2 2 31 1 1 423 01 1 3 23 2 0 02 1 4 1 p t p p t p p t p p t p p p p t t p p                                                                     B D K B I D K B D M C B D M C D M C K (7) 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 đầu (initial stiffness). 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, C0, K0 và 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 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 đó, do vậy cần có một bước đệm khi tính toán với 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 tiên. Với các bước tiếp theo, giá trị hệ số p được lựa chọn tùy theo yêu cầu tính toán. Quy trình tính toán với hệ nhiều bậc tự do được thực hiện như sau: Bước 1: Tính giá trị vec-tơ chuyển vị di+1 từ dòng thứ hai của công thức (6); Bước 2: Thế giá trị vec-tơ chuyển vị di+1 vừa tìm được và giá trị vec-tơ vận tốc vi+1 ở dòng thứ ba vào dòng thứ nhất cùng của công thức (6), giải và tìm vec- tơ gia tốc ai+1; Bước 3: Giá trị vec-tơ vận tốc được tính bằng dòng thứ ba của công thức (6) sau khi đã có giá trị vec-tơ vận tốc vi+1. 5. Một số ví dụ tính toán Để làm rõ hơn các đặc điểm của phương pháp này, một số ví dụ sẽ được trình bày ở đây. Các ví dụ 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 được biết đến rất rộng rãi là NEM [13], đại diện cho phương pháp ngoại hiển thức và AAM [13], tiêu biểu cho phương pháp nội ẩn thức. Nhìn chung, với tất cả các phương pháp phân tích động phi tuyến đều cần sử dụng máy tính do khối lượng tính toán lớn. Với những bài toán đơn giản, người dùng có thể tự lập trình bằng phần mềm như Excel hoặc viết các chương trình dựa trên ngôn ngữ lập trình Fortran, Matlab, C++ Với các bài toán phức tạp hơn thì nên sử dụng những phần mềm chuyên dụng như Sap, Etabs hay OpenSees. Các ví dụ tính toán ở đây được tác giả tự lập trình bằng Matlab. 5.1 Ví dụ 1: Hệ một bậc tự do đàn dẻo tuyệt đối m =104 kg k =106 N/m di ag Hình 5. Mô hình thí nghiệm trên bàn rung với hệ một bậc tự do Giả thiết có một hệ một bậc tự do được thí nghiệm trên bàn rung như hình 5. Tải trọng tập trung ở đầu cột bằng m = 104kg, cột giả thiết như một thanh đàn dẻo tuyệt đối, độ cứng k = 106N/m, do đó chu kỳ dao động tự nhiên ban đầu của hệ ω0 = 10 rad/s. Cường độ chịu kéo và chịu nén của vật liệu thanh giả thiết bằng Rt = Rc = 50kN. Bỏ qua hệ số cản nhớt vật lý của hệ. Gia tốc nền tác dụng lên hệ được điều 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 Chi-Chi xảy ra ở Đài Loan vào năm 1999 (tên phổ ghi gia tốc này theo ký hiệu quốc tế thường dùng là CHY 028), đỉnh gia tốc nền được lấy bằng 0,5g. Dùng phương pháp PFM1 để dự đoán chuyển vị của hệ đồng thời so sánh kết quả với hai phương pháp NEM và AAM. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 8 Tạp chí KHCN Xây dựng - số 4/2015 Hình 6 thể hiện kết quả tính toán chuyển vị theo thời gian của hệ và biểu đồ quan hệ giữa chuyển vị và nội lực thanh. Giá trị bước thời gian tính toán theo 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, làm chuẩn để so sánh với các phương pháp khác, trong khi bước thời gian tính toán theo AAM và PFM1 là Δt =0,02 s. 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, mỗi tầng được mô hình hóa như một tấm cứng tuyệt đối, toàn bộ hệ kết cấu 5 tầng được coi như có 5 bậc tự do. Các giá trị tải trọng và độ cứng mỗi tầng được 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ể hiện trong hình vẽ. Độ cứng của mỗi tầng được thể hiện qua công thức:  1/50 11 , 1~ 5j i i i i ik k p u u i        (8) trong đó: kj-i là độ cứng tức thời của tầng thứ i tại thời điểm cuối cùng của bước thứ j; k0-i là độ cứng ban đầu (trước khi biến dạng) của tầng thứ i; |ui – ui-1| là chuyển vị lệch tầng của tầng thứ i so với tầng thứ (i-1) bên dưới; Hệ số pi đặc trưng cho tính phi tuyến của mỗi tầng, pi = 0 tương ứng với thanh đàn hồi, pi < 0 tương ứng với hệ mềm hóa sau biến dạng (softening), pi > 0 tương ứng với hệ cứng hóa sau biến dạng (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 giữa chuyển vị lệch tầng và nội lực trong cột, đường nét đứt thể hiện mối quan hệ phi tuyến, ở đây là quan hệ mềm hóa. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 9 Hình 8. Quan hệ giữa chuyển vị lệch tầng và nội lực trong cột Trong ví dụ này, mô hình công trình chịu tác dụng của lực động đất, chuyển động với gia tốc nền CHY 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 tính từ trạng thái nghỉ d(0) = 0 được in trong hình 9a để so sánh. Tất cả các trường hợp này đều tính với bước thời gian Δt = 0,01s. Kết quả từ hình 9 cho thấy chuyển vị tại tầng 5 của hệ kết cấu bị ảnh hưởng rất nhiều của các dạng dao động bậc cao khi tính bằng phương pháp AAM và PFM1. Các giá trị dao động trở nên rất “nhiễu” (không chính xác) và gần như không thể sử dụng được để phân tích kết cấu. So sánh hình 9a và 9b cho thấy, độ “nhiễu” theo tính toán bằng phương pháp PFM1 thậm chí còn lớn hơn độ “nhiễu” theo tính toán bằng AAM. Tuy nhiên, độ “nhiễu” này đã hoàn toàn bị loại bỏ trong kết quả của phương pháp PFM2. Kết quả chuyển vị thu được từ hình 9c cho thấy đường cong biểu diễn chuyển vị rất trơn và trùng với kết quả thu được từ phương pháp AAM trong trường hợp hệ dao động từ trạng thái nghỉ. 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 có thể điều chỉnh được hệ số cản nhớt số để loại bỏ ảnh hưởng của các dạng dao động bậc cao mà không làm ảnh hưởng đến độ chính xác của kết quả tính toán. Với giá trị p = 1, phương pháp này không có hệ số cản nhớt số, giá trị p giảm cho hệ số cản nhớt số càng tăng, và hệ số cản nhớt số lớn nhất khi p = 0,5. Hình 9. Dao động tầng 5 của hệ kết cấu dưới tác dụng của tải trọng động đất 5.3 Ví dụ 3: So sánh về hiệu quả tính toán Hình 10. Mô hình hệ nhiều bậc tự do Xem xét mô hình một hệ phi tuyến nhiều bậc tự do như trong hình 10. Hệ gồm n vật nặng, mỗi vật có khối lượng 1kg, nối với nhau bằng lò xo có độ cứng ki, trong đó giá trị ki thay đổi phụ thuộc vào biến dạng của lò xo, có công thức tính được ghi trong hình 10. Toàn bộ hệ chịu tác dụng của chuyển vị nền, dao động theo dạng hình sin như trong hình vẽ. Xem xét hai trường hợp, hệ có 500 bậc tự do (n = 500) và hệ có 1000 bậc tự do (n = 1000). Về nguyên tắc thì mô hình hệ lò xo này cũng giống mô hình hệ 5 tầng 5 bậc tự do trong ví dụ 2 (khác số bậc tự do). Chu kỳ dao động tự nhiên nhỏ nhất được tính toán dựa vào độ cứng ban đầu trước khi biến dạng với hệ 500 bậc tự do có kết quả bằng ω0 = 31,38 rad/s, trong đó với hệ 1000 bậc tự do bằng ω0 = 15,70 rad/s, trong khi chu kỳ dao động tự nhiên lớn nhất với cả hai trường hợp đề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. Bước thời gian tính toán Δt đối với phương pháp NEM được lựa chọn căn cứ theo điều kiện ổn định Ω0 = ω0(Δt)≤2, theo đó với chu kỳ dao động lớn nhất ω0 = 20000 rad/s thì bước thời gian được lựa chọn Δt ≤ 2/20000 = 0,0001s. Giá trị tính toán của bước thời gian này cũng đủ nhỏ để có thể coi kết quả tính toán thu được là chính xác, làm cơ sở để so sánh với KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 10 Tạp chí KHCN Xây dựng - số 4/2015 các phương pháp khác. Trong khi đó, do hai phương pháp AAM và PFM2 đều là những phương pháp không có điều kiện ổn định nên bước thời gian Δt được lựa chọn cho hai phương pháp AAM và PFM2 căn cứ vào yêu cầu chính xác của kết quả tính toán với Δt càng nhỏ thì độ chính xác càng cao, nhưng ngược lại khối lượng tính toán càng lớn. Ở đây lựa chọn Δt = 0,005s giây cho cả hai phương pháp AAM và PFM2. Hình 11 thể hiện kết quả kết quả tính toán chuyển vị của lò xo thứ 500 và thứ 1000 tương ứng với hai hệ đã nêu ở trên. Biểu đồ cho thấy đường cong biểu diễn chuyển vị gần như trùng nhau, cho thấy kết quả tính toán từ phương pháp AAM và PFM2 là chính xác. Tiếp theo so sánh về thời gian tính toán, cả ba phương pháp này đều được tính toán bằng chương trình viết bằng Matlab, chạy trên máy tính cá nhân dùng chip Intel®CoreTMi5 CPU M460 @2.53GHz, RAM 4.00 GB. Thời gian tính toán (tính bằng giây) cho từng hệ 500 và 1000 bậc tự do tương ứng với mỗi phương pháp được trình bày trong bảng 1. Kết quả so sánh từ bảng 1 cho thấy, thời gian tính toán của phương pháp PFM2 chỉ bằng khoảng 5% so với phương pháp NEM và bằng dưới 2% so với phương pháp AAM. Điều này được lý giải là do trong một bước tính toán thì tính bằng phương pháp PFM2 mất nhiều công sức hơn (PFM2 phải giải hệ phương trình ma trận hai lần, NEM chỉ cần giải một lần). Tuy nhiên, do NEM bị giới hạn về điều kiện ổn định nên số bước tính toán bị tăng lên rất nhiều, ở đây số bước tính toán của phương pháp NEM lớn hơn phương pháp PFM2 tới 50 lần. Do vậy, xét về tổng thể thì thời gian thời gian tính toán bằng NEM vẫn nhiều hơn PFM2. Mặt khác, khi so sánh với phương pháp AAM, tuy PFM2 và AAM có số bước tính toán bằng nhau 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 pháp Newton Raphson thường mất trung bình khoảng 10-25 lần lặp trong một bước). Việc tính lặp này rất mất thời gian, công sức và kết quả không phải lúc nào cũng hội tụ hoặc hội tụ đến kết quả không mong muốn. Do vậy, trong một bước tính toán thì tính bằng AAM mất thời gian và công sức hơn PFM2 khá nhiều, tùy thuộc vào tốc độ hội tụ và độ chính xác yêu cầu. Trong phân tích ứng xử phi tuyến cho công trình xây dựng, việc giảm thời gian tính toán có ý nghĩa rất lớn. Để giảm thời gian tính toán có thể sử dụng những hệ máy tính mạnh hoặc chạy song song nhiều máy tính, tuy nhiên cách này thường gây tốn kém chi phí, do đó sử dụng một phương pháp có hiệu quả tính toán cao như phương pháp đề xuất ở đây cũng là một cách tốt. Trong ví dụ này chỉ tính toán với hệ 1000 bậc tự do, với hệ có số bậc tự do càng lớn, chẳng hạn như mô hình một công trình xây dựng có thể lên đến hàng triệu bậc tự do, thì việc tiết kiệm này sẽ càng có ý nghĩa. Hình 11. Kết quả tính toán chuyển vị KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 11 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 Trong bài báo này, một họ phương pháp phân tích phi tuyến động theo lịch sử thời gian được đề xuất, trong đó các đặc trưng số học của phương pháp được điều chỉnh chỉ thông qua một hệ số p duy nhất. Hệ số này được lựa chọn trong khoảng 0,5 ≤ p ≤ 1, với p = 1 cho phép tính toán không có hệ số cản nhớt số, p giảm cho hệ số cản nhớt số tăng. Phương pháp thuộc họ ngoại hiển thức không có điều kiện ổn định, do vậy tiết kiệm được rất nhiều công sức tính toán so với hai phương pháp truyền thống được so sánh là NEM và AAM, trong khi độ chính xác trong kết quả thu được là tương đương. Thông qua các ví dụ đã chỉ ra rằng phương pháp này là phù hợp để giải các bài toán phi tuyến trong xây dựng. Lời cảm ơn Nhóm tác giả xin cảm ơn Hội đồng Khoa học Quốc gia Đài Loan đã tài trợ kinh phí cho nghiên cứu này thông qua hợp đồng số NSC-99-2221-E-027-029. Bài báo được gửi đăng với sự đồng ý của ©2015 Taylor & Francis Group, London, UK. TÀI LIỆU THAM KHẢO [1]. Bathe, K.J. and Wilson, E.L. (1973). “Stability and accuracy analysis of direct integration methods”, Earthquake Engineering and Structural Dynamics. Vol. 1: 283-291. [2]. Belytschko, T. and Schoeberle, D.F. (1975). “On the unconditional stability of an implicit algorithm for nonlinear structural dynamics”, Journal of Applied Mechanics, Vol. 17: 865-869. [3]. Belytschko, T. and Hughes, T.J.R. (1983). “Computational methods for transient analysis”, Elsevier Science Publishers B.V., North-Holland. [4]. Chang, S.Y. (2002). “Explicit pseudodynamic algorithm with unconditional stability”. Journal of Engineering Mechanics, ASCE, Vol. 128, No.9: 935-947. [5]. Chang, S.Y. (2007). “Improved explicit method for structural dynamics”, Journal of Engineering Mechanics, ASCE, Vol. 133, No. 7: 748-760. [6]. Chang, S.Y. (2009). “An explicit method with improved stability property”, International Journal for Numerical Method in Engineering, Vol. 77, No. 8:1100-1120. [7]. Chang, S.Y. (2010). “A new family of explicit method for linear structural dynamics”, Computers & Structures, Vol. 88, No.11-12:755-772. [8]. Chang, S.Y., N.C. Tran. (2014). “A two-step unconditionally stable explicit method with controllable numerical dissipation for structural dynamics”, Advances in Civil Engineering and Building Materials IV. 379-383. [9]. Dobbs, M.W. (1974). “Comments on ‘stability and accuracy analysis of direct integration methods’ by Bathe & Wilson”, Earthquake Engineering and Structural Dynamics, Vol. 2: 295-299. [10]. Hilber, H.M., Hughes, T.J.R. and Taylor, R.L. (1977). “Improved numerical dissipation for time integration algorithms in structural dynamics”, Earthquake Engineering and Structural Dynamics, Vol. 5: 283-292. [11]. Hughes, T.J.R. (1987). “The Finite element method”, Prentice-Hall, Inc., Englewood Cliffs, N.J., U.S.A. [12]. Krieg, R.D. (1973). “Unconditional stability in numerical time integration methods”, Journal of Applied Mechanics, Vol. 40: 417-421. [13]. Newmark, N.M. (1959). “A method of computation for structural dynamics”, Journal of Engineering Mechanics Division, ASCE, Vol. 85: 67-94. [14]. Zhou X, Tamma KK. (2004). “Design, analysis and synthesis of generalized single step single solve and optimal algorithms for structural dynamics”, International Journal for Numerical Methods in Engineering; Vol. 59: 597-668. [15]. Zhou X, Tamma KK. (2006). “Algorithms by design with illustrations to solid and structural mechanics/ dynamics”, International Journal for Numerical Methods in Engineering; Vol. 66: 1841-1870. [16]. Zienkiewicz, O.C. (1977). “The Finite Element Method”, McGraw-Hill Book Co (UK) Ltd. Third edition. [17]. W.L. Wood, M. Bossak, and O.C. Zienkiewicz. “An alpha modification of Newmark’s method”, International Journal for Numerical Methods in Engineering, Vol. 15: 1562-1566, 1981. Ngày nhận bài: 03/8/2015. Ngày nhận bài sửa lần cuối: 06/11/2015. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 12 Tạp chí KHCN Xây dựng - số 4/2015 Ứ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 Engineering Simulation) là một phần mềm mã nguồn mở với thư viện mã code được viết chủ yếu bằng C++, một ngôn ngữ lập trình hướng đối tượng. Nó cho phép người dùng tạo ra các mô hình phần tử hữu hạn để mô phỏng phản ứng của hệ kết cấu và đất nền dưới tác dụng của động đất. Bài báo này sẽ giới thiệu về OpenSees từ đó xây dựng thuật toán ứng dụng OpenSees vào trong lập trình mô phỏng một ví dụ kết cấu cầu chịu tác dụng của động đất. Từ khóa: OpenSees, động đất, mô phỏng. 1. Mở đầu OpenSees (Open System of Earthquake Engineering Simulation) là một phần mềm mã nguồn mở được phát triển bởi nhóm nghiên cứu ở PEER (Pacific Earthquake Engineering Research) của trường đại học University of California, Berkeley. Thư viện mã code của OpenSees được viết chủ yếu bằng C++ là ngôn ngữ lập trình hướng đối tượng. OpenSees là một phần mềm mã nguồn mở miễn phí. Phần mềm có khả năng tính toán mạnh với một thư viện hỗ trợ nhiều loại vật liệu, người dùng có thể tự phát triển các loại vật liệu khác tùy thuộc vào nhu cầu của người sử dụng. OpenSees hỗ trợ nhiều phương pháp tính khác nhau, do đó tùy vào mục đích và đặc điểm của bài toán kết cấu mà người dùng có thể linh hoạt lựa chọn phương pháp thích hợp để 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 nghiên cứu lớn trên thế giới, các chuyên gia từ các trường đại học, viện nghiên cứu, đây là phần mềm có lượng người dùng lớn do đó việc trao đổi, thảo luận là tương đối dễ dàng [1]. Tuy nhiên, OpenSees có nhược điểm là một phần mềm không có giao diện đồ họa, vì vậy gây khó khăn cho những người dùng thông thường, hiện nay có một số phần mềm đã được viết như OpenSees Navigator, GiD, TclBuilder hoặc toolbox trong Matlab như OpenSees Pre- and Post- Processing [2] để xuất và nhập kết quả. Tuy nhiên, việc sử dụng các phần mềm hỗ trợ này vẫn chưa đầy đủ, hoàn chỉnh như các phần mềm thương mại hiện nay. Bởi vì OpenSees hỗ trợ nhiều phương pháp giải khác nhau, do đó đối với người dùng thông thường sẽ rất khó khăn trong việc lựa chọn, thậm chí dẫn đến tính toán đưa ra kết quả sai nếu người dùng không hiểu hoặc 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 công tác nghiên cứu và phát triển chứ chưa thực sự hướng đến những người dùng phổ thông. Sử dụng OpenSees trong phân tích kết cấu xây dựng nói chung và đặc biệt trong kết cấu cầu nói riêng được nhiều nhóm nghiên cứu quan tâm, Jinchi Lu, Kevin Mackie, and Ahmed Elgamal [3] đã xây dựng chương trình phân tích Pushover cho kết cấu nhịp giản đơn; Pallavi Gavali, Mahesh S. Shah, Gouri Kadam và Kranti Meher [4] đã xây dựng chương trình phân tích 3D kết cấu cầu và đất nền chịu động đất. OpenSees hiện tại và trong tương lai gần vẫn là một công cụ mạnh trong phát triển và mô phỏng kết cấu chịu động đất. Bài báo này sẽ giới thiệu về OpenSees từ đó xây dựng sơ đồ thuật toán ứng dụng vào trong lập trình mô phỏng một ví dụ kết cấu cầu cụ thể chịu tác dụng của động đất. 2. Tổng quan về OpenSees 2.1 Những tính năng của OpenSees Như đã đề cập trong phần mở đầu, OpenSees là một phần mềm mã nguồn mở phục vụ cho phân tích tính toán kết cấu theo phương pháp phần tử hữu hạn. Tính năng mạnh nhất của OpenSees là mô phỏng phản ứng của kết cấu chịu tác động của động đất, các bài toán phân tích phi tuyến và phân tích kết cấu tương tác với đất nền. Sơ đồ hình 1 thể hiện các ưu điểm mà OpenSees mang lại: KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 13 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. 2.2 Mô hình trong OpenSees OpenSees bao gồm một tập hợp các mô-đun để thực hiện công việc tạo ra các mô hình phần tử hữu hạn. Trong một mô hình phân tích phần tử hữu hạn của kết cấu trong OpenSees, các lệnh được sử dụng với mục đích để tạo ra 4 loại đối tượng chính (hình 3). Hình 3. Lớp trừu tượng chính (main abstractions) trong OpenSees [7] Các loại đối tượng chính này thực hiện các nhiệm vụ khác nhau, cụ thể như: ModelBuilder: Là một đối tượng trong chương trình chịu trách nhiệm tạo ra mô hình kết cấu như nút, 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. KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 14 Tạp chí KHCN Xây dựng - số 4/2015 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 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 15 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 10 00 0 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 đất Hướng Trạm đo Cấp động đất (Mw) Bước thời gian (s) PGA (g) El Centro12 Phương dọc Cầu (Trục X hệ tọa độ tổng thể) Imperial valley 10/15/79 EL CENTRO ARRAY #12 6.5 0.005 0.116 El Centro01 Theo phương ngang Cầu (Trục Y hệ tọa độ tổng thể) Imperial valley 10/15/79 EL CENTRO ARRAY #1 6.5 0.005 0.139 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 16 Tạp chí KHCN Xây dựng - số 4/2015 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 Dầm chủ Mô hình 3D phần tử dầm – cột đàn hồi Diện tích tiết diện: 6.9 m 2 Mô men quán tính: 1.31 m4 Trụ cầu 3Mô hình 3D phần tử dầm – cột phi tuyến, phân tích ảnh hưởng hiệu ứng P – Delta. Diện tích tiết diện: 1.7 m2 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σ = Kf 2 -c c ε εco co                 ε εc o (1)  'σ = Kf 1- Z(ε - ε )c c c co co c uε ε ε  (2) 'σ = 0.2Kfc c cu cε ε (3) ρ fs yh K = 1+ 'fc coε = 0.002K (4) Hình 10. Mô hình bê tông [10] 0.5 Z = ' '3 + 0.0284f hc + 0.75ρ - 0.002Ks' s14.21f - 1000c h (5) Ở đâ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 chế; Z là độ dốc đường biến dạng; 'cf là cường độ chịu nén của bê tông kg/cm 2; fyh 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 Bê tông Ec (Mpa) f ’ c (Mpa) γ (kg/m3) ε0 εu K Z H.lượng thép dọc ρ (%) Cốt thép đ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 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 17 Mô hình cốt thép: Quan hệ ứng suất – biến dạng của cốt thép được mô hình đơn giản là các đoạn thẳng tuyến tính (hình 11). Lý do sử dụng mô hình xấp xỉ này là tiện lợi cho mô hình tính toán. Mô hình này đã được OpenSees đưa vào hàm Steel01 [7]. Số liệu cốt thép đưa vào ví dụ như sau: Cường độ chảy của cốt thép: fy = 420 MPa; mô đun đàn hồi của thép: E0 = 200000 MPa; tỷ số biến cứng: b = 0.01; trọng lượng riêng thép: γ = 7500 kg/m3. Hình 11. Mô hình vật liệu thép [7] 3.3 Cài đặt các thông số phân tích cho mô hình trong OpenSees Để thực hiện phân tích động phi tuyến, phương pháp tích phân trực tiếp Newmark với Δt = 0.01s trên cơ sở thuật toán gia tốc trung bình được sử dụng, phương pháp này đã được lập sẵn trong hàm integrator Newmark $gamma $beta. Thuật giải lặp Newton-Raphson được khai báo trong OpenSees bằng hàm algorithm Newton. Để thiết lập ma trận cản tỉ số cản được giả thiết là ξ = 5% cho 2 tần số dao động đầu tiên [11]. Để tiết kiệm bộ nhớ máy tính khi phân tích sử dụng phương pháp chuyển ma trận dải của ma trận độ cứng về một ma trận dạng chữ nhật bằng lệnh system BandGeneral. Số bước lặp lớn nhất 100 bước và độ hội tụ 10-8 được đưa vào để dừng phân tích khi vượt quá số bước lặp mà kết quả chưa hội tụ, sử dụng hàm test NormUnbalance $tol $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 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 quả. Midas/Civil là phần mềm thương mại của hãng Midas IT (Hàn Quốc). Đây là phần mềm phân tích kết cấu có độ chính xác cao cùng với các phần mềm như Sap2000, ADINA và ANSYS. Midas/Civil hiện nay được nhiều kỹ sư và nhà khoa học sử dụng để phân tích kết cấu đặc biệt là trong kết cấu cầu. Các thông số đầu vào đưa vào phần mềm Midas/Civil tương tự như các thông số được sử dụng lập trình trong OpenSees. Hình 12. Mô hình kết cấu trong Midas/Civil KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 18 Tạp chí KHCN Xây dựng - số 4/2015 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 KẾT CẤU – CÔNG NGHỆ XÂY DỰNG Tạp chí KHCN Xây dựng - số 4/2015 19 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 Phần mềm thương mại phải trả phí bản quyền. Phần mềm mã nguồn mở không mất phí Tính linh động Đối với ví dụ trên khi sử dụng chức năng mô hình tiết diện bê tông cốt thép Midas/Civil chỉ cung cấp 7 loại mô hình bê tông và 4 loại mô hình thép. Trong quá trình chia thớ của tiết diện trụ Midas/Civil chỉ cho phép chia 1000 phần tử. Chỉ cho phép liên hợp 6 loại vật liệu trong cùng một tiết diện. Trong khi đó OpenSees cung cấp hơn 200 loại mô hình vật liệu khác nhau và cho phép người dùng tạo ra mô hình vật liệu mới. OpenSees cho phép người sử dụng chia không giới hạn. Cho phép mô hình liên hợp không giới hạn số loại vật liệu trong cùng 1 tiết diện. Một số ưu điểm của OpenSees so với Midas/Civil Phân tích Midas/Civil không cho phép người dùng lựa chọn phương pháp giải phù hợp với từng loại bài toán. Thời gian phân tích trong Midas/Civil trong cùng một ví dụ trên là 301s. OpenSees cho phép người dùng lựa chọn phương pháp giải và độ hội tụ yêu cầu với từng bài toán cụ thể của người sử dụng. Trong ví dụ trên thời gian phân tích trong OpenSees là 61s. Sử dụng Có giao diện, không mất nhiều thời gian để sử dụng. Không có giao diện, tuy nhiên người dùng có thể lập trình tạo giao diện người dùng Một số nhược điểm của OpenSees so với Midas/Civil Tính đơn giản Người dùng không cần biết quá nhiều về kiến thức phần tử hữu hạn cũng có thể giải chính xác được các bài toán phức tạp. OpenSees yêu cầu người dùng phải lựa chọn các phương pháp giải cụ thể do đó người dùng cần phải nắm vững đượ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 Sử dụng OpenSees, việc mô hình, mô phỏng kết cấu rất linh động đặc biệt là các bài toán kết cấu chịu động đất, phân tích phi tuyến và mô hình tương tác 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 đến phương pháp phân tích. Đặc biệt OpenSees là 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 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 người dùng, hiện nay đã có một số chương trình trợ giúp xử lý nhập và xuất dữ liệu nhưng vẫn còn nhiều hạn chế. TÀI LIỆU THAM KHẢO [1]. Website The OpenSees Community Forum: KẾT CẤU – CÔNG NGHỆ XÂY DỰNG 20 Tạp chí KHCN Xây dựng - số 4/2015 [2]. WaiChing Sun (2004): ange/6113-opensees-pre-and-post-processing. [3]. Jinchi Lu, Kevin Mackie, and Ahmed Elgamal. (2011). “OpenSees 3D Pushover and Earthquake Analysis of Single-Column 2-span Bridges”, UC, Berkeley ( [4]. P. Gavali, M S. Shah, G. Kadam và K. Meher. (2013). ”Seismic response and simulations of reinforced concrete bridge using OpenSees on high performance computing”, CSI Transactions on ICT, Volume 1, Issue 3, pp 215-220. [5]. Website: hnjs.htm. [6]. Website: https://www.csiamerica.com/application- programming-interface. [7]. McKenna F, Fenves GL. (2001). The OpenSees command language manual, Version 1.2, Pacific Earthquake Engineering Research Center, University of California at Berkeley. [8]. Website: er/download.php. [9]. Website: _motion_db.html. [10]. B.D.Scott, R.Park, and M.J.Priestley. (1982).. “Stress-strain Behavior of concrete confined by overlapping hoop at low and high strain rates”, ACI Journal, January-February 1982, title no. 79-2. [11]. Finley A. Charney. (April 2008). "Unintended Consequences of Modeling Damping in Structures", J. Struct. Engrg. Volume 134, Issue 4, pp. 581-592. Ngày nhận bài: 05/8/2015. Ngày nhận bài sửa lần cuối: 03/9/2015.