Phần tử MITC3+ được làm trơn trên cạnh dùng phân tích tĩnh tấm Reissner mindlin

Tạp chí Khoa học Công nghệ Xây dựng NUCE 2019. 13 (4V): 139–150<br /> <br /> <br /> <br /> PHẦN TỬ MITC3+ ĐƯỢC LÀM TRƠN TRÊN CẠNH DÙNG<br /> PHÂN TÍCH TĨNH TẤM REISSNER-MINDLIN<br /> <br /> Châu Đình Thànha,∗, Trần Văn Chơnb , Tôn Thất Hoàng Lâna<br /> a<br /> Khoa Xây dựng, Trường Đại học Sư phạm Kỹ thuật Tp. Hồ Chí Minh,<br /> Số 01 đường Võ Văn Ngân, quận Thủ Đức, Tp. Hồ Chí Minh, Việt Nam<br /> b<br /> Công ty Quản lý Dự án Shin Yeong, Tòa nhà SFC,<br /> Tầng 4, 146E Nguyễn Đình Chín, quận Phú Nhuận, Tp. Hồ Chí Minh, Việt Nam<br /> Nhận ngày 09/08/2019, Sửa xong 06/09/2019, Chấp nhận đăng 09/09/2019<br /> <br /> <br /> Tóm tắt<br /> Trong bài báo này, công thức phần tử hữu hạn tấm tam giác 3 nút mới được đề xuất. So với phần tử tấm tam<br /> giác 3 nút truyền thống, xấp xỉ chuyển vị của phần tử đề xuất được bổ sung thêm hàm dạng nổi bậc ba (cubic<br /> bubble shape function) tại nút nổi (bubble node) ở vị trí trọng tâm phần tử. Biến dạng uốn của phần tử được<br /> làm trơn trên trên miền chung cạnh (ES) xác định bởi các đoạn thẳng nối nút nổi của 2 phần tử chung cạnh với<br /> 2 nút của cạnh chung này. Nhờ vào kỹ thuật làm trơn trên cạnh, tích phân trên miền làm trơn của độ cứng uốn<br /> được chuyển sang tích phân trên biên của miền làm trơn và sẽ ít bị ảnh hưởng bởi hình dạng phần tử. Để khử<br /> hiện tượng khóa cắt khi phân tích tấm mỏng, biến dạng cắt ngoài mặt phẳng của phần tử được xấp xỉ lại theo<br /> kỹ thuật khử khóa cắt MITC3+. Phần tử đề xuất, gọi là ES-MITC3+, được sử dụng để phân tích tĩnh một số<br /> bài toán tấm điển hình nhằm đánh giá mức độ chính xác và hội tụ. Thông qua các kết quả số đạt được, phần tử<br /> ES-MITC3+ có khả năng phân tích tĩnh cho cả tấm mỏng và tấm dày với độ chính xác tương đương hoặc tốt<br /> hơn một số loại phần tử khác.<br /> Từ khoá: tấm Reissner-Mindlin; khử khóa cắt MITC3+; phần tử hữu hạn trơn trên cạnh (ES-FEM); phân<br /> tích tĩnh.<br /> AN EDGE-BASED SMOOTHED MITC3+ ELEMENT FOR STATIC ANALYSIS OF REISSNER-MINDLIN<br /> PLATES<br /> Abstract<br /> In this paper, a novel formula of 3-node triangular plate finite element is proposed. In comparison with the<br /> traditional 3-node triangular plate finite elements, the displacements of the proposed element are added the<br /> cubic bubble shape function at the bubble node located at the centroid of the element. The bending strains of<br /> the suggested element are averaged on edge-based smoothed (ES) domains which are determined by straight<br /> lines connecting 2 bubble nodes of 2 adjacent elements with 2 common nodes of them. Thanks to the edge-<br /> based smoothed technique, the integration of the bending stiffness is transformed from the smoothed domain<br /> into its boundary and thus reduces errors due to element shapes. To remove the shear-locking phenomenon,<br /> the transverse shear strains are separately interpolated by using the MITC3+ technique. The proposed element,<br /> namely ES-MITC3+, is employed to statically analyze some benchmark plates for evaluation of the accuracy<br /> and robustness. Numerical results show that the ES-MITC3+ element can analyze both thin and thick plates in<br /> good agreement, or better than other reference elements.<br /> Keywords: Reissner-Mindlin plates; shear-locking removal MITC3+; edge-based smoothed FEM; static analy-<br /> sis.<br /> c 2019 Trường Đại học Xây dựng (NUCE)<br /> https://doi.org/10.31814/stce.nuce2019-13(4V)-13 <br /> <br /> <br /> ∗<br /> Tác giả chính. Địa chỉ e-mail: chdthanh@hcmute.edu.vn (Thành, C. Đ.)<br /> <br /> 139<br />  Thành, C. Đ., và cs. / Tạp chí Khoa học Công nghệ Xây dựng<br /> <br /> 1. Giới thiệu<br /> <br /> Kết cấu tấm là một trong những kết cấu phổ biến được ứng dụng vào nhiều bộ phận công trình xây<br /> dựng như: sàn, mái, tường, vách, . . . do đặc trưng mỏng, nhẹ, khả năng chịu uốn, vượt nhịp lớn. Ứng<br /> xử tấm đồng nhất được tính toán dựa trên (1) lý thuyết tấm mỏng Kirchhoff-Love hoặc (2) lý thuyết<br /> tấm dày Reissner-Mindlin [1]. Ứng xử kết cấu tấm có thể được phân tích bằng các phương pháp giải<br /> tích hoặc các phương pháp số. Các phương pháp giải tích [1–3] cho lời giải có độ chính xác cao, là<br /> kết quả để so sánh đánh giá các nghiên cứu bằng phương pháp số nhưng chỉ áp dụng cho các kết cấu<br /> tấm có hình học, điều kiện biên và tải trọng đơn giản. Vì vậy, để phân tích ứng xử các kết cấu tấm<br /> bất kỳ, các phương pháp số đã và đang được nghiên cứu phát triển nhằm cải thiện độ chính xác, tốc<br /> độ hội tụ và thời gian tính toán. Trong đó, các phương pháp số dựa trên tiếp cận phương pháp phần<br /> tử hữu hạn (PTHH) đang thu hút sự quan tâm của nhiều nhà nghiên cứu trong và ngoài nước. Với giả<br /> thuyết bỏ qua biến dạng cắt ngoài mặt phẳng, lý thuyết Kirchhoff-Love đòi hỏi xấp xỉ chuyển vị của<br /> phương pháp PTHH có đạo hàm cấp 1 liên tục, tức là dạng C 1 . Trong khi đó, xấp xỉ chuyển vị dạng<br /> C 0 đủ để thiết lập công thức PTHH tấm theo lý thuyết tấm dày Reissner-Mindlin. Việc xây dựng hàm<br /> xấp xỉ chuyển vị dạng C 0 , đặc biệt là đối với các phần tử đẳng tham số, dễ hơn rất nhiều so với xây<br /> dựng hàm xấp xỉ chuyển vị dạng C 1 . Tuy nhiên, xấp xỉ chuyển vị dạng C 0 thuần túy sẽ không loại<br /> bỏ được biến dạng cắt ngoài mặt phẳng khi phân tích tấm mỏng và dẫn đến kết quả chuyển vị giảm<br /> khi chiều dày tấm giảm, hay còn gọi là hiện tượng khóa cắt. Cùng với khả năng chia lưới dễ dàng của<br /> phần tử tam giác so với phần tử tứ giác, một số nghiên cứu tập trung vào phát triển phần tử tấm tam<br /> giác 3 nút dựa trên lý thuyết tấm dày Reissner-Mindlin sử dụng xấp xỉ chuyển vị dạng C 0 kết hợp với<br /> các kỹ thuật khử khóa cắt như các phần tử MIN3 [4], DSG3 [5] hoặc MITC3 [6]. Các loại phần tử<br /> này có biến dạng uốn là hằng số trên phần tử và biến dạng cắt ngoài mặt phẳng được xấp xỉ lại để có<br /> thể phân tích ứng xử tấm dày và mỏng. Bằng cách sử dụng thêm hàm nổi bậc 3 ứng với nút nổi đặt tại<br /> trọng tâm phần tử tam giác 3 nút trong xấp xỉ chuyển vị, phần tử MITC3+ [7] có biến dạng uốn tuyến<br /> tính trên phần tử. Do đó, phần tử MITC3+ có độ chính xác và hội tụ tốt hơn phần tử MITC3 [6].<br /> Để giảm chênh lệch biến dạng giữa các phần tử trong kết cấu tấm phân tích bằng các loại phần tử<br /> tam giác 3 nút, các kỹ thuật xấp xỉ lại biến dạng bằng cách trung bình biến dạng giữa các phần tử có<br /> chung cạnh hoặc chung nút, hay còn gọi là phương pháp PTHH trơn, đã được đề xuất [8]. Kỹ thuật<br /> làm trơn trên cạnh (ES) hoặc làm trơn trên nút (NS) đã được áp dụng cho các phần tử DSG3 hoặc<br /> MITC3 để phân tích tấm Reissner-Mindlin [9–12]. Trong các phần tử tấm ES-DSG3 [9], ES-MITC3<br /> [11], NS-DSG3 [10], NS-MITC3 [12], các biến dạng được tính trên miền chung cạnh hoặc chung<br /> nút phần tử là trung bình các biến dạng của các phần tử này. Nhờ đó, biến dạng chênh lệch giữa các<br /> phần tử được làm trơn trên miền nhiều phần tử. Kết quả, các phần tử làm trơn ES-DSG3, ES-MITC3,<br /> NS-DSG3, NS-MITC3 cải thiện được khả năng tính toán so với các phần tử không làm trơn DSG3,<br /> MITC3, DSG3, MITC3.<br /> Trong nghiên cứu này, kỹ thuật làm trơn trên cạnh sẽ được phát triển cho phần tử tấm tam giác 3<br /> nút MITC3+. Với tên gọi ES-MITC3+, phần tử đề xuất có biến dạng uốn được trung bình trên miền<br /> xác định bởi các đoạn thẳng nối nút nổi của 2 phần tử chung cạnh với 2 nút của cạnh chung này, và<br /> biến dạng cắt ngoài mặt phẳng được xấp xỉ theo kỹ thuật khử khóa cắt MITC3+ [7]. Khác với các phần<br /> tử ES-DSG3 và ES-MITC3, biến dạng uốn của phần tử đề xuất không là hằng số trên miền phần tử<br /> nên tích phân trên miền làm trơn của biến dạng uốn được chuyển thành tích phân trên đường biên của<br /> miền làm trơn bằng cách áp dụng định lý phân kỳ Gauss-Ostrogradsky. Nhờ đó, phần tử ES-MITC3+<br /> có thể giảm được sai số do tính tích phân Gauss, đặc biệt trong các trường hợp lưới phần tử không<br /> đều. Trong phần tiếp theo, công thức phần tử ES-MITC3+ được thiết lập. Độ chính xác và tính hiệu<br /> quả của phần tử đề xuất được trình bày ở phần 3 thông qua đánh giá kết quả phân tích chuyển vị và<br /> <br /> 140<br />  Thành, C. Đ., và cs. / Tạp chí Khoa học Công nghệ Xây dựng<br /> <br /> mô-men của một số kết cấu tấm điển hình. Cuối cùng, một số kết luận được tổng kết ở phần 4.<br /> <br /> 2. Công thức phần tử tấm Reissner-Mindlin ES-MITC3+<br /> <br /> 2.1. Công thức phần tử tấm MITC3+<br /> Xét tấm có diện tích mặt trung bình Ω chịu uốn bởi lực pz tác dụng vuông góc với mặt trung bình<br /> Ω như Hình 1. Các chuyển vị thẳng u, v, w tương ứng các phương x, y, z của tấm được xác định bởi [1]<br /> <br /> u (x, y, z) = zβ x (x, y) ; v (x, y, z) = zβy (x, y) ; w (x, y, z) = w0 (x, y) (1)<br /> Tạp chí Khoa học Công nghệ Xây dựng NUCE 2019<br /> trong đó w0 , β x , βyTạp<br /> lầnchí<br /> lượt là chuyển<br /> Khoa học Côngvịnghệ<br /> thẳngXây theo<br /> dựngphương z và góc xoay của pháp tuyến mặt trung<br /> NUCE 2019<br /> bình quanh trục y và trục x với chiều dương qui ước như Hình 1.<br /> <br /> <br /> <br /> <br /> HìnhHình<br /> h 1. Tấm 1. Tấm<br /> chịu chịuvà<br /> 1.uốn<br /> Tấm uốn<br /> chịu và định<br /> định<br /> uốn nghĩa<br /> vànghĩa<br /> định cáccác<br /> các<br /> nghĩa chuyển<br /> chuyểnvịvị<br /> chuyển vị Hình 2.Hình<br /> Hình Phần tử<br /> 2. 2.Phần tấmtửtử<br /> Phần 3tấm<br /> nút3với<br /> tấm 3 nút<br /> nút nútnútvới<br /> với nổi nút<br /> của tấm và của mặt trung bình và chiều dương qui ước<br /> của tấmcủavàtấm<br /> củavàmặt<br /> của mặt trung<br /> trung bình<br /> bình nổi nổi<br /> và chiều dương dương<br /> và chiều qui ước qui ước<br /> Xét<br /> Mặttấm<br /> trung diệnΩtích<br /> có bình củamặt<br /> tấm trung<br /> được rời rạcWbằng<br /> bình chịuN uốn bởi<br /> phần tử lực<br /> tấm ptam<br /> tácgiác<br /> dụng vuông<br /> 3 nút góctích Ω . Trường<br /> có diện<br /> Xét tấm có vị<br /> diện tích tửmặt trung bình W chịu e<br /> uốn z<br /> bởi lực pz tác dụng vuônge góc<br /> với mặt trung bình W như Hình 1. Các chuyển vị thẳng u, v, w tương ứng các phương<br /> chuyển của phần tấm tam giác được xấp xỉ như sau [7]<br /> i mặt trung bình W như Hình 1. Các chuyển vị thẳng u, v, w tương ứng các phương<br /> x, y, z của tấm được xác định bởiX [1]<br /> 3 4 4<br /> y, z của tấm được xác định bởi [1]<br /> X X<br /> w0 = Ni wi ; β x = Ni θyi ; βy = − Ni θ xi (2)<br /> u ( x, y, z ) = zb x ( x, y ) ; v ( x, yi=1<br /> , z ) = zb y ( x, y ) ; i=1<br /> w ( x, y, z ) = w0 ( x, yi=1<br /> ) (1)<br /> u ( x, y, z ) = zb x ( x, y ) ; v ( x, y, z ) = z b y ( x, y ) ; w ( x, y, z ) = w0 ( x, y ) (1)<br /> trong đó,đó<br /> trong , bxθ,xib, yθyilần<br /> w0wi, lầnlượt<br /> lượtlàlàchuyển<br /> chuyển vị thẳng và<br /> vị thẳng theo<br /> gócphương<br /> xoay quanhz và trục<br /> góc xxoay củaypháp<br /> và trục của nút i với chiều<br /> ong đó,tuyến bx, qui<br /> wdương<br /> 0, mặt by ước<br /> lần<br /> trung lượt<br /> được<br /> bình là nghĩa<br /> định<br /> quanh chuyển<br /> trục ytrong vị xthẳng<br /> Hình<br /> và trục 2. N<br /> với i theo<br /> là các<br /> chiều phương<br /> hàm<br /> dương dạng<br /> qui ước z như<br /> vàhệHình<br /> trong góc xoay<br /> tọa độ<br /> 1. tự nhiên η) được<br /> của(ξ,pháp<br /> xác định ứng với các nút đỉnh (i = 1, 2, 3) và nút nổi (i = 4) đặt tại trọng tâm phần tử như sau<br /> yến mặt trung<br /> Mặtbình<br /> trung quanh<br /> bình W trục y và<br /> của tấm trục<br /> được rờixrạc<br /> vớibằng<br /> chiều dương<br /> Ne phần quitam<br /> tử tấm ướcgiác<br /> như Hình<br /> 3 nút có 1.<br /> N1 = 1 − ξ − η − 9ξη (1 − ξ − η) ; N2 = ξ − 9ξη (1 − ξ − η)<br /> Mặtdiện tích W<br /> trung bình W của<br /> e. Trường chuyển<br /> tấmvịđược<br /> của phần<br /> rờitửrạc<br /> tấmbằng<br /> tam giác<br /> Neđược<br /> phần xấptửxỉtấm<br /> như sau<br /> tam[7]giác 3 nút có(3)<br /> N3 = η − 9ξη (1 − ξ − η) ; N4 = 27ξη (1 − ξ − η)<br /> ện tích We. Trường 3<br /> chuyển 4vị của phần tử 4<br /> tấm tam giác được xấp xỉ như sau (2) [7]<br /> w0 = å Ni wi ; b x = å Niq yi ; b y = -å Niq xi<br /> 3<br /> Từ (1)<br /> i =1<br /> và (2), quan<br /> 4<br /> hệ giữa biến dạng<br /> i =1 4i =1<br /> và chuyển vị nút phần tử được thiết lập<br /> w0 = å Ni wi ; b x = å Niqyi ; b y = -å  Niq xi  (2)<br /> ở đây, wi, qxi, q ε xi =lượt<br /> yi lần  là chuyểnβ x,xvị thẳng<br /> 1 vàX 4góc<br />  0xoay 0quanhNi,x trục x và<br /> wi trục y của<br /> 4 nút<br />  <br /> i =1 1 i=<br />    <br />   <br />  X<br /> εy ước = z được βy,y nghĩa = trong θ xi  =z<br />  <br /> z 0 −N2. i,yNi là 0các hàm Bhệ<br /> bi dei (4)<br />        <br /> i với chiều dương  γqui<br /> <br />  β định<br />   Hình<br />  θ dạng<br /> <br />  trong<br /> , qxiđộ, qtự<br />     <br /> đây, witọa yi lần lượt xy là chuyển + β<br /> x,y vị y,xthẳng và góc xoay i,x quanh trụcyi x và trục y của nút<br />      <br />   i=1 0 −N N i=1<br /> nhiên (x,h) được| xác {z<br /> i,y<br /> định ứng } với các | nút đỉnh<br /> {z (i =}1,| 2, {z3)} và nút nổi<br /> ới chiều<br /> (i =dương<br /> 4) đặt tạiqui<br /> trọng ước tâmđượcphần tửđịnh<br /> như nghĩa<br /> κ<br /> sau trong HìnhBbi2. Ni là các dei hàm dạng trong hệ<br /> <br /> a độ tự nhiên (x,h) được xác định ứng với các nút đỉnh (i = 1, 2, 3) và nút nổi<br /> N1 = 1 - x - h - 9xh (1 - x - h ) ; N 2 = x - 9xh (1 - x - h ) ; N 3 = h - 9xh (1 - x - h )<br /> = 4) đặt tại Ntrọng tâm phần tử như sau 141 (3)<br /> (<br /> = 27xh 1 - x - h<br /> 4 )<br /> N1 = 1 - x - h - 9xh (1 - x - h ) ; N 2 = x - 9xh (1 - x - h ) ; N 3 = h - 9xh (1 - x - h )<br /> Từ (1) và (2), quan hệ giữa biến dạng và chuyển vị nút phần tử được thiết lập (3)<br /> N 4 = 27xh (1 - x - h )<br />  Thành, C. Đ., và cs. / Tạp chí Khoa học Công nghệ Xây dựng<br /> <br /> 4 " # wi <br /> <br /> 4<br /> γ xz β x + w0,x<br /> ( ) ( ) <br /> X Ni,x 0 Ni   <br />  X<br /> = = θ =<br />  <br /> xi  B si dei (5)<br /> γyz βy + w0,y Ni,x −Ni 0 <br />  <br /> } θyi<br />  <br /> i=1 | i=1<br />  <br /> {z<br /> B si | {z }<br /> dei<br /> <br /> Do xấp xỉ (2) không sử dụng độ võng tại nút nổi (i = 4) nên trong các công thức (4) và (5) w4 = 0.<br /> Trong bài báo này ký hiệu ,  là đạo hàm của hàm  đối với biến .<br /> Vì biến dạng cắt ngoài mặt phẳng xấp xỉ bởi (5) không thể tiến đến 0, nghĩa là khi phân tích tấm<br /> mỏng mỏng biến dạng cắt ngoài mặt phẳng vẫn tồn tại. Điều này không phù hợp với ứng xử thực tế<br /> của tấm mỏng và dẫn đến kết quả, tấm càng mỏng thì độ võng càng lớn. Đây chính là hiện tượng khóa<br /> cắt xảy ra với các phần tử dùng hàm xấp xỉ chuyển vị bậc thấp. Do đó, biến dạng cắt ngoài mặt phẳng<br /> được xấp xỉ lại theo kỹ thuật nội suy các thành phần ten-xơ hỗn hợp MITC3+ [7]. Cụ thể, biến dạng<br /> cắt ngoài mặt phẳng trong hệ tọa độ tự nhiên được xấp xỉ lại như sau<br /> !<br /> 2 B 1 B 1 C  1<br /> γˆ ξζ = γξζ − γηζ + γξζ + γCηζ + cˆ (3η − 1)<br /> 3 2 3 3<br /> ! (6)<br /> 2 A 1 A 1 C  1<br /> γˆ ηζ = γ − γ + γ + γηζ + cˆ (3ξ − 1)<br /> C<br /> 3 ξζ 2 ηζ 3 ξζ 3<br />    <br /> với cˆ = γξζ<br /> F<br /> − γξζ<br /> D<br /> − γηζ<br /> F<br /> + γηζ<br /> E<br /> . Ở đây, γξζ<br /> I<br /> , γηζ<br /> I<br /> là các giá trị của biến dạng cắt ngoài mặt phẳng tính<br /> tại các điểm buộc I = A, B, C, D, E, F có tọa độ cho trong Bảng 1.<br /> Bảng 1. Tọa độ các điểm buộc của kỹ thuật khử khóa cắt MITC3+ với d = 1/10000<br /> <br /> Điểm buộc ξ η<br /> A 1/6 2/3<br /> B 2/3 1/6<br /> C 1/6 1/6<br /> D 1/3 + d 1/3 − 2d<br /> E 1/3 − 2d 1/3 + d<br /> F 1/3 + d 1/3 + d<br /> <br /> Từ các tọa độ điểm buộc cho trong Bảng 1, các giá trị biến dạng cắt ngoài mặt phẳng trong hệ tọa<br /> độ (x, y) được xác định theo (5). Sử dụng công thức biến đổi biến dạng cắt ngoài mặt phẳng từ hệ tọa<br /> độ (x, y) sang hệ tọa độ (ξ, η) và thế vào xấp xỉ biến dạng cắt cho bởi (6), ta có thể thiết lập được quan<br /> hệ giữa biến dạng cắt ngoài mặt phẳng và chuyển vị nút phần tử theo kỹ thuật khử khóa cắt MITC3+<br /> như sau<br /> 4<br /> γˆ xz<br /> ( ) X<br /> = ˆ si dei<br /> B (7)<br /> γˆ yz<br /> i=1<br /> Quan hệ giữa ứng suất và biến dạng của tấm đồng nhất đẳng hướng xác định<br /> σx   1 ν εx   1 ν<br />       <br /> 0 0 4<br /> E  Ez<br /> <br />       X<br /> σy  =  ν 1 εy  =  ν 1<br />    <br /> 0 0 Bbi dei (8)<br />   <br />    <br />  <br />  1 − ν2  0 0 (1 − ν)/2    1 − ν2  0 0 (1 − ν)/2<br /> <br />    <br />  τ <br />    γ <br />  <br /> i=1<br /> xy xy<br /> 4<br /> τ xz γˆ xz<br /> ( ) ( )<br /> E E X<br /> = = Bˆ si dei (9)<br /> τyz 2(1 − ν) γˆ yz 2(1 − ν) i=1<br /> <br /> 142<br />  Thành, C. Đ., và cs. / Tạp chí Khoa học Công nghệ Xây dựng<br /> <br /> trong đó E là mô-đun đàn hồi và ν là hệ số Poisson của vật liệu.<br /> Nguyên lý công ảo của tấm có diện tích mặt trung bình Ω chịu tải trọng pz tác dụng vuông góc<br /> với mặt trung bình [13]<br /> <br /> Z Zt/2 h σx  Z Zt/2 h<br />  <br /> i<br />   kt2 i ( τ xz ) Z<br /> δε x δεy δγ xy  σy  dzdΩ + δˆγ xz δˆγyz dzdΩ = δwpz dΩ<br />  <br /> (10)<br />  <br /> <br />  τ <br />  <br />  t2 + αh2e τyz<br /> Ω −t/2 xy Ω −t/2 Ω<br /> <br /> trong đó δ chỉ đại lượng ảo, k = 5/6 là hằng số hiệu chỉnh kể đến sự phân bố không đều theo chiều<br /> dày của các ứng suất cắt ngoài mặt phẳng τ xz , τyz , he là cạnh dài nhất của phần tử và α = 0,1 là hệ số<br /> ổn định [14].<br /> Thế các quan hệ giữa ứng suất – biến dạng (8), (9) vào (10), ta có<br />    <br /> XNe Z  XNe Z  XNe Z<br /> δdTe  BTb Db Bb dΩ de + δdTe  B  d = δd<br />    ˆT ˆ T<br /> D B<br /> s s s dΩ  e e Npz dΩ (11)<br /> e=1<br />   e=1<br />  e=1<br /> Ωe Ωe Ωe<br /> <br /> <br /> trong đó Bb = [Bb1 Bb2 Bb3 Bb4 ], B s = [B s1 B s2 B s3 B s4 ], de = [dTe1 dTe2 dTe3 dTe4 ]T , N =<br /> [N1 0 0 N2 0 0 N3 0 0 N4 0 0]T và<br />  <br />  1 v 0<br /> Et3<br /> <br /> Db = D  v 1 0 với D = (12)<br />  <br /> 12 1 − ν2<br /> 
<br /> 0 0 (1 − ν)/2<br />  <br /> <br /> <br /> kEt3<br /> " #<br /> 1 0<br /> Ds =  (13)<br /> 0 1<br /> <br /> t2 + αh2e 2(1 + v)<br /> <br /> Từ (11), phương trình cân bằng rời rạc PTHH của tấm chịu tải trọng tĩnh pz được viết dưới dạng<br /> <br /> Kd = F (14)<br /> <br /> trong đó d, K, F lần lượt là véc-tơ chuyển vị, ma trận độ cứng và véc-tơ tải trọng của tấm. Ma trận K<br /> và véc-tơ F tương ứng được lắp ghép từ các ma trận độ cứng ke và véc-tơ tải trọng fe của phần tử có<br /> công thức như sau<br /> Z Z<br /> ke = Bb Db Bb dΩ + B<br /> T ˆ Ts D s B<br /> ˆ s dΩ (15)<br /> Ωe Ωe<br /> Z<br /> fe = Npz dΩ (16)<br /> Ωe<br /> <br /> <br /> 2.2. Công thức phần tử tấm ES-MITC3+<br /> Trong nghiên cứu này, biến dạng uốn sẽ được trung bình trên miền làm trơn có diện tích Ωk được<br /> giới hạn bởi các đoạn thẳng nối 2 nút nổi của 2 phần tử chung cạnh với 2 nút đỉnh của cạnh chung này<br /> như Hình 3.<br /> <br /> <br /> <br /> 143<br />  Tạp chí Khoa học Công nghệ Xây dựng NUCE 2019<br /> <br /> <br /> <br /> (<br /> fe = ò Np dWz<br /> Thành, C. Đ., và cs. / Tạp chíW Khoa học Công nghệ Xây dựng<br /> e<br /> <br /> <br /> <br /> <br /> Do đó, biến dạng uốn được xấp xỉ lại2.2.<br /> như Công<br /> sauthức phần tử tấm ES-MITC3+<br />  Trong nghiên cứunày, biến dạng uốn sẽ được trung bình trên miền làm trơn<br /> ε εx <br /> <br /> ˜<br />  tích 1Wk được<br /> x<br />  diện  giới hạn<br />  bởi các đoạn thẳng nối 2 nút nổi của 2 phần tử chung cạ<br /> <br />  <br />  Z   <br /> ε = ε<br />    <br /> ˜ y<br /> với 2 nút đỉnh của y<br /> cạnh<br /> dΩ (17)<br />  γ˜  Ωk chung này như Hình 3.<br /> <br />  <br />  <br />  <br />    γ <br />  <br /> xy Ωk xy<br /> Do đó, biến dạng uốn được xấp xỉ lại như sau<br /> Thế (4) vào (17) ta được quan hệ giữa biến ì e!dạng<br /> x<br /> ü trơnìvà e x chuyển<br /> ü vị nút phần tử<br /> ï ï 1 ï ï (<br /> í e! y ý = Z ò í e y ý dW<br /> ïg! ï W k W ïg ï<br />  <br /> k<br />  0<br />  î<br /> 0 xy þ î<br /> N þ<br /> xy dΩ<br /> i,x<br /> <br /> <br /> Thế (4) vào (17)<br /> Ω ta được quan  hệ giữa biến dạng trơn và chuyển vị nút phần tử<br />  <br /> k<br /> ε<br />   <br /> ˜ wi <br />  Z<br /> x  4 4<br /> 1  0 − Ni,y dΩ<br />    <br /> 0é =Nz dW ùúB<br />  X    X<br /> ε˜ y  = z 0 θ xi <br /> ˜ bi dei<br />     <br /> (18)<br />    <br /> Ω ê0 ò<br />  <br />   <br />  i,x<br />  γ˜ <br />  <br />  i=1 k  Ω<br /> Zk ì e"x ü ê  <br />  θyi W<br /> <br />  <br /> <br />  i=1ú ì wi ü<br /> xy  k<br /> <br /> ï ï 4<br /> 1 ê úï ï 4 (<br /> Z <br /> ý = z å Ni,y<br /> e" ydΩ - òN i , y dW q xi ý = z å B " d<br />  | {z }<br />  0 − Níi,x ê0dΩ 0 ú í<br /> <br /> bi ei<br /> i =1 W k ê<br /> ïg" ï W k<br />  dei ú ïq ï i =1<br /> î xy þ<br /> ê ú!î yi þ<br /> <br /> Ωk Ω<br /> {z k ê0 - } ò Ni , xdW ò Ni , y dW ú d ei<br /> <br /> ë &û<br /> |<br /> Wk Wk<br /> ˜ bi<br /> B #$$$$ $%$$$$$<br /> "<br /> B bi<br /> <br /> <br /> <br /> Áp dụng định lý phân kỳ Gauss-Ostrogradsky,<br /> tích phân trên miền làm trơn Ωk của các đạo hàm<br /> hàm dạng trong (18) sẽ được chuyển thành tích<br /> phân trên đường biên Γk của Ωk như sau<br /> Z Z Z Z<br /> Ni,x dΩ = Ni n x dΓ; Ni,y dΩ = Ni ny dΓ<br /> Ωk Γk Ωk Γk<br /> (19)<br /> trong đó n x và ny lần lượt là hình chiếu theo<br /> phương x và y của véc-tơ n pháp tuyếnHình với 3. Miền làm trơn trên cạnh (ES) cho biến dạng uốn của các phần tử ES-MITC<br /> biên<br /> Hình 3. Miền làm trơn trên cạnh (ES) cho biến<br /> Γk như biểu diễn ở Hình 3. và định nghĩa véc-tơ pháp tuyến của biên miền làm trơn<br /> dạng uốn của các phần tử ES-MITC3+ và định<br /> Với các hàm dạng cho bởi (3), tích phân đường<br /> Áp dụng địnhnghĩalý phân kỳ Gauss-Ostrogradsky,<br /> véc-tơ pháp tuyến của biên tích<br /> miềnphân<br /> làmtrên<br /> trơnmiền làm trơn<br /> ở (19) được tính chính xác bằng cách sử dụng 2<br /> của các đạo hàm hàm dạng trong (18) sẽ được chuyển thành tích phân trên đường b<br /> điểm Gauss trên mỗi đoạn thẳng của biên Γk . WDo<br /> Gk của k như sau<br /> đó, tích phân ở (19) có thể được tính như sau<br /> òN i,x dW = ò N n dG; ò N<br /> i x i, y dW = ò N n dG<br /> i y<br /> (<br /> Wk Gk Wk Gk<br /> Z Z Ned X<br /> X 2  <br /> Ni,x dΩ = n dΓ<br /> Ntrong = i ξgp<br /> i x đó, nx và ny lầnNlượt là, η<br /> ed ed<br /> hình wed ed<br /> gp n xtheo phương x và y của véc-tơ n pháp tuy<br /> gp chiếu<br /> Ωk Γk ed=1 gp=1<br /> Z Z Ned X<br /> 2 7 (20)<br /> X  <br /> Ni,y dΩ = Ni ny dΓ = Ni ξgp , ηgp<br /> ed ed<br /> wed ed<br /> gp ny<br /> Ωk Γk ed=1 gp=1<br /> <br /> <br /> trong đó Ned = 3 đối với miền Ωk nằm ở biên tấm và Ned = 4 đối với miền Ωk không nằm ở biên tấm.<br /> ξgp , ηgp , wed<br /> ed ed<br /> gp , n x , ny lần lượt là tọa độ tự nhiên, trọng số của điểm Gauss và thành phần của véc-tơ<br /> ed ed<br /> <br /> pháp tuyến n trên cạnh ed của biên Γk .<br /> <br /> <br /> <br /> <br /> 144<br />  Thành, C. Đ., và cs. / Tạp chí Khoa học Công nghệ Xây dựng<br /> <br /> Vì vậy, biến dạng uốn sau khi được làm trơn có thể được xác định bởi<br /> <br />  XNed X 2  <br /> <br /> <br />  0<br />  0 Ni ξ ,<br /> gp gpη<br /> ed ed<br /> w ed ed <br /> n<br /> gp x <br /> ed=1 gp=1<br /> <br /> ε<br />   <br />    <br />  ˜<br />  x <br />   4<br /> X 1   XNed X 2      wi <br /> ε = ξ , η<br /> ed ed ed ed   θ xi <br />   <br /> <br /> ˜ y<br /> <br /> z  0 − Ni w n 0  <br /> gp gp gp y<br /> Ωk <br /> <br />  <br />  <br />  γ˜ <br />  <br /> i=1<br /> <br /> <br /> ed=1 gp=1<br /> <br />  <br />   θ <br />  <br /> <br /> xy  yi<br /> N N<br />  2 2<br /> <br /> ed ed<br />  X X  X X  <br /> <br />  | {z }<br />  <br />  0 −<br />  Ni ξgp , η<br /> ed ed<br /> gp w ed ed<br /> gp n x Ni ξgp , η<br /> ed ed<br /> gp w ed ed <br /> gp ny  dei<br /> ed=1 gp=1 ed=1 gp=1<br /> | {z }<br /> ˜ bi<br /> B<br /> 4<br /> X<br /> =z ˜ bi dei<br /> B<br /> i=1<br /> (21)<br /> Tạptrơn<br /> Thế biến dạng uốn được xấp xỉ lại theo kỹ thuật làm chí Khoa<br /> trên học<br /> cạnhCông<br /> chonghệ<br /> bởi Xây<br /> (21)dựng NUCE 2019<br /> vào nguyên lý<br /> công ảo (10), phương trình cân bằng rời rạc PTHH của tấm chịu uốn được viết lại<br /> <br /> k! e = Β ˜ != F<br /> ! T Kd<br /> b Db B b W k + ò Βˆ ˆ dW (22)<br /> T<br /> s Ds B s<br /> We<br /> trong đó K ˜ được lắp ghép từ các ma trận độ cứng phần tử ES-MITC3+ có biến dạng uốn được làm<br /> trơn trên cạnh 3. VíZ dụ số Z<br /> e = Bài b dΩ +test B<br /> ˜k3.1. ˜ T ˜ ˆ Ts D s B<br /> ˆ s dΩ<br /> Btoán<br /> b Db Bpatch (23)<br /> Ωe Ωe<br /> h i Để kiểm tra khả năng xấp xỉ trường biến dạng và ứng suất của ph<br /> ˜b = B<br /> trong đó B ˜ b1 B˜ b2 B<br /> ˜ b3 B˜ b4<br /> MITC3+,<br /> . xét tấm chữ nhật dày t = 0,01 m có tọa độ nút và lưới phần tử n<br /> Theo cách tính B˜ bi cho bởi (21)<br /> Tấm làm bằngB˜ bvật<br /> thì B˜ bi hay là hằng<br /> liệu cósố.mô-đun<br /> Do đó, công<br /> đàn thức<br /> hồi Etính ma7 trận<br /> = 10 kN/m độvà<br /> cứng<br /> hệ số Poisson<br /> phần tử ES-MITC3+ có thể được viết lại 2 2<br /> Tấm chịu độ võng cưỡng bức w = (1 + x + 2y + x + xy + y ) / 200 m [9].<br /> Z<br /> ˜ke = Theo<br /> ˜BT Db Blý thuyết<br /> ˜ bΩ k + tấm ˆ Ts Dmỏng,<br />