Luận án Phân tích tĩnh và động tấm, vỏ thoải hai độ cong composite nano carbon - Áp điện

BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC XÂY DỰNG 
Vũ Văn Thẩm 
PHÂN TÍCH TĨNH VÀ ĐỘNG TẤM, VỎ THOẢI HAI ĐỘ 
CONG COMPOSITE NANO CARBON - ÁP ĐIỆN 
Chuyên ngành: Cơ kỹ thuật 
Mã số: 9520101 
LUẬN ÁN TIẾN SĨ 
Hà Nội - Năm 2021 
V
Ũ
 V
Ă
N
 T
H
Ẩ
M
 *
 L
U
Ậ
N
 Á
N
 T
IẾ
N
 S
Ĩ *
 M
Ã
 S
Ố
 9
5
2
0
1
0
1
 *
 N
Ă
M
 2
0
2
1
BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC XÂY DỰNG 
Vũ Văn Thẩm 
PHÂN TÍCH TĨNH VÀ ĐỘNG TẤM, VỎ THOẢI HAI ĐỘ 
CONG COMPOSITE NANO CARBON - ÁP ĐIỆN 
Chuyên ngành: Cơ kỹ thuật 
Mã số: 9520101 
 NGƯỜI HƯỚNG DẪN KHOA HỌC 
 1. PGS. TS Trần Hữu Quốc 
 2. GS. TS Trần Minh Tú 
Hà Nội - Năm 2021
i 
LỜI CAM ĐOAN 
Tên tôi là: Vũ Văn Thẩm 
Tôi xin cam đoan đây là công trình nghiên cứu của riêng cá nhân tôi. Các số 
liệu và kết quả được trình bày trong luận án là trung thực, đáng tin cậy và không trùng 
lặp với bất kỳ nghiên cứu nào khác đã thực hiện. 
 Hà Nội, ngày 05 tháng 06 năm 2020 
Tác giả 
Vũ Văn Thẩm 
ii 
LỜI CẢM ƠN 
Tác giả xin bày tỏ lòng biết ơn sâu sắc đến hai thầy giáo hướng dẫn là PGS. TS 
Trần Hữu Quốc và GS. TS Trần Minh Tú đã tận tình hướng dẫn, giúp đỡ, động viên 
trong suốt quá trình học tập, nghiên cứu và hoàn thành luận án. 
Tác giả chân thành cảm ơn tập thể các thầy cô giáo, các bạn đồng nghiệp Bộ 
môn Sức bền Vật liệu, Trường Đại học Xây dựng - Nơi tác giả đang công tác đã luôn 
quan tâm, giúp đỡ và tạo các điều kiện thuận lợi nhất để tác giả có thể hoàn thành tốt 
nhiệm vụ được giao và học tập, nghiên cứu, hoàn thành luận án. 
Tác giả xin cảm ơn tập thể các thầy cô giáo, cán bộ Khoa Đào tạo Sau đại học, 
Trường Đại học Xây dựng đã tạo điều kiện thuận lợi, giúp đỡ trong quá trình thực 
hiện luận án. 
Tác giả xin bày tỏ lòng biết ơn các bạn bè, đồng nghiệp đã tận tình giúp đỡ và 
động viên tác giả học tập, nghiên cứu làm luận án. Cuối cùng tác giả xin bày tỏ lòng 
biết ơn sâu sắc đến các thành viên trong gia đình đã luôn tạo điều kiện, chia sẻ những 
khó khăn trong suốt quá trình học tập, nghiên cứu và hoàn thành luận án. 
Tác giả: Vũ Văn Thẩm 
iii 
MỤC LỤC 
Nội dung Trang 
LỜI CAM ĐOAN ....................................................................................................... i 
LỜI CẢM ƠN ............................................................................................................ ii 
MỤC LỤC ........................................................................................................... iii 
DANH MỤC CÁC KÝ HIỆU.................................................................................. vii 
DANH MỤC CÁC CHỮ VIẾT TẮT ....................................................................... ix 
DANH MỤC BẢNG BIỂU ........................................................................................ x 
DANH MỤC HÌNH VẼ ........................................................................................... xii 
MỞ ĐẦU ............................................................................................................. 1 
1. Lý do lựa chọn đề tài ........................................................................................ 1 
2. Mục tiêu nghiên cứu ......................................................................................... 2 
3. Đối tượng, phạm vi nghiên cứu ........................................................................ 3 
4. Cơ sở khoa học của luận án .............................................................................. 3 
5. Nội dung nghiên cứu của luận án ..................................................................... 3 
6. Phương pháp nghiên cứu .................................................................................. 4 
7. Ý nghĩa khoa học và thực tiễn của luận án ....................................................... 4 
8. Những đóng góp mới của luận án..................................................................... 5 
9. Cấu trúc của luận án ......................................................................................... 5 
CHƯƠNG 1 TỔNG QUAN VẤN ĐỀ NGHIÊN CỨU ........................................... 6 
1.1. Vật liệu composite ....................................................................................... 6 
1.2. Ống nano carbon .......................................................................................... 8 
1.3. Vật liệu áp điện ............................................................................................ 9 
1.4. Vật liệu composite nano carbon - áp điện ................................................. 10 
1.5. Tổng quan các nghiên cứu liên quan đến nội dung đề tài luận án ............. 11 
1.5.1. Các nghiên cứu về kết cấu tấm, vỏ composite nano carbon .............. 11 
iv 
 tấm, vỏ composite nano carbon – áp điện ................................................... 16 
