The GMC micromechanical model was employed successfully to calculate the thermal residual stress of the fiber composites with different fiber arrays, i.e., square edge packing, square diagonal packing, and hexagonal packing, during the cooling process. Based on the micromechanical analysis, the nonlinear mechanical behaviors of the fiber composites in the presence of the thermal residual stress effect were determined. Results indicated that for the composites with square edge packing, the constitutive behaviors are affected appreciably by the thermal residual stress.
However, for the composites with hexagonal packing and square diagonal packing, the thermal residual stress exhibits little effects on their properties.
With regard to the damping properties, it was observed that the composite structures constructed based on square diagonal packing fiber array demonstrate superior vibration damping properties than the other two cases. This phenomenon could be due to the fact that in the composites with square diagonal packing, the matrix phase with higher damping capacity dissipates more strain energy under one cycle of vibration.
In addition, from the comparison of the GMC and SCMC micromechanical models, it was suggested that although more constrains were implemented in the GMC model inducing the shear coupling effect, GMC model still exhibit the advantage of converging efficiency during the micromechanical analysis. As a result, the GMC micromechanical model still can be employed efficiently for characterizing the mechanical behaviors of fiber composites in terms of the geometry of microstructure.
Reference
[1] Paley, M. and Aboudi, J. 1992 “Micromechanical Analysis of Composites by the Generalized Cells Model,” Mechanics of Materials, Vol. 14, No. 2, pp.
127-139.
[2] Sun, C. T. and Vaidya, R. S. 1996 “Prediction of Composite Properties from A Representative Volume element,” Composites Science and Technology, Vol. 56, No. 2, pp. 171-179.
[3] Zhu, C. and Sun, C. T. 2003 “Micromechanical Modeling of Fiber Composites under Off-Axis Loading,” Journal of Thermoplastic Composite Materials, Vol.
16, No. 4, pp. 333-344.
[4] Hsu, S.Y., Vogler, T. J. and Kyriakides, S. 1999 “Inelastic Behavior of An AS4/PEEK Composite under Combined Transverse Compression and Shear Park II: Modeling,” International Journal of Plasticity, Vol. 15, No. 8, pp.
807-836.
[5] Orozco, C. E. and Pindera, M. J. 1999 “Plastic Analysis of Complex Microstructure Composites Using the Generalized Method of Cells,” AIAA Journal, Vol. 37, No. 4, pp. 482-488.
[6] Pindera, M. J. and Bednarcyk, B. A. 1999 “An Efficient Implementation of the Generalized Method of Cells for Unidirectional, Multi-phased Composites with Complex Microstructures,” Composites Part B, Vol. 30, No. 1, pp. 87-105.
[7] Pindera, M. J., Herakovich, C. T., Becker, W. and Aboudi, J. 1990 “Nonlinear Response of Unidirectional Boron/Aluminum,” Journal of Composite Materials, Vol. 24, No. 1, pp. 2-21.
[8] Zhou, F. H., Hashimoto, R., Ogawa, A. and Sofue, Y. 2006 “Residual Stress and Its Effect on Yielding in SiC/Ti Plate,” JSME International Journal Series A - Solid Mechanics & Material Engineering, Vol. 49, No. 1, pp. 25-31.
[9] Aghdam, M. M., Pavier, M. J. and Smith, D. J. 2001 “Micromechanics of Off-Axis Loading of Metal Matrix Composites Using Finite Element Analysis,”
International Journal of Solids and Structures, Vol. 38, No. 22-23, pp.
3905-3925.
[10] Arnold, S. M., Pindera, M. J. and Wilt, T. E. 1996 “Influence Fiber Architecture on the Inelastic Response of Metal Matrix Composites,” International Journal of Plasticity, Vol. 12, No. 4, pp. 507-545.
[11] Saravanos, D. A. and Chamis, C. C. 1990 “Unified Micromechanics of Damping for Unidirectional and Off-Axis Fiber Composites,” Journal of Composites Technology and Research, Vol. 12, No. 1, pp. 31-40.
[12] Hwang, S. J. and Gibson, R. F. 1993 “Prediction of Fiber-Matrix Interphase Effects on Damping of Composites Using A Micromechanical Strain Energy/Finite Element Approach,” Composites Engineering, Vol. 3, No. 10, pp.
975-984.
[13] Kaliske, M. and Rothert, H. 1995 “Damping Characterization of Unidirectional Fiber Reinforced Polymer Composites,” Composites Engineering, Vol. 5, No.
5, pp. 551-567.
[14] ANSYS Element Reference. 000855. English Edition. SAS IP, Inc. 1997.
[15] Gan, H., Orozco, C. E. and Herakovich, C. T. 2000 “A Strain-Ccompatible Method for Micromechanical Analysis of Multi-Phase Composites,”
International Journal of Solids and Structure, Vol. 37, No. 37, pp. 5097-5122.
[16] Orozco, C. E. and Gan, H. 2002 “Viscoplastic Response of Multiphase Composites Using A Strain-Ccompatible Volume-Aaverage Method,”
Composites: Part B, Vol. 33, No. 4, pp. 301-313.
[17] Hanselka, H. 1993 Ein Beitrag Zur Charakterisierung Des Dämpfungsverhaltens Polymer Faserverbound-Werksofft, Dissertation, Technische Universität Clausthal, Germany.
Appendix A MATLAB code for calculating the thermal residual stresses using GMC analysis
clear all
% =============================== something to input
==============================
% material properties (transversely isotropic fiber + isotropic matrix) Ef1=234000; % 11 Young's modulus of fiber
Ef2=14000; % 22 Young's modulus of fiber Ef3=Ef2;
nuf12=0.2; % 12 Poisson's ratio of fiber nuf13=0.2;
nuf23=0.25;
Gf12=27600;%Ef1/(2*(1+nuf12)); % 12 shear modulus of fiberGf13=Gf12;
Gf23=5500;
Em=3400; % Young's Modulus of matrix num=0.37; % Poisson's ratio of matrix Gm=Em/(2*(1+num)); % Shear modulus of matrix
T=-200; % temperature change af=0; % thermal coefficient of fiber am=1.18e-4; % thermal coefficient of matrix
% Geometry information
n=3 ; %================= (1)2*2 (2)square random (3)hexagonal======================
if n==1 h=1;
l=1;
volf=0.8;
Nb=2; % divided portions in the beta direction Nr=2; % divided portions in the gama direction bb=sqrt(volf);
hh=1-bb;
hb=[bb hh]; % h_beta hr=hb; % h_gama Nfiber=1;
regionf=[1]; % denote what cells are fiber
%end of code elseif n==2 h=1;
l=1;
load hb.txt;
hb=hb';
hr=hb;
Nb=size(hb);
Nb=Nb(1,2);
Nr=Nb;
load regionf.txt;
regionf=regionf';
Nfiber=size(regionf);
Nfiber=Nfiber(1,1);
elseif n==3 h=sqrt(3);
l=1;
load hb.txt;
load hr.txt;
hb=hb';
hr=hr';
Nb=size(hb);
Nr=size(hr);
Nb=Nb(1,2);
Nr=Nr(1,2);
load regionf.txt;
regionf=regionf';
Nfiber=size(regionf);
Nfiber=Nfiber(1,1);
end
% ================================== matrix deploy
==================================
Sm=zeros(6,6);
AG=zeros(2*(Nb+Nr)+Nb*Nr+1,6*Nb*Nr);
J=zeros(2*(Nb+Nr)+Nb*Nr+1,6);
AVPM=zeros(5*Nb*Nr-2*(Nb+Nr)-1,6*Nb*Nr);
Cs=zeros(6*Nb*Nr,6); % sub-cell stiffness matrix AVP=zeros(6*Nb*Nr,6);
cs=zeros(6*Nb*Nr,6);
aT=zeros(6*Nb*Nr,1);
displ=0;
% =================== current stress & strain independent matrix ====================
% fiber stiffness matrix nuf21=Ef2*nuf12/Ef1;
nuf31=Ef3*nuf13/Ef1;
nuf32=Ef3*nuf23/Ef2;
Sf=[1/Ef1 -nuf21/Ef2 -nuf31/Ef3 0 0 0; -nuf12/Ef1 1/Ef2 -nuf32/Ef3 0 0 0;...
-nuf13/Ef1 -nuf23/Ef2 1/Ef3 0 0 0; 0 0 0 1/Gf23 0 0; 0 0 0 0 1/Gf13 0; 0 0 0 0 0 1/Gf12];
Cf=inv(Sf);
Sm=[1/Em -num/Em -num/Em 0 0 0; -num/Em 1/Em -num/Em 0 0 0;...
-num/Em -num/Em 1/Em 0 0 0; 0 0 0 1/Gm 0 0; 0 0 0 0 1/Gm 0; 0 0 0 0 0 1/Gm];
Cm=inv(Sm);
% ============================= main loop and code
===================================
% AG matrix index=-5;
for i=1:Nb*Nr % 11 displ. continuity
index=index+6;
AG(i,index)=1;
J(i,1)=1;
end
for gama=1:Nr % loop to apply displ. continuity for beta=1:Nb
index=Nr*beta-Nr+gama;
AG(Nb*Nr+gama,(index-1)*6+2)=hb(1,beta); % 22 J(Nb*Nr+gama,2)=h;
AG(Nb*Nr+Nr+beta,(index-1)*6+3)=hr(1,gama); %33 J(Nb*Nr+Nr+beta,3)=l;
AG(Nb*Nr+Nr+Nb+1,(index-1)*6+4)=hb(1,beta)*hr(1,gama)/2; %23
J(Nb*Nr+Nr+Nb+1,4)=h*l/2;
AG(Nb*Nr+Nr+Nb+1+beta,(index-1)*6+5)=hr(1,gama)/2; %13 J(Nb*Nr+Nr+Nb+1+beta,5)=l/2;
AG(Nb*Nr+Nr+2*Nb+1+gama,(index-1)*6+6)=hb(1,beta)/2; %12 J(Nb*Nr+Nr+2*Nb+1+gama,6)=h/2;
end end
% K matrix
K=[zeros(5*Nb*Nr-2*(Nb+Nr)-1,6); J];
for i=1:Nfiber;
cs(regionf(i,1)*6-5,1)=1;
end
for beta=1:Nb;
for gama=1:Nr;
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
a1=6*(index-1)+1;
a2=6*(index-1)+3;
if cs(index*6-5,1)==1;
Cs(p1:p2,:)=Cf;
aT(a1:a2,:)=af*T;
else
Cs(p1:p2,:)=Cm;
aT(a1:a2,:)=am*T;
end
end end count=0;
for beta=1:Nb-1 % construct AVPM matrix for gama=1:Nr
betahat=beta+1;
index=(beta-1)*Nr+gama-1;
index2=(betahat-1)*Nr+gama-1;
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+2,j); % 22 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+2,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+6,j); % 12 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+6,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+4,j); % 23 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end end
for gama=1:Nr-1 % construct AVPM matrix for beta=1:Nb
gamahat=gama+1;
index=(beta-1)*Nr+gama-1;
index2=(beta-1)*Nr+gamahat-1;
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+3,j); % 33 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+3,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+5,j); % 13 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+5,j);
end end end
for gama=1:Nr-1 count=count+1;
beta=Nb;
gamahat=gama+1;
index=(beta-1)*Nr+gama-1;
index2=(beta-1)*Nr+gamahat-1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+4,j); % 23 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end
AtVP=[AVPM; AG];
AVP=AtVP\K;
Bvp=zeros(6,6);
APh=[AVPM;zeros(2*(Nb+Nr)+Nb*Nr+1,6*Nb*Nr)];
AP=AtVP\APh;
APaT=AP*aT;
for beta=1:Nb for gama=1:Nr
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
Bvp=Bvp+hb(1,beta)*hr(1,gama)*Cs(p1:p2,:)*AVP(p1:p2,:);
end end
Bvp=Bvp/h/l;
S=inv(Bvp);
E1=1/S(1,1);
E2=1/S(2,2);
G12=1/S(6,6);
G23=1/S(4,4);
nu12=S(2,1)*(-E1);
%end
%=========================== thermal strain =======================
afT=zeros(6,1);
for beta=1:Nb;
for gama=1:Nr;
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
afT=-S*hb(1,beta)*hr(1,gama)*Cs(p1:p2,:)*(APaT(p1:p2,:)-aT(p1:p2,:))+afT;
end end
avT=-S*afT/(h*l);
deps=avT;
%========================== thermal residual stresses=========================
AtVP=[AVPM; AG];
APh=[AVPM;zeros(2*(Nb+Nr)+Nb*Nr+1,6*Nb*Nr)];
AVP=AtVP\K;
AP=AtVP\APh;
APaT=AP*aT;
for beta=1:Nb;
for gama=1:Nr;
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
sigma(p1:p2,:)=Cs(p1:p2,:)*(AVP(p1:p2,:)*deps+APaT(p1:p2,:)-aT(p1:p2,:));
end end
Appendix B MATLAB code for GMC analysis
clear all tic
% =============================== something to input
==============================
% material properties (transversely isotropic fiber + isotropic matrix) Ef1=234000; % 11 Young's modulus of fiber
Ef2=14000; % 22 Young's modulus of fiber Ef3=Ef2;
nuf12=0.2; % 12 Poisson's ratio of fiber nuf13=0.2;
nuf23=0.25;
Gf12=27600;%Ef1/(2*(1+nuf12)); % 12 shear modulus of fiber, Ef/(2*(1+nuf)) - if fibers are isotropic
Gf13=Gf12;
Gf23=5500;
Em=3400; % Young's Modulus of matrix num=0.37; % Poisson's ratio of matrix Gm=Em/(2*(1+num)); % Shear modulus of matrix
% Power law coefficient
b=6.42e-11; % power law coefficient beta pn=4.11; % power law coefficient n
% Geometry information
n=1 ; %================= (1)2*2 (2)square random (3)hexagonal======================
if n==1 h=1;
l=1;
volf=0.6;
Nb=2; % divided portions in the beta direction Nr=2; % divided portions in the gama direction bb=sqrt(volf);
hh=1-bb;
hb=[bb hh]; % h_beta hr=hb; % h_gama
Nfiber=1;
regionf=[1]; % denote what cells are fiber
%end of code elseif n==2 h=1;
l=1;
load hb.txt;
hb=hb';
hr=hb;
Nb=size(hb);
Nb=Nb(1,2);
Nr=Nb;
load regionf.txt;
regionf=regionf';
Nfiber=size(regionf);
Nfiber=Nfiber(1,1);
elseif n==3 h=1;
l=sqrt(3);
load hb.txt;
load hr.txt;
hb=hb';
hr=hr';
Nb=size(hb);
Nr=size(hr);
Nb=Nb(1,2);
Nr=Nr(1,2);
load regionf.txt;
regionf=regionf';
Nfiber=size(regionf);
Nfiber=Nfiber(1,1);
end
%=================================Step====================================
= step=80;
Data=zeros(step+1,12);
bb=0;
for aaa=7;
if aaa==1 | 2 | 3 | 4 | 5 | 7 bb=bb+2;
bbb=bb-1;
%====================================Angle===============================
==========
% Loop & angle Nstep=step;
delta=2; % stress increment angle=(aaa-1)*15;
rad=angle/180*pi;
% ================================== matrix deploy
==================================
Sm=zeros(6,6);
AG=zeros(2*(Nb+Nr)+Nb*Nr+1,6*Nb*Nr);
J=zeros(2*(Nb+Nr)+Nb*Nr+1,6);
subStrain=zeros(6*Nb*Nr,1);
subStress=zeros(6*Nb*Nr,1);
deps=zeros(6,1);
AVPM=zeros(5*Nb*Nr-2*(Nb+Nr)-1,6*Nb*Nr);
sigam=zeros(6,1); % current stress
Cs=zeros(6*Nb*Nr,6); % sub-cell stiffness matrix cs=zeros(6*Nb*Nr,6);
Tstress=zeros(Nstep+1,1);
Tstrain=zeros(Nstep+1,1);
AVP=zeros(6*Nb*Nr,6);
% =================== current stress & strain independent matrix ====================
% fiber stiffness matrix nuf21=Ef2*nuf12/Ef1;
nuf31=Ef3*nuf13/Ef1;
nuf32=Ef3*nuf23/Ef2;
Sf=[1/Ef1 -nuf21/Ef2 -nuf31/Ef3 0 0 0;...
-nuf12/Ef1 1/Ef2 -nuf32/Ef3 0 0 0;...
-nuf13/Ef1 -nuf23/Ef2 1/Ef3 0 0 0;...
0 0 0 1/Gf23 0 0;...
0 0 0 0 1/Gf13 0;...
0 0 0 0 0 1/Gf12];
Cf=inv(Sf);
Sme=[1/Em -num/Em -num/Em 0 0 0;...
-num/Em 1/Em -num/Em 0 0 0;...
-num/Em -num/Em 1/Em 0 0 0;...
0 0 0 1/Gm 0 0;...
0 0 0 0 1/Gm 0;...
0 0 0 0 0 1/Gm];
% AG matrix index=-5;
for i=1:Nb*Nr % 11 displ. continuity index=index+6;
AG(i,index)=1;
J(i,1)=1;
end
for gama=1:Nr % loop to apply displ. continuity for beta=1:Nb
index=Nr*beta-Nr+gama;
AG(Nb*Nr+gama,(index-1)*6+2)=hb(1,beta); % 22 J(Nb*Nr+gama,2)=h;
AG(Nb*Nr+Nr+beta,(index-1)*6+3)=hr(1,gama); %33 J(Nb*Nr+Nr+beta,3)=l;
AG(Nb*Nr+Nr+Nb+1,(index-1)*6+4)=hb(1,beta)*hr(1,gama)/2; %23 J(Nb*Nr+Nr+Nb+1,4)=h*l/2;
AG(Nb*Nr+Nr+Nb+1+beta,(index-1)*6+5)=hr(1,gama)/2; %13 J(Nb*Nr+Nr+Nb+1+beta,5)=l/2;
AG(Nb*Nr+Nr+2*Nb+1+gama,(index-1)*6+6)=hb(1,beta)/2; %12 J(Nb*Nr+Nr+2*Nb+1+gama,6)=h/2;
end end
% Tsig & Teps
Tsig=[cos(rad)^2 sin(rad)^2 0 0 0 2*cos(rad)*sin(rad);...
sin(rad)^2 cos(rad)^2 0 0 0 -2*cos(rad)*sin(rad);...
0 0 1 0 0 0;...
0 0 0 cos(rad) -sin(rad) 0;...
0 0 0 sin(rad) cos(rad) 0;...
-cos(rad)*sin(rad) cos(rad)*sin(rad) 0 0 0 cos(rad)^2-sin(rad)^2];
Teps=[cos(rad)^2 sin(rad)^2 0 0 0 -cos(rad)*sin(rad);...
sin(rad)^2 cos(rad)^2 0 0 0 cos(rad)*sin(rad);...
0 0 1 0 0 0;...
0 0 0 cos(rad) sin(rad) 0;...
0 0 0 -sin(rad) cos(rad) 0;...
2*cos(rad)*sin(rad) -2*cos(rad)*sin(rad) 0 0 0 cos(rad)^2-sin(rad)^2];
% K matrix
K=[zeros(5*Nb*Nr-2*(Nb+Nr)-1,6); J];
for i=1:Nfiber;
cs(regionf(i,1)*6-5,1)=1;
end
for beta=1:Nb;
for gama=1:Nr;
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
if cs(index*6-5,1)==1;
Cs(p1:p2,:)=Cf;
else
Cs(p1:p2,:)=inv(Sme);
end
end end
%==========================thermal stress in GMC============================
load residualStress.txt;
subStress=residualStress;
clear residualStress;
% ============================= main loop and code
===================================
dsig1=[ delta; 0; 0; 0; 0; 0]; % !!!!!!!
dsig=Tsig*dsig1;
%dsig=dsig1;
for i=1:Nstep now=i
Tstress(i+1,1)=Tstress(i,1)+dsig1(1,1); %!!!!!!
count=0;
for beta=1:Nb-1 % construct AVPM matrix for gama=1:Nr
betahat=beta+1;
index=(beta-1)*Nr+gama-1;
index2=(betahat-1)*Nr+gama-1;
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+2,j); % 22 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+2,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+6,j); % 12 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+6,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+4,j); % 23 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end end
for gama=1:Nr-1 % construct AVPM matrix for beta=1:Nb
gamahat=gama+1;
index=(beta-1)*Nr+gama-1;
index2=(beta-1)*Nr+gamahat-1;
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+3,j); % 33 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+3,j);
end
count=count+1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+5,j); % 13 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+5,j);
end end end
for gama=1:Nr-1 count=count+1;
beta=Nb;
gamahat=gama+1;
index=(beta-1)*Nr+gama-1;
index2=(beta-1)*Nr+gamahat-1;
for j=1:6
AVPM(count,6*index+j)=Cs(6*index+4,j); % 23 traction continuity AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end
AtVP=[AVPM; AG];
AVP=AtVP\K;
Bvp=zeros(6,6);
for beta=1:Nb for gama=1:Nr
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
Bvp=Bvp+hb(1,beta)*hr(1,gama)*Cs(p1:p2,:)*AVP(p1:p2,:);
end end
Bvp=Bvp/h/l;
deps=Bvp\dsig; % engineering strain increment of composite in the material axis deps1=Teps*deps; % strain increment in the loading axis
Tstrain(i+1,1)=Tstrain(i,1)+deps1(1,1); %!!!!!!!!
%epss(i,:)=deps(:,1)';
%if Tstrain(i+1,1) >= 0.04 % break
%else %end
for i=1:Nfiber;
cs(regionf(i,1)*6-5,1)=1;
end
for beta=1:Nb;
for gama=1:Nr;
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
if cs(index*6-5,1)==1;
Cs(p1:p2,:)=Cf;
else
%leps(p1:p2,1)=AVP(p1:p2,:)*deps;
sigam=Cs(p1:p2,:)*AVP(p1:p2,:)*deps; % evaluate local stress increment %check(p1:p2,1)=sigam;
subStress(p1:p2,1)=subStress(p1:p2,1)+sigam; % current stress state sigam=subStress(p1:p2,1);
s1=1/3*(2*sigam(1,1)-sigam(2,1)-sigam(3,1));
s2=1/3*(-sigam(1,1)+2*sigam(2,1)-sigam(3,1));
s3=1/3*(-sigam(1,1)-sigam(2,1)+2*sigam(3,1));
s4=2*sigam(4,1);
s5=2*sigam(5,1);
s6=2*sigam(6,1);
efsig=sqrt((sigam(1,1)+sigam(2,1)+sigam(3,1))^2-3*(sigam(2,1)*sigam(3,1)-sigam(4,1)^2+sigam(1,1)
*sigam(3,1)-...
sigam(5,1)^2+sigam(1,1)*sigam(2,1)-sigam(6,1)^2));
%if index == 2 % test(i,1)=efsig;
%else %end
Smp=3/4*b*pn*efsig^(pn-3)*[s1; s2; s3; s4; s5; s6]*3*[s1; s2; s3; s4; s5; s6]';
%Smp(4:6,:)=Smp(4:6,:)*2;
Sm=Sme+Smp;
Cs(p1:p2,:)=inv(Sm);
end end
end end
sho=[Tstrain, Tstress];
Data(:,bbb:bb)=sho;
end end
%sho=[Tstrain, Tstress];
plot(Tstrain,Tstress,'b');
toc
Appendix C MATLAB code for calculating the modal shapes and damping capacity of composite structures
clear all
Array=1; %different fiber array (1)SEP (2)SDP (3)HP lo=10; %structure density
modal=7; %save number of modal
%================================================================
load MT.txt; %material properties
load element3d.txt; % nodal number of each element load nposition3d.txt; % nodal location
load dampinge.txt; %damping capacity of six direction load dampingd.txt;
load dampingh.txt;
if Array==1
damping=dampinge;
elseif Array==2
damping=dampingd;
else Array==3
damping=dampingh;
end
c=size(element3d);
d=size(nposition3d);
d=d(1,1);
G=zeros(3*3,8*3);
gk=zeros(d(1,1)*3,d(1,1)*3);
GK=zeros(d(1,1)*3,d(1,1)*3); % global K gkcapa=zeros(d(1,1)*3,d(1,1)*3);
GKcapa=zeros(d(1,1)*3,d(1,1)*3);
gm=zeros(d(1,1)*3,d(1,1)*3);
GM=zeros(d(1,1)*3,d(1,1)*3); % global M
%==================material property================
for i=Array
Ef1=MT(1,i); % zz Young's modulus of fiber Ef2=MT(2,i); % xx Young's modulus of fiber Ef3=Ef2; % yy Young's modulus of fiber nuf12=MT(5,i);
nuf13=nuf12;
nuf23=MT(6,i);
Gf12=MT(3,i); %Ef1/(2*(1+nuf12)) - if fibers are isotropic Gf13=Gf12;
Gf23=MT(4,i); %Ef2/(2*(1+nuf23));%Gf12;
end
nuf21=Ef2*nuf12/Ef1;
nuf31=Ef3*nuf13/Ef1;
nuf32=Ef3*nuf23/Ef2;
Sf=[1/Ef1 -nuf21/Ef2 -nuf31/Ef3 0 0 0; -nuf12/Ef1 1/Ef2 -nuf32/Ef3 0 0 0;...
-nuf13/Ef1 -nuf23/Ef2 1/Ef3 0 0 0; 0 0 0 1/Gf12 0 0; 0 0 0 0 1/Gf23 0; 0 0 0 0 0 1/Gf13];
E=inv(Sf);
%======================= B matrix========================
for a=1:c(1,1);
KE=zeros(24,24);
KM=zeros(24,24);
KEcapa=zeros(24,24);
for k=0:2:2 for j=0:2:2 for i=0:2:2
XI=(i-1)/sqrt(3);
YI=(j-1)/sqrt(3);
ZI=(k-1)/sqrt(3);
L=[1 0 0 0 0 0 0 0 0 % L matrix 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0];
DN(5,1)=-1/8*(1-YI)*(1+ZI);
DN(6,1)= 1/8*(1-YI)*(1+ZI);
DN(7,1)= 1/8*(1+YI)*(1+ZI);
DN(8,1)=-1/8*(1+YI)*(1+ZI);
DN(1,1)=-1/8*(1-YI)*(1-ZI);
DN(2,1)= 1/8*(1-YI)*(1-ZI);
DN(3,1)= 1/8*(1+YI)*(1-ZI);
DN(4,1)=-1/8*(1+YI)*(1-ZI);
DN(5,2)=-1/8*(1-XI)*(1+ZI);
DN(6,2)=-1/8*(1+XI)*(1+ZI);
DN(7,2)= 1/8*(1+XI)*(1+ZI);
DN(8,2)= 1/8*(1-XI)*(1+ZI);
DN(1,2)=-1/8*(1-XI)*(1-ZI);
DN(2,2)=-1/8*(1+XI)*(1-ZI);
DN(3,2)= 1/8*(1+XI)*(1-ZI);
DN(4,2)= 1/8*(1-XI)*(1-ZI);
DN(5,3)= 1/8*(1-XI)*(1-YI);
DN(6,3)= 1/8*(1+XI)*(1-YI);
DN(7,3)= 1/8*(1+XI)*(1+YI);
DN(8,3)= 1/8*(1-XI)*(1+YI);
DN(1,3)=-1/8*(1-XI)*(1-YI);
DN(2,3)=-1/8*(1+XI)*(1-YI);
DN(3,3)=-1/8*(1+XI)*(1+YI);
DN(4,3)=-1/8*(1-XI)*(1+YI);
NN(1,1)=1/8*(1-XI)*(1-YI)*(1-ZI);
NN(2,1)=1/8*(1+XI)*(1-YI)*(1-ZI);
NN(3,1)=1/8*(1+XI)*(1+YI)*(1-ZI);
NN(4,1)=1/8*(1-XI)*(1+YI)*(1-ZI);
NN(5,1)=1/8*(1-XI)*(1-YI)*(1+ZI);
NN(6,1)=1/8*(1+XI)*(1-YI)*(1+ZI);
NN(7,1)=1/8*(1+XI)*(1+YI)*(1+ZI);
NN(8,1)=1/8*(1-XI)*(1+YI)*(1+ZI);
index=0; % for I=1:8
index=index+3;
for J=1:3
G(J,index-2)=DN(I,J);
G(J+3,index-1)=DN(I,J);
G(J+6,index)=DN(I,J);
end
end
index2=0; % N(shape function) matrix for J=1:8
index2=index2+3;
NNN(1,index2-2)=NN(J,1);
NNN(2,index2-1)=NN(J,1);
NNN(3,index2)=NN(J,1);
end
%=============================Jacobian Matrix=========================
JACO=zeros(3,3);
for M=1:3 for N=1:3
for K=1:8
JACO(M,N)=JACO(M,N)+DN(K,M)*nposition3d(element3d(a,K+6),N+1);
end end end
INVJ=inv(JACO);
T=[INVJ(1,1) INVJ(1,2) INVJ(1,3) 0 0 0 0 0 0 INVJ(2,1) INVJ(2,2) INVJ(2,3) 0 0 0 0 0 0 INVJ(3,1) INVJ(3,2) INVJ(3,3) 0 0 0 0 0 0 0 0 0 INVJ(1,1) INVJ(1,2) INVJ(1,3) 0 0 0 0 0 0 INVJ(2,1) INVJ(2,2) INVJ(2,3) 0 0 0 0 0 0 INVJ(3,1) INVJ(3,2) INVJ(3,3) 0 0 0 0 0 0 0 0 0 INVJ(1,1) INVJ(1,2) INVJ(1,3) 0 0 0 0 0 0 INVJ(2,1) INVJ(2,2) INVJ(2,3) 0 0 0 0 0 0 INVJ(3,1) INVJ(3,2) INVJ(3,3)];
B=L*T*G; % B
%======================= K, M Kfi matrix=====================
kk=B'*E*B;
kkcapa=B'*E*damping*B;
KE=KE+det(JACO)*kk; % K matrix of element
KEcapa=KEcapa+det(JACO)*kkcapa; %kcapa matrix of element
km=lo*NNN'*NNN; %mass matrix of gauss integral point KM=KM+det(JACO)*km; %mass matrix of element
end end end
%======================Global K==========================
for n=1:8 for m=1:8
gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=KE(n*3-2,m*3-2);
GK(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=GK(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3-2)+gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2);
gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=KE(n*3-1,m*3-2);
GK(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=GK(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3-2)+gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2);
gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=KE(n*3,m*3-2);
GK(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=GK(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3-2)+gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2);
gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=KE(n*3-2,m*3-1);
GK(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=GK(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3-1)+gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1);
gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=KE(n*3-1,m*3-1);
GK(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=GK(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3-1)+gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1);
gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=KE(n*3,m*3-1);
GK(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=GK(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3-1)+gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1);
gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=KE(n*3-2,m*3);
GK(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=GK(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3)+gk(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3);
gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=KE(n*3-1,m*3);
GK(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=GK(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3)+gk(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3);
gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=KE(n*3,m*3);
GK(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=GK(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3)+gk(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3);
%===============================Global
Kcapa======================================
gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=KEcapa(n*3-2,m*3-2);
GKcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=GKcapa(element3d(a,n+1+5)*3-2,eleme nt3d(a,m+1+5)*3-2)+gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2);
gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=KEcapa(n*3-1,m*3-2);
GKcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=GKcapa(element3d(a,n+1+5)*3-1 ,e lement3d(a,m+1+5)*3-2)+gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2);
gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=KEcapa(n*3,m*3-2);
GKcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=GKcapa(element3d(a,n+1+5)*3,ele
ment3d(a,m+1+5)*3-2)+gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2);
gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=KEcapa(n*3-2,m*3-1);
GKcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=GKcapa(element3d(a,n+1+5)*3-2 ,e lement3d(a,m+1+5)*3-1)+gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1);
gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=KEcapa(n*3-1,m*3-1);
GKcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=GKcapa(element3d(a,n+1+5)*3-1 ,e lement3d(a,m+1+5)*3-1)+gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1);
gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=KEcapa(n*3,m*3-1);
GKcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=GKcapa(element3d(a,n+1+5)*3, ele ment3d(a,m+1+5)*3-1)+gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1);
gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=KEcapa(n*3-2,m*3);
GKcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=GKcapa(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3)+gkcapa(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3);
gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=KEcapa(n*3-1,m*3);
GKcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=GKcapa(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3)+gkcapa(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3);
gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=KEcapa(n*3,m*3);
GKcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=GKcapa(element3d(a,n+1+5)*3, element3d(a,m+1+5)*3)+gkcapa(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3);
%==========================Global mass matrix================================
gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=KM(n*3-2,m*3-2);
GM(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2)=GM(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3-2)+gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-2);
gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=KM(n*3-1,m*3-2);
GM(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2)=GM(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3-2)+gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-2);
gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=KM(n*3,m*3-2);
GM(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2)=GM(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3-2)+gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-2);
gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=KM(n*3-2,m*3-1);
GM(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1)=GM(element3d(a,n+1+5)*3-2,
element3d(a,m+1+5)*3-1)+gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3-1);
gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=KM(n*3-1,m*3-1);
GM(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1)=GM(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3-1)+gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3-1);
gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=KM(n*3,m*3-1);
GM(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1)=GM(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3-1)+gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3-1);
gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=KM(n*3-2,m*3);
GM(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3)=GM(element3d(a,n+1+5)*3-2, element3d(a,m+1+5)*3)+gm(element3d(a,n+1+5)*3-2,element3d(a,m+1+5)*3);
gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=KM(n*3-1,m*3);
GM(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3)=GM(element3d(a,n+1+5)*3-1, element3d(a,m+1+5)*3)+gm(element3d(a,n+1+5)*3-1,element3d(a,m+1+5)*3);
gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=KM(n*3,m*3);
GM(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3)=GM(element3d(a,n+1+5)*3,
element3d(a,m+1+5)*3)+gm(element3d(a,n+1+5)*3,element3d(a,m+1+5)*3);
end end
end
%============================ boundary conditions============================
load BC.txt; % nodal boundary condition at x direction x=size(BC);
z=z(1,1);
GKnew=zeros(3*d-(3*x),3*d-(3*x));
GKcapanew=zeros(3*d-(3*x),3*d-(3*x));
GMnew=zeros(3*d-(3*x),3*d-(3*x));
GK1=zeros(1,(3*d-(3*x))^2);
GKcapa1=zeros(1,(3*d-(3*x))^2);
GM1=zeros(1,(3*d-(3*x))^2);
pk=923456;
for i=1:x
NBC=BCx(i,1)*3-2;
GK(:,NBC)=pk;
GK(:,NBC+1)=pk;
GK(:,NBC+2)=pk;
GK(NBC,:)=pk;
GK(NBC+1,:)=pk;
GK(NBC+2,:)=pk;
end index=0;
for i=1:3*d for j=1:3*d if GK(i,j)==pk;
else
index=index+1;
GK1(1,index)=GK(i,j);
GM1(1,index)=GM(i,j);
GKcapa1(1,index)=GKcapa(i,j);
end end end
for i=1:(3*d-(3*x))
GKnew(i,:)=GK1(1,(i-1)*(3*d-(3*x))+1:i*(3*d-(3*x)));
GMnew(i,:)=GM1(1,(i-1)*(3*d-(3*x))+1:i*(3*d-(3*x)));
GKcapanew(i,:)=GKcapa1(1,(i-1)*(3*d-(3*x))+1:i*(3*d-(3*x)));
end
%===================eigenvalue & eigenvector=============================
[eigenvector,eigenvalue]=eig(GKnew,GMnew);
eigenvalue=eig(GKnew,GMnew);
Esort=sort(eigenvalue);
qq=size(Esort);
qq=qq(1,1);
FF=sqrt(Esort)/(2*pi);
index2=0;
for i=1:modal index=0;
es=Esort(i,1);
for j=1:qq
index=index+1;
ee=eigenvalue(j,1);
if ee==es %pp=index index2=index2+1;
evector(:,index2)=eigenvector(:,index);
end end end
Mdamp=zeros(modal,1);
for i=1:modal
Mdamp(i,1)=(1/2*evector(:,i)'*GKcapanew*evector(:,i))/(1/2*evector(:,i)'*GKnew*evector(:,i));
end yy=size(evector);
Dvector=zeros(yy(1,1)/3,yy(1,2)*3);
DDvector=zeros(d,modal*3);
for k=1:modal
mr=evector(:,k)'*GMnew*evector(:,k);
index3=0;
index2=0;
index=1;
for i=1:d
BBCx=BCx(index,1);
if index==x index=1;
index3=index3+1;
elseif nposition3d(i,1)==BBCx index=index+1;
index3=index3+1;
DDvector(index3,k*3-2)=nposition3d(i,2);
DDvector(index3,k*3-1)=nposition3d(i,3);
DDvector(index3,k*3)=nposition3d(i,4);
else
index2=index2+1;
index3=index3+1;
DDvector(index3,k*3-2)=(evector(index2*3-2,k)/sqrt(mr)+nposition3d(i,2));
DDvector(index3,k*3-1)=(evector(index2*3-1,k)/sqrt(mr)+nposition3d(i,3));
DDvector(index3,k*3)=(evector(index2*3,k)/sqrt(mr)+nposition3d(i,4));
end end end
Table 2.1 Mechanical properties and thermal properties of the ingredients of fiber composites used in the GMC analysis.
4.11 n
6.42×10-11 A
1.18×10-4 α(1/ oC)
0.25 ν23
0.37 0.2
ν12
5.5 G23(GPa)
27.6 G12(GPa)
14 E2(GPa)
3.4 234
E1(GPa)
Matrix Fiber
4.11 n
6.42×10-11 A
1.18×10-4 α(1/ oC)
0.25 ν23
0.37 0.2
ν12
5.5 G23(GPa)
27.6 G12(GPa)
14 E2(GPa)
3.4 234
E1(GPa)
Matrix Fiber
Table 3.1 Four parameters used in ANSYS to simulate the nonlinear behavior of matrix materials
620 48
2000 29
Parameter
b R∞
R0 K
620 48
2000 29
Parameter
b R∞
R0 K
Table 4.1 Mechanical properties and damping capacities of fiber and matrix used in GMC analysis [13]
0.06537
Table 4.2 Damping property of fiber composites with SEP packing obtained by using the GMC and FEM analysis
1.5
Table 4.3 Damping property of fiber composites with SDP packing obtained by using the GMC and FEM analysis
3.4
Table 4.4 Damping property of fiber composites with HP packing obtained by using the GMC and FEM analysis
3
Table 4.5 Fiber array effect on the first three modal damping capacities of composite rod with fiber extended in x-direction under free-free boundary condition
0.02163
Table 4.6 Fiber array effect on the first three modal damping capacities of composite rod with fiber extended in z-direction under free-free boundary condition
0.02008
Table 4.7 Fiber array effect on the first three modal damping capacities of composite plate with fiber extended in x-direction under free-free boundary condition
0.01911
Table 4.8 Fiber array effect on the first three modal damping capacities of composite rod with fiber extended in x-direction under clamp-free boundary condition
0.02167
Table 4.9 Fiber array effect on the first three modal damping capacities of composite rod with fiber extended in z-direction under clamp-free boundary condition
0.01638
Table 4.10 Fiber array effect on the first three modal damping capacities of composite plate with fiber extended in x-direction under one side clamped boundary condition
0.02015
Table 4.11 Fiber array effect on the first three modal damping capacities of composite plate with fiber extended in z-direction under one side clamped boundary
0.01777 0.01966
0.01641 Third mode
0.02004 0.02092
0.01932 Second mode
0.01538 0.01749
0.01367 First mode
HP SDP
SEP
0.01777 0.01966
0.01641 Third mode
0.02004 0.02092
0.01932 Second mode
0.01538 0.01749
0.01367 First mode
HP SDP
SEP
x
1RVE
x
3x
2Fig. 2.1 Representative volume element, (RVE)
=1
γ γ =2 γ =Nγ
=1 β
=2 β β =Nβ
h1
h2 Nβ
h
l1 l2 lNγ
l
h
x1
x
2x
3=1
γ γ =2 γ =Nγ
=1 β
=2 β β =Nβ
h1
h2 Nβ
h
l1 l2 lNγ
l
h
x1
x
2x
3Fig. 2.2 Geometry and coordinate system of representative volume element [1]
( )
βˆFig. 2.3 Local coordinate systems of the representative volume element [1]
Square edge array Square diagonal array Square hexagonal array
Fig. 2.4 Fiber composites with three different fiber arrangements
(a) (b)
Fig. 2.5 RVE with square edge packing portioned into (a) 28
×
28 subcells and (b) 39×
39 subcells(a) (b)
Fig. 2.6 RVE with square diagonal packing portioned into (a) 27
×
27 subcells and (b) 39×
39 subcells(a) (b)
Fig. 2.7 RVE with square edge packing portioned into (a) 20
×
35 subcells and (b) 31×
49 subcellsStrain
St re s s (M Pa )
0 0.1 0.2 0.3
0 50 100 150 200 250 300 350
Square edge (28x28) Square edge (39x39) Square diagonal (27x27) Square diagonal (39x39) Hexagonal (20x35) Hexagonal (31x49)
Fig. 2.8 Comparison of stress and strain curves obtained from the RVEs with coarse subclls and refined subcells
Strain
S tr e ss( M P a )
0 0.0025 0.005 0.0075 0.01
0 50 100 150 200 250 300
Square edge (Thermal) Square edge
Square diagonal (Thermal) Square diagonal
Hexagonal (Thermal) Hexagonal
Fig. 2.9 Thermal residual stress effects on the stress and strain curves of 15° off-axis fiber composites with three different fiber arrays. (fiber volume fraction 60%)
Strain
S tr e ss( M P a )
0 0.01 0.02
0 40 80 120 160 200
Square edge (Thermal) Square edge
Square diagonal (Thermal) Square diagonal
Hexagonal (Thermal) Hexagonal
Fig. 2.10 Thermal residual stress effects on the stress and strain curves of 30° off-axis fiber composites with three different fiber arrays. (fiber volume fraction 60%)
Strain
S tr e ss( M P a )
0 0.01 0.02 0.03 0.04
0 50 100 150 200
Square edge (Thermal) Square edge
Square diagonal (Thermal) Square diagonal
Hexagonal (Thermal) Hexagonal
Fig. 2.11 Thermal residual stress effects on the stress and strain curves of 45° off-axis fiber composites with three different fiber arrays. (fiber volume fraction 60%)
Strain
S tr e ss( M P a )
0 0.01 0.02 0.03 0.04
0 50 100 150 200
Square edge (Thermal) Square edge
Square diagonal (Thermal) Square diagonal
Hexagonal (Thermal) Hexagonal
Fig. 2.12 Thermal residual stress effects on the stress and strain curves of 60° off-axis fiber composites with three different fiber arrays. (fiber volume fraction 60%)
)
Fig. 3.1 Unit cell with 3
×
3 subcells to illustrate finite difference procedure(u)
Fig. 3.2 Coordinate system of the RVE with square edge packing fiber
Fig. 3.3 Finite element mesh of the RVE with square edge packing fiber
Fig. 3.4 Finite element mesh of the RVE with square diagonal packing fiber
Fig. 3.5 Finite element mesh of the RVE with hexagonal packing fiber
Fig. 3.6 Finite element mesh for a quadrant of the RVE with square edge packing fiber
Fig. 3.7 Finite element mesh for a quadrant of the RVE with square diagonal packing fiber
Fig. 3.8 Finite element mesh for a quadrant of the RVE with hexagonal packing fiber
(u)
x
1b
a c
(v) x
2(w)
x
3(u)
x
1b
a c
(v) x
2(w)
x
3Fig. 3.9 Coordinate system and dimension for a quadrant of the RVE with square edge packing fiber
Strain
St re s s (M Pa )
0 0.01 0.02 0.03 0.04 0.05 0.06
0 30 60 90 120
GMC FEM
Fig. 3.10 Stress strain curve of matrix employed in GMC, SCMC and FEM analysis
ε
22σ
220 0.005 0.01 0.015 0.02
0 50 100 150
Square edge (GMC) Square diagonal (GMC) Hexagonal (GMC) Square edge (FEM) Square diagonal (FEM) Hexagonal (FEM)
(MPa)
Fig. 3.11 Comparison of stress and strain curves obtained from GMC and FEM analysis under transverse loading σ22
ε
22σ
220 0.005 0.01 0.015 0.02
0 50 100 150
Square edge (SCMC) Square diagonal (SCMC) Hexagonal (SCMC) Square edge (FEM) Square diagonal (FEM) Hexagonal (FEM)
(MPa)
Fig. 3.12 Comparison of stress and strain curves obtained from SCMC and FEM analysis under transverse loading σ22
γ
12σ
120 0.005 0.01 0.015 0.02
0 20 40 60
Hexagonal (16x28) Hexagonal (26x45)
(MPa)
Fig. 3.13 Converge test for the stress and strain curves of composites with hexagonal square edge packing fiber obtained from SCMC model under shear loading σ12
γ
12σ
120 0.005 0.01 0.015 0.02
0 20 40 60
Hexagonal (16x28) Hexagonal (26x45)
(MPa)
Fig. 3.14 Converge test for the stress and strain curves of composites with hexagonal square edge packing fiber obtained from GMC model under shear loading σ12
Representatively volume
element Equivalent element E, Ψ
Composite structures Representatively volume
element Equivalent element E, Ψ
Composite structures
Fig. 4.1 Modeling procedure for characterizing the damping properties of composite structures.
10cm 2cm z x
y
2cm 10cm 2cm
z x y
z x y
2cm (a)
z x y
10cm
2cm
10cm z x
y
z x y
10cm
2cm
10cm
(b)
Fig. 4.2 The dimension of composite structures (a) composite rod with fiber extended in x-direction (b) composite plate with fiber extended in x direction