1.5.3 Các mô hình tính toán kết cấu tấm, vỏ composite áp điện ................. 16 
1.5.4 Các nghiên cứu về tấm, vỏ composite áp điện ở Việt Nam ................ 20 
1.6. Nhận xét chương 1 ..................................................................................... 21 
CHƯƠNG 2 PHÂN TÍCH TĨNH VÀ ĐỘNG TẤM, VỎ THOẢI HAI ĐỘ CONG 
COMPOSITE NANO CARBON – ÁP ĐIỆN THEO TIẾP CẬN GIẢI TÍCH ........ 22 
2.1 Mở đầu ....................................................................................................... 22 
2.2 Mô hình vỏ thoải hai độ cong composite nano carbon – áp điện .............. 22 
2.2.1. Dạng hình học của vỏ thoải hai độ cong PFG-CNTRC .................... 22 
2.2.2. Cơ học vật liệu PFG-CNTRC ........................................................... 23 
2.3 Phân tích vỏ thoải hai độ cong PFG-CNTRC theo HSDST-4 ................... 25 
2.3.1 Các giả thiết....................................................................................... 25 
2.3.2 Các thành phần chuyển vị cơ học ..................................................... 26 
2.3.3 Các thành phần biến dạng cơ học ..................................................... 27 
2.3.4 Chuyển dịch điện tích trong lớp áp điện ........................................... 30 
2.3.5 Trường ứng suất ................................................................................ 32 
2.3.6 Phương trình chuyển động ................................................................ 34 
2.3.7 Điều kiện biên ................................................................................... 42 
2.4 Lời giải giải tích ......................................................................................... 43 
2.5 Nhận xét chương 2 ..................................................................................... 48 
CHƯƠNG 3 PHÂN TÍCH TĨNH VÀ ĐỘNG TẤM, VỎ THOẢI HAI ĐỘ CONG 
COMPOSITE NANO CARBON – ÁP ĐIỆN BẰNG PHƯƠNG PHÁP PHẦN TỬ 
HỮU HẠN ............................................................................................................. 50 
3.1 Mở đầu ....................................................................................................... 50 
3.2 Lựa chọn mô hình phần tử ......................................................................... 50 
3.3 Các phương trình cơ bản ............................................................................ 52 
3.3.1. Trường chuyển vị .............................................................................. 52 
v 
3.3.2. Trường biến dạng .............................................................................. 53 
3.3.3. Trường ứng suất ................................................................................ 54 
3.4 Mô hình phần tử hữu hạn ........................................................................... 55 
3.4.1. Các hàm nội suy ................................................................................ 55 
3.4.2. Các liên hệ tọa độ .............................................................................. 59 
3.4.3. Phương trình chuyển động ................................................................ 63 
3.5 Nhận xét chương 3 ..................................................................................... 70 
CHƯƠNG 4 KHẢO SÁT ẢNH HƯỞNG CỦA MỘT SỐ YẾU TỐ ĐẾN ĐẶC 
TRƯNG TĨNH VÀ ĐỘNG CỦA TẤM, VỎ THOẢI HAI ĐỘ CONG COMPOSITE 
NANO CARBON – ÁP ĐIỆN .................................................................................. 72 
4.1 Mở đầu ....................................................................................................... 72 
4.2 Các ví dụ kiểm chứng mô hình HSDST-4 ................................................. 73 
4.2.1. Kiểm chứng bài toán dao động riêng ................................................ 73 
4.2.2. Kiểm chứng bài toán uốn .................................................................. 80 
4.2.3. Nhận xét các ví dụ kiểm chứng ......................................................... 81 
4.3 Khảo sát bài toán uốn ................................................................................. 82 
4.3.1 Độ võng và ứng suất của kết cấu tấm, vỏ thoải composite PFG-
CNTRC. ......................................................................................................... 83 
4.3.2 Ảnh hưởng của điện thế áp đặt đến độ võng ..................................... 92 
4.3.3 Ảnh hưởng của điện thế áp đặt đến sự phân bố ứng suất ................. 93 
4.3.4 Ảnh hưởng của điều kiện biên .......................................................... 96 
4.4 Khảo sát bài toán dao động riêng............................................................. 100 
4.4.1 Tần số dao động riêng của vỏ thoải composite PFG-CNTRC ........ 100 
4.4.2 Ảnh hưởng của kiểu phân bố CNT ................................................. 107 
4.4.3 Ảnh hưởng của trạng thái mạch ...................................................... 107 
4.4.4 Ảnh hưởng của độ thoải vỏ PFG-CNTRC ...................................... 108 
4.4.5 Ảnh hưởng của số lớp vật liệu FG-CNTRC ................................... 109 
vi 
4.4.6 Ảnh hưởng của chiều dày lớp áp điện ............................................. 109 
4.4.7 Ảnh hưởng của điều kiện biên ........................................................ 110 
4.5 Bài toán đáp ứng chuyển vị theo thời gian của vỏ thoải composite PFG-
CNTRC ............................................................................................................. 113 
4.5.1. Ảnh hưởng của dạng tải trọng cưỡng bức tác dụng theo thời gian . 113 
4.5.2. Ảnh hưởng của V*CNT ...................................................................... 115 
4.5.3. Ảnh hưởng của kiểu phân bố CNT ................................................. 116 
4.6 Nhận xét chương 4 ................................................................................... 118 
KẾT LUẬN ......................................................................................................... 119 
KIẾN NGHỊ ......................................................................................................... 120 
DANH MỤC CÁC CÔNG TRÌNH ĐÃ CÔNG BỐ............................................... 121 
TÀI LIỆU THAM KHẢO ....................................................................................... 122 
PHỤ LỤC ........................................................................................................ PL1 
vii 
DANH MỤC CÁC KÝ HIỆU 
a, b Kích thước các cạnh hình chiếu bằng của vỏ lần lượt theo 
các phương x, y 
hc Chiều dày của lớp lõi composite 
hp Chiều dày của lớp áp điện 
Rx, Ry Bán kính cong theo các phương trục x, y của vỏ 
k Số thứ tự lớp 
hk Chiều dày lớp thứ k 
 1,2,3 Hệ toạ độ vật liệu 
 , ,x y z Hệ toạ độ toàn cục 
 , ,l l lx y z Hệ toạ độ phần tử 
 ,zq x y
 Tải trọng cơ học phân bố tác dụng lên mặt trên của vỏ 
 , ; ,t bq x y q x y Tải trọng điện phân bố tác dụng lên mặt trên của vỏ 
, , , b su v w w Chuyển vị theo các phương x, y, z 
  Véc tơ các thành phần biến dạng 
  Véc tơ các thành phần ứng suất 
[ ]ijQ Ma trận độ cứng vật liệu composite trong hệ (1,2,3) 
[ ]ijQ Ma trận độ cứng vật liệu composite trong hệ (x,y,z) 
[ ]ijC Ma trận độ cứng vật liệu áp điện trong hệ (1,2,3) 
[ ]ijC Ma trận độ cứng vật liệu áp điện trong hệ (x,y,z) 
 kD Véc tơ chuyển dịch điện tích trong lớp áp điện thứ k 
viii 
 kE Véc tơ cường độ điện trường của lớp áp điện thứ k 
[ ]ije Ma trận các hệ số ứng suất áp điện 
[ ]ijp Ma trận các hệ số điện môi 
(x, y,z)k Điện thế tại điểm bất kỳ trong lớp áp điện thứ k 
(x, y,z)k Điện thế tại mặt trung bình của lớp áp điện thứ k 
 Biến phân của các đại lượng 
   ; ; M N Q Véc tơ các thành phần nội lực màng, mô men, lực cắt 
U; W; T Thế năng biến dạng đàn hồi; Công của ngoại lực; Động năng 
[ ]uuK Ma trận độ cứng cơ học của vỏ 
[ ]uK  Ma trận độ cứng tương tác cơ – điện 
[ ]uK Ma trận độ cứng tương tác điện – cơ 
[ ]K Ma trận độ cứng điện 
[ ]M Ma trận khối lượng 
 F Véc tơ lực nút 
 d Véc tơ các thành phần chuyển vị 
  Véc tơ các chuyển dịch điện tích 
; Gd vG Hệ số hồi tiếp chuyển vị và hệ số hồi tiếp vận tốc 
ix 
DANH MỤC CÁC CHỮ VIẾT TẮT 
CNT Carbon nanotube (ống nano carbon) 
FG-CNTRC 
Functionally graded carbon nanotube reinforced composite (vật 
liệu composite có cơ tính biến thiên được gia cường bởi ống nano 
carbon) 
PFG-CNTRC Vật liệu composite áp điện được gia cường bởi ống nano carbon 
FGM 
Functionally Graded Material (vật liệu có cơ tính biến thiên hay 
vật liệu biến đổi chức năng) 
3D Three-dimensional elasticity theory (lý thuyết đàn hồi ba chiều) 
CST Classical shell theory (lý thuyết vỏ cổ điển) 
ESL Lý thuyết đơn lớp tương đương 
FSDT First-order shear deformation theory (lý thuyết biến dạng cắt bậc 
nhất) 
HSDT Higher-order shear deformation theory (lý thuyết biến dạng cắt bậc 
cao) 
TSDT Third-order shear deformation theory (lý thuyết biến dạng cắt bậc 
ba của Reddy) 
HSDST-4 Lý thuyết vỏ bậc cao bốn ẩn chuyển vị 
GT Giải tích 
PTHH Phần tử hữu hạn 
SSSS Điều kiện biên bốn cạnh liên kết khớp 
SCSC Điều kiện biên: Khớp – Ngàm – Khớp – Ngàm 
CCCC Điều kiện biên bốn cạnh liên kết ngàm 
CFFF Điều kiện biên ngàm một cạnh (con-xon) 
CYL Vỏ trụ 
SPH Vỏ cầu 
HPR Vỏ yên ngựa 
PLA Tấm 
x 
DANH MỤC BẢNG BIỂU 
Bảng 4.1. Tần số dao động tự do không thứ nguyên thứ nhất (m=1, n=1) của mảnh 
vỏ trụ thoải composite [90/0] (a=b, Rx/a= ) ............................................................ 73 
Bảng 4.2. Tần số dao động tự do không thứ nguyên 
mn mnh G  của tấm vuông 
đẳng hướng bốn biên tựa khớp (a=b, a/h=10) .......................................................... 75 
Bảng 4.3. Tần số dao động riêng thứ nhất (m=1, n=1)  (Hz) của vỏ thoải đẳng hướng 
áp điện khi mạch kín (a/b=1, Ry/a=5) ....................................................................... 77 
Bảng 4.4. Tần số dao động riêng thứ nhất (m=1,n=1)  (Hz) của vỏ thoải đẳng hướng 
áp điện khi mạch hở (a/b=1, Ry/a=5) ................................................................. ... 
 % 2 - FG-V 
 % 3 - FG-A 
 % 4 - FG-X 
 % 5 - FG-O 
%------------------------------------------------------------------------ 
CNTP=[E11_cnt,E22_cnt,G12_cnt,nuy12_cnt,ro_cnt,V_cnts,nheta1,nheta2,nheta
3]; 
%------------------------------------------------------------------------ 
nuym=0.34; % He so Poisson-metal 
rom=1150; % kg/m3 
Em=(3.52-0.0034*300)*1e9; % Pa 
METP=[Em,nuym,rom]; 
%------------------------------------------------------------------------ 
PL23 
% Piezoelectric 
%------------------------------------------------------------------------ 
E_pie=63e9; % E-piezo 
nuy_pie=0.35; 
G12_pie=23.2e9; 
ro_pie=7750; 
e31e=-7.209; e33e=15.12; e15e=12.322; p11=1.53e-8; p33=1.5e-8; p22=p11; 
PIEP=[E_pie,nuy_pie,G12_pie,ro_pie]; 
%------------------------------------------------------------------------ 
Zpt=[2*z-hp 0 0; 
 0 2*z-hp 0; 
 0 0 1/hp]; 
Zpb=[2*z+hp 0 0; 
 0 2*z+hp 0; 
 0 0 1/hp]; 
%------------------------------------------------------------------------ 
fz=z*(-1/8+3/2*(z/(h+2*hp))^2); 
gz=1-diff(fz,z); 
Nlayer=length(Angleply); 
hz=linspace(-h/2,h/2,Nlayer+1); 
%------------------------------------------------------------------------ 
% Piezo 
Delta=(1-3*nuy_pie^2-2*nuy_pie^3)/E_pie^3; 
c11e=(1-nuy_pie^2)/(E_pie^2*Delta); 
c12e=(nuy_pie+nuy_pie^2)/(E_pie^2*Delta); 
c13e=c12e; 
c33e=(1-nuy_pie^2)/(E_pie^2*Delta); 
c55e=c33e; 
%------------------------------------------------------------------------ 
% Cac he so dan hoi 2 lop ap dien 
%------------------------------------------------------------------------ 
C11_ng=c11e-c13e^2/c33e; 
C12_ng=c12e-c13e^2/c33e; 
e31_ng=e31e-c13e/c33e*e33e; 
p33_ng=p33+e33e^2/c33e; 
mt_e=[0 0 e31_ng; 
 0 0 e31_ng; 
 0 0 0]; 
mt_es=[-e15e 0 0; 
 0 -e15e 0]; 
mt_p=[p11 0 0; 
 0 p22 0; 
 0 0 p33_ng]; 
%------------------------------------------------------------------------ 
Apt=double(int(mt_e*Zpt,z,h/2,h/2+hp)); 
Apb=double(int(mt_e*Zpb,z,-h/2-hp,-h/2)); 
Bpt=double(int(mt_e*Zpt*z,z,h/2,h/2+hp)); 
Bpb=double(int(mt_e*Zpb*z,z,-h/2-hp,-h/2)); 
Cpt=double(int(mt_e*Zpt*fz,z,h/2,h/2+hp)); 
Cpb=double(int(mt_e*Zpb*fz,z,-h/2-hp,-h/2)); 
%------------------------------------------------------------------------ 
mt_Dpb=zeros(9,6); 
mt_Dpb(1:3,1:3)=Apt; 
mt_Dpb(1:3,4:6)=Apb; 
mt_Dpb(4:6,1:3)=Bpt; 
mt_Dpb(4:6,4:6)=Bpb; 
mt_Dpb(7:9,1:3)=Cpt; 
mt_Dpb(7:9,4:6)=Cpb; 
PL24 
%------------------------------------------------------------------------ 
DDpt=double(int(mt_es*Zpt*gz,z,h/2,h/2+hp)); 
DDpb=double(int(mt_es*Zpb*gz,z,-h/2-hp,-h/2)); 
mt_Dps=zeros(2,6); 
mt_Dps(1:2,1:3)=DDpt; 
mt_Dps(1:2,4:6)=DDpb; 
%------------------------------------------------------------------------ 
ADpt=double(int(Zpt*mt_p*Zpt,z,h/2,h/2+hp)); 
ADpb=double(int(Zpb*mt_p*Zpb,z,-h/2-hp,-h/2)); 
mt_Dpp=zeros(6,6); 
mt_Dpp(1:3,1:3)=ADpt; 
mt_Dpp(4:6,4:6)=ADpb; 
%------------------------------------------------------------------------ 
NodeDof=8; % So bac tu do tren 1 nut 
nnel=4; %So nut cua phan tu 
%------------------------------------------------------------------------ 
Me=zeros(NodeDof*nnel,NodeDof*nnel); %Ma tran khoi luong phan tu 
Ke=zeros(NodeDof*nnel,NodeDof*nnel); %Ma tran do cung phan tu 1 
%------------------------------------------------------------------------ 
fxy=-(1/(2*Rxs)*(Xs-as/2)^2+1/(2*Rys)*(Ys-bs/2)^2); 
%------------------------------------------------------------------------ 
Xnodes = Xelement+1; % Number of nodes at X direction 
Ynodes = Yelement+1; % Number of nodes at Y direction 
nH = linspace(0,a,Xnodes); 
nT = linspace(0,b,Xnodes); 
[XX YY] = meshgrid(nH,nT); 
[m,n]=size(XX); 
for i=1:m 
 for j=1:n 
ZZ(i,j)=double(subs(fxy,{as,bs,Xs,Ys,Rxs,Rys},{a,b,XX(i,j),YY(i,j),Rx,Ry}
)); 
 end 
end 
%------------------------------------------------------------------------ 
[Ex,Ey,Ez,Node_number,Element_number,Edof,Dof,Coord,Tson]=Meshing4NodeShe
llFxy(XX,YY,ZZ,NodeDof); 
%------------------------------------------------------------------------ 
[AA,BB,DD,BBs,DDs,HHs,QQs,mm]=StressResultantFGCNT(CNTP,METP,TYPE,Anglepl
y,fz,hz); 
[Apie,Bpie,Dpie,Bspie,Dspie,Hspie,Qspie,mpie]=StressResultantPIEZO(PIEP,f
z,h,hp); 
AA=AA+Apie; BB=BB+Bpie; DD=DD+Dpie; BBs=BBs+Bspie; 
DDs=DDs+Dspie; HHs=HHs+Hspie; QQs=QQs+Qspie; 
Db=[AA BB BBs; 
 BB DD DDs; 
 BBs DDs HHs]; 
Ds=[QQs]; 
mm=mm+mpie; 
%------------------------------------------------------------------------ 
Kg=zeros(Node_number*NodeDof,Node_number*NodeDof); % MT do cung tong the 
Kpg=zeros(Node_number*NodeDof,Node_number*NodeDof); % MT do cung tong the 
Mg=zeros(Node_number*NodeDof,Node_number*NodeDof); % MT KL tong the 
%------------------------------------------------------------------------ 
% DIEU KIEN BIEN 
%------------------------------------------------------------------------
% Bc=SSSS(a,b,Coord); % Dieu kien bien 
Bc=CCCC(a,b,Coord); % Dieu kien bien 
PL25 
Bcc=Bc(:,1); 
%------------------------------------------------------------------------ 
ir=2; 
irs=2; 
ExLocal=zeros(Element_number,4); 
EyLocal=zeros(Element_number,4); 
for i=1:Element_number 
 dai=sqrt((Ex(i,2)-Ex(i,1))^2+(Ey(i,2)-Ey(i,1))^2+(Ez(i,2)-
Ez(i,1))^2); 
 ExLocal(i,:)=[0 dai dai 0]; 
 rong=sqrt((Ex(i,4)-Ex(i,1))^2+(Ey(i,4)-Ey(i,1))^2+(Ez(i,4)-
Ez(i,1))^2); 
 EyLocal(i,:)=[0 0 rong rong]; 
end 
%------------------------------------------------------------------------ 
for iel=1:Element_number 
[Ke,Me]=Stiffness_Mass4(ExLocal(iel,:),EyLocal(iel,:),Db,Ds,ir,irs,mm); 
 Ke=Tson(:,:,iel)'*Ke*Tson(:,:,iel); 
 Me=Tson(:,:,iel)'*Me*Tson(:,:,iel); 
[Ke_uphi,Ke_phiphi]=StiffnessPie4(ExLocal(iel,:),EyLocal(iel,:),mt_Dpb,mt
_Dps,mt_Dpp,ir,irs); 
 Kepie=Ke_uphi*inv(Ke_phiphi)*Ke_uphi'; 
 Kepie=Tson(:,:,iel)'*Kepie*Tson(:,:,iel); 
 [Kg]=assem(Edof(iel,:),Kg,Ke); 
 [Kpg]=assem(Edof(iel,:),Kpg,Kepie); 
 [Mg]=assem(Edof(iel,:),Mg,Me); 
end 
%------------------------------------------------------------------------ 
Kg=Kg+Kpg; 
Kg=sparse(Kg); 
Mg=sparse(Mg); 
%------------------------------------------------------------------------ 
Sotanso=1; 
[La,Egv]=eigen(Kg,Mg,Bcc,Sotanso); 
%------------------------------------------------------------------------ 
% Ket qua tan so goc dao dong rieng 
%------------------------------------------------------------------------ 
Omega =sqrt(La)/2/pi 
%------------------------------------------------------------------------ 
Phụ lục B5: Chương trình tính để giải bài toán đáp ứng chuyển vị theo thời gian 
bằng phương pháp PTHH 
VibControl_FE_DCurvedShell_PFG_CNTRC 
% Problem: Dieu khien dao dong vo thoai hai do cong composite PFG-CNTRC 
% PTHH 
% Clear memory 
close all; clear all; format long; clc; 
%------------------------------------------------------------------------ 
syms z Xs Ys Zs fxy Rxs Rys as bs 
%------------------------------------------------------------------------ 
tt=1e-3; 
hc=3*tt; %Chieu day tong cac lop CNT 
PL26 
hp=250e-6; %Chieu day lop Peizo 
ha=20; % Ty so a/h 
ba=1; 
h_tam=hc+2*hp; 
a=hc*ha; % Canh a 
b=a*ba; % Canh b 
Rxa=10; 
Ryb=10; 
Angleply=[-45 45 -45]*pi/180; %Goc soi 
fz=z*(-1/8+3/2*(z/h_tam)^2); 
%------------------------------------------------------------------------ 
Xelement=20; 
Yelement=20; 
%------------------------------------------------------------------------ 
Rx=Rxa*a; 
Ry=Ryb*b; 
%------------------------------------------------------------------------ 
% Co tinh CNT 
%------------------------------------------------------------------------ 
TYPE=5; % Kieu phan bo CNT 
 % 1 - UD 
 % 2 - FG-V 
 % 3 - FG-A 
 % 4 - FG-X 
 % 5 - FG-O 
%------------------------------------------------------------------------ 
nuy12_cnt=0.175; % Table 1 
ro_cnt=1400; % BB Piezo CNT Reddy 
E11_cnt=5.64e12; E22_cnt=7.0800e12; G12_cnt=1.9455e12; 
%------------------------------------------------------------------------ 
V_cnts=0.12; nheta1=0.137; nheta2=1.022; nheta3=0.7*nheta2; 
%V_cnts=0.17; nheta1=0.142; nheta2=1.626; nheta3=0.7*nheta2; 
%V_cnts=0.28; nheta1=0.141; nheta2=1.585; nheta3=0.7*nheta2; 
%------------------------------------------------------------------------ 
CNTP=[E11_cnt,E22_cnt,G12_cnt,nuy12_cnt,ro_cnt,V_cnts,nheta1,nheta2,nheta
3]; 
%------------------------------------------------------------------------ 
nuym=0.34; % He so Poisson-metal 
rom=1150; % kg/m3 
Em=(3.52-0.0034*300)*1e9; % Pa 
METP=[Em,nuym,rom]; 
%------------------------------------------------------------------------ 
% Vat lieu Piezoelectric 
E_pie=63e9; % E-piezo 
nuy_pie=0.35; 
G12_pie=23.2e9; 
ro_pie=7750; 
e31e=-7.209; e33e=15.12; e15e=12.322; p11=1.53e-8; p33=1.5e-8; p22=p11; 
PIEP=[E_pie,nuy_pie,G12_pie,ro_pie]; 
%------------------------------------------------------------------------ 
Zpt=[2*z-hp 0 0; 
 0 2*z-hp 0; 
 0 0 1/hp]; 
Zpb=[2*z+hp 0 0; 
 0 2*z+hp 0; 
 0 0 1/hp]; 
%------------------------------------------------------------------------ 
fz=z*(-1/8+3/2*(z/(h+2*hp))^2); 
PL27 
gz=1-diff(fz,z); 
Nlayer=length(Angleply); 
hz=linspace(-h/2,h/2,Nlayer+1); 
%------------------------------------------------------------------------ 
% Piezo 
Delta=(1-3*nuy_pie^2-2*nuy_pie^3)/E_pie^3; 
c11e=(1-nuy_pie^2)/(E_pie^2*Delta); 
c12e=(nuy_pie+nuy_pie^2)/(E_pie^2*Delta); 
c13e=c12e; 
c33e=(1-nuy_pie^2)/(E_pie^2*Delta); 
c55e=c33e; 
%------------------------------------------------------------------------ 
% Cac he so dan hoi 2 lop ap dien 
%------------------------------------------------------------------------ 
C11_ng=c11e-c13e^2/c33e; 
C12_ng=c12e-c13e^2/c33e; 
e31_ng=e31e-c13e/c33e*e33e; 
p33_ng=p33+e33e^2/c33e; 
mt_e=[0 0 e31_ng; 
 0 0 e31_ng; 
 0 0 0]; 
mt_es=[-e15e 0 0; 
 0 -e15e 0]; 
mt_p=[p11 0 0; 
 0 p22 0; 
 0 0 p33_ng]; 
%------------------------------------------------------------------------ 
Apt=double(int(mt_e*Zpt,z,h/2,h/2+hp)); 
Apb=double(int(mt_e*Zpb,z,-h/2-hp,-h/2)); 
Bpt=double(int(mt_e*Zpt*z,z,h/2,h/2+hp)); 
Bpb=double(int(mt_e*Zpb*z,z,-h/2-hp,-h/2)); 
Cpt=double(int(mt_e*Zpt*fz,z,h/2,h/2+hp)); 
Cpb=double(int(mt_e*Zpb*fz,z,-h/2-hp,-h/2)); 
%------------------------------------------------------------------------ 
mt_Dpb=zeros(9,6); 
mt_Dpb(1:3,1:3)=Apt; 
mt_Dpb(1:3,4:6)=Apb; 
mt_Dpb(4:6,1:3)=Bpt; 
mt_Dpb(4:6,4:6)=Bpb; 
mt_Dpb(7:9,1:3)=Cpt; 
mt_Dpb(7:9,4:6)=Cpb; 
%------------------------------------------------------------------------ 
DDpt=double(int(mt_es*Zpt*gz,z,h/2,h/2+hp)); 
DDpb=double(int(mt_es*Zpb*gz,z,-h/2-hp,-h/2)); 
mt_Dps=zeros(2,6); 
mt_Dps(1:2,1:3)=DDpt; 
mt_Dps(1:2,4:6)=DDpb; 
%------------------------------------------------------------------------ 
ADpt=double(int(Zpt*mt_p*Zpt,z,h/2,h/2+hp)); 
ADpb=double(int(Zpb*mt_p*Zpb,z,-h/2-hp,-h/2)); 
mt_Dpp=zeros(6,6); 
mt_Dpp(1:3,1:3)=ADpt; 
mt_Dpp(4:6,4:6)=ADpb; 
%------------------------------------------------------------------------ 
NodeDof=8; % So bac tu do tren 1 nut 
nnel=4; % So nut cua phan tu 
%------------------------------------------------------------------------ 
Me=zeros(NodeDof*nnel,NodeDof*nnel); %Ma tran khoi luong phan tu 
PL28 
Ke_uu=zeros(NodeDof*nnel,NodeDof*nnel); %Ma tran do cung phan tu 1 
%------------------------------------------------------------------------ 
% Bat dau chia luoi tu dong cho loai phan tu 4 nut 
%------------------------------------------------------------------------ 
fxy=-(1/(2*Rxs)*(Xs-as/2)^2+1/(2*Rys)*(Ys-bs/2)^2); 
%------------------------------------------------------------------------ 
Xnodes = Xelement+1; % Number of nodes at X direction 
Ynodes = Yelement+1; % Number of nodes at Y direction 
nH = linspace(0,a,Xnodes); 
nT = linspace(0,b,Xnodes); 
[XX YY] = meshgrid(nH,nT); 
[m,n]=size(XX); 
for i=1:m 
 for j=1:n 
ZZ(i,j)=double(subs(fxy,{as,bs,Xs,Ys,Rxs,Rys},{a,b,XX(i,j),YY(i,j),Rx,Ry}
)); 
 end 
end 
[Ex,Ey,Ez,Node_number,Element_number,Edof,~,Coord,Tson]=Meshing4NodeShell
Fxy(XX,YY,ZZ,NodeDof); 
%------------------------------------------------------------------------ 
[AA,BB,DD,BBs,DDs,HHs,QQs,mm]=StressResultantFGCNT(CNTP,METP,TYPE,Anglepl
y,fz,hz); 
[Apie,Bpie,Dpie,Bspie,Dspie,Hspie,Qspie,mpie]=StressResultantPIEZO(PIEP,f
z,h,hp); 
AA=AA+Apie; BB=BB+Bpie; DD=DD+Dpie; BBs=BBs+Bspie; 
DDs=DDs+Dspie; HHs=HHs+Hspie; QQs=QQs+Qspie; 
Db=[AA BB BBs; 
 BB DD DDs; 
 BBs DDs HHs]; 
Ds=[QQs]; 
mm=mm+mpie; 
%------------------------------------------------------------------------ 
Kg=zeros(Node_number*NodeDof,Node_number*NodeDof); % Ma tran do cung 
tong the 
Kpg=zeros(Node_number*NodeDof,Node_number*NodeDof); % Ma tran do cung 
tong the 
Mg=zeros(Node_number*NodeDof,Node_number*NodeDof); % Ma tran khoi luong 
tong the 
Fg=zeros(Node_number*NodeDof,1); % Ma tran khoi luong tong the 
%------------------------------------------------------------------------ 
% DIEU KIEN BIEN 
%------------------------------------------------------------------------ 
Bc=CCCC(a,b,Coord); % Dieu kien bien 
Bcc=Bc(:,1); 
%------------------------------------------------------------------------ 
q0=1e4; 
LoadType=1; % Dang phan bo 
 % 1 - phan bo hinh sin 
 % 2 - phan bo deu 
Taitrongcohoc=[q0 LoadType a b]; 
ir=2; 
irs=2; 
ExLocal=zeros(Element_number,4); 
EyLocal=zeros(Element_number,4); 
for i=1:Element_number 
PL29 
 dai=sqrt((Ex(i,2)-Ex(i,1))^2+(Ey(i,2)-Ey(i,1))^2+(Ez(i,2)-
Ez(i,1))^2); 
 ExLocal(i,:)=[0 dai dai 0]; 
 rong=sqrt((Ex(i,4)-Ex(i,1))^2+(Ey(i,4)-Ey(i,1))^2+(Ez(i,4)-
Ez(i,1))^2); 
 EyLocal(i,:)=[0 0 rong rong]; 
end 
%------------------------------------------------------------------------ 
for iel=1:Element_number 
[Ke_uu,Me,Fe]=Stiffness_Mass4(ExLocal(iel,:),EyLocal(iel,:),Db,Ds,ir,irs,
mm,Taitrongcohoc); 
 Ke=Tson(:,:,iel)'*Ke_uu*Tson(:,:,iel); 
 Me=Tson(:,:,iel)'*Me*Tson(:,:,iel); 
C1=StiffnessPie4as(ExLocal(iel,:),EyLocal(iel,:),mt_Dpb,mt_Dps,mt_Dpp,ir,
irs); 
 C1=Tson(:,:,iel)'*C1*Tson(:,:,iel); 
 [Kg,Fg]=assem(Edof(iel,:),Kg,Ke,Fg,Fe); 
 [Kpg]=assem(Edof(iel,:),Kpg,C1); 
 [Mg]=assem(Edof(iel,:),Mg,Me); 
end 
Ksao=Kg+Gd*Kpg; 
CC=-Gv*Kpg; 
Kg=sparse(Kg); 
Mg=sparse(Mg); 
Ksao=sparse(Ksao); 
%---------------------------------------------------------------------- 
Sotanso=2; 
[La,~]=eigen(Kg,Mg,Bcc,Sotanso); 
w1=sqrt(La(1))/2/pi; 
w2=sqrt(La(2))/2/pi; 
syms a1 a2 
Eq1=a1+w1^2*a2-2*Tilecan*w1; 
Eq2=a1+w2^2*a2-2*Tilecan*w2; 
S11=solve(Eq1,Eq2,a1,a2); 
a11=double(S11.a1); 
a22=double(S11.a2); 
Can=[a11 a22]; 
DetalT=0.00001; 
T1=0.2e-2; 
tf=1.4e-2; 
nt1=round(T1/DetalT); 
nt=round(tf/DetalT); 
%---------------------------------------------------------------------- 
nutgiua=(Node_number+1)/2; 
BTD=(nutgiua-1)*8+3; 
%---------------------------------------------------------------------- 
[u,~,~,t] = Newmark4(Mg,Ksao,CC,Can,Fg,tf,nt,Bc,nt1); 
TT=t; 
UU=u(BTD,:); 
%----------------------------------------------------------------------