The thermal stress effect, fiber shape effect and fiber arrangement effect on nonlinear rate dependent behavior of fiber composites was investigated using square fiber model, generalized method of cells and finite element analysis. It was shown that with the presence of thermal residual stress, the fiber composites become stiffer than those without thermal stress. Regarding the fiber arrangement effects, the RVE with square edge packing (SEP) array is stiffer than that with square diagonal packing (SDP). In addition, the RVEs containing either round fiber or square fiber result in the similar material behavior of fiber composites.
In order to validate the model predictions, the experimental results of the graphite/epoxy composites tested at various strain rates were employed for the comparison. The model predictions were obtained from the SFM and GMC analysis on the RVE with SEP fiber array containing the transverse isotropic fiber and elastic-viscoplastic matrix. The viscoplastic behaviors of epoxy resin were described using a three parameters viscoplasticity model written in the form of power law. Comparing the model predictions and experimental data reveals that the GMC is better than the SFM to predict the nonlinear rate sensitivity of off-axis specimens at strain rate up to 550/s if the current material properties were applied. Comparison of model predictions obtained from GMC and SFM analysis with the experimental results revealed that the micromechanical approaches are capable of predicting the nonlinear rate sensitivity of off-axis specimens although there are still distinctions between the model and the experimental results.
Reference
[1] Sun, C. T. and Chen, J. L. 1989 “A Simple Flow Rule for Characterizing Nonlinear Behavior of Fiber Composites,” Journal of Composite Materials, Vol.
23, No. 10, pp. 1009-1020.
[2] Gates, T. S. and Sun, C. T. 1991 “Elastic/Viscoplastic Constitutive Model for Fiber Reinforced Thermoplastic Composites,” AIAA Journal, Vol. 29, No. 3, pp.
457-463.
[3] Eisenberg, M. A. and Yen, C. F. 1981 “Theory of Multiaxial Anisotropic Viscoplasticity,” ASME Journal of Applied Mechanics, Vol. 48, No.2, pp.
276-284.
[4] Yoon, K. J. and Sun, C. T. 1991 “Characterization of Elastic-Viscoplastic Properties of an AS4/PEEK Thermoplastic Composite,” Journal of Composite Materials, Vol. 25, No. 10, pp. 1277-1296.
[5] Stouffer, D. C. and Dame, L. T. 1996 Inelastic Deformation of Metals/Models, Mechanical Properies, and Metallurgy, John Wiley & Sons, Inc., New York.
[6] Weeks, C. A. and Sun, C. T. 1998 “Modeling Non-Linear Rate-Dependent Behavior in Fiber-Reinforced Composites,” Composites Science and Technology, Vol. 58, No. 3-4, pp. 603-611.
[7] Meyers, M. A. 1994 Dynamic Behavior of Materials, John Wiley & Sons, New York.
[8] Thiruppukuzhi, S. V. and Sun, C. T. 2001 “Models for the Strain-Rate-Dependent Behavior of Polymer Composites,” Composites Science and Technology, Vol. 61, No. 1, pp. 1-12.
[9] Eshelby, J. D. 1957 “The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems,” Proceedings of the Royal Society of London, A241, No. 1226, pp. 376-396.
[10] Mori, T. and Tanaka, K. 1973 “Average Stress in Matrix and Average Elastic Energy of Materials with Misfitting Inclusions,” Acta Metallurgica, Vol. 21, No. 5, pp. 571-574.
[11] Benveniste, Y. 1987 “A New Approach to the Application of Mori-Tanaka Theory in Composite Materials,” Mechanics of Materials, Vol. 6, No. 2, pp.
147-157.
[12] Sun, C. T. and Chen, J. L. 1991 “A Micromechanical Model for Plastic Behavior of Fibrous Composites,” Composite Science and Technology, Vol. 40, No. 2, pp.
115-129.
[13] Aboudi, J. 1982 “A Continuum Theory for Fiber-Reinforced Elastic-Viscoplastic Composites,” International Journal of Engineering Science, Vol. 20, No. 5, pp.
605-621.
[14] Aboudi, J. 1991 Mechanics of Composite Materials – A Unified Micromechanical Approach, Elsevier Science Publishing Company, Inc., New York, pp. 35-53.
[15] 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.
[16] Mura, T. 1987 Micromechanics of Defects in Solids, Kluwer Academic, Dordrecht.
[17] Pettermann, H. E. and Plankensteiner, A. F. and Bohm, H. J. and Rammerstorfer, F. G. 1999 “A Thermo-Elastic-Plastic Constitutive Law for Inhomogeneous Materials Based on an Incremental Mori-Tanaka Approach, ” Computers &
Structures, Vol. 71, pp. 197-214.
[18] Goldberg, R. K. and Stouffer, D. C. 2002 “Strain Rate Dependent Analysis of a Polymer Matrix Composite Utilizing a Micromechanics Approach,” Journal of Composite Materials, Vol. 36, No. 7, pp. 773-793.
[19] Orozco, C. E. 1997 “Computational Aspects of Modeling Complex Microstructure Composites Using GMC,” Composites Part B, Vol. 28, No. 1-2, pp.
167-175.
[20] 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.
[21] 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.
[22] Ogihara, S., Kobayashi, S. and Reifsnider, K. L. 2003 “Characterization of Nonlinear Behavior of Carbon/Epoxy Unidirectional and Angle-Ply Laminates,”
Advanced Composite Materials, Vol. 11, No. 3, pp. 239-254.
[23] Kawai, M., Masuko, Y., Kawase, Y. and Negishi, R. 2001 “Micromechanical Analysis of the Off Axis Rate-Dependent Inelastic Behavior of Unidirectional AS4/PEEK at High Temperature,” International Journal of Mechanical Sciences, Vol. 43, No. 9, pp. 2069-2090.
[24] 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.
[25] Tuttle, M. E. 1989 “Fundamental Strain-Gage Technology,” Manual on Experiment Methods for Mechanical Testing of Composites, edited by Pendleton, R. L. and Tuttle, M. E., Society for Experimental Mechanics, Bethel, CT.
[26] 1987 “Measurement of Thermal Expansion Coefficient Using Strain Gages,”
Micro-Measurement Technical Note TN-513, Measurement Group, Inc.
[27] Scalea, F. L. D. 1998 “Measurement of Thermal Expansion Coefficient of Composites Using Strain Gages,” Experimental Mechanics, Vol. 38, No. 4, pp.
233-241.
[28] ANSYS Element Reference. 000855. Eighth Edition. SAS IP, Inc. 1997.
[29] ANSYS Theory Reference. 000855. Eighth Edition. SAS IP, Inc. 1997.
[30] 郭濬清, “應變率對纖維複合材料側向抗壓強度影響之探討,” 國立交通大學 碩士論文, 2005.
[31] Tsai, S. W. and Hahn, H. T. 1980 Introduction to Composite Materials, Technomic Publishing Co., Westport, CT, USA.
Appendix A. A MATLAB Code for Square Fiber Model
Input Symbol List
Ef1, Ef2, nuf12, Gf12, volf – Young’s modulus 1 (and 2), Poisson’s ratio ,shear modulus and volume fraction of the fiber.
Em, num, Gm – Young’s modulus, Poisson’s ratio and shear modulus of the matrix.
pn, xi, pm – denote n, χ and m in the three parameters model.
time – A matrix to store the total loading time. The Row represents different strain rates while the column denotes various off-axis angles.
ratetype – An index to point what strain rate was in calculation.
angle, Nstep – off-axis angle and total steps of loop.
Matrix and Scalar Symbol List saf – Fiber compliance matrix.
eSm – Matrix compliance matrix.
v1, v2 – Volume fraction of region AF and AM relative to region A.
Tsig, Teps – Transformation matrix relatively for stress and strain to translate between loading coordinates and material principle axis.
sigma, sigb – Current stress states in region AM and B, respectively.
Tstress, Tstrain – Overall stress and strain in the loading coordinate system.
depsam, depsaf – Total strain rates in region AM and AF.
dpepsam, dpepsb – Plastic strain rates in region AM and B.
dsigaf, dsigam – Stress rate in region AF and AM.
t, sigrate – Current time and associated stress rate.
dsig, dsig1 – Applied stress rate in the material principle axis and loading coordinate system, respectively.
efsig, efeps – Effective stress and effective plastic strain rate.
pSm, sa, sb, ca, cb, c, s -
[ ]
, Sp[ ]
, SA[ ]
, SB[ ]
, CA[ ]
, CB[ ]
C ,[ ]
Sdeps, deps1 – Strain rate associated to applied stress rate in the material coordinate system and loading coordinate system.
Note: Subroutine “stressrate(t,ratetype,angle/15)”, which involves stress rate – time polynomial equation from curve-fitting of experimental data, have three input parameters to output corresponding stress rate at each time increment. The subroutine can be designed by users and was not demonstrated here.
Code
% Square Fiber Model Combined with Three Parameters Model clear all
%---Material Properties---
% Fiber
Ef1=235000; % E11 of fiber Ef2=18000; % E22 of fier
nuf12=0.2; % Poisson’s ratio 12 of fiber Gf12=35000; % Shear modulus 12 of fiber volf=0.69; % fiber volume fraction
% Matrix
Em=3400; % Young’s modulus of the matrix num=0.373; % Poisson’s ratio of the matrix Gm=Em/(2*(1+num)); % Shear modulus of the matrix
pn=5.62; % Coefficient n in the 3 parameters model xi=1.23e-14; % Coefficient xi in the 3 parameters model pm=-0.168; % Coefficient m in the 3 parameters model
%---
% Loop & angle
time=[208 694 695 675 380;
1.83 7.25 10.9 6.9 4.54;
0.0323 0.0725 0.0981 0.0801 0.05;
5.1e-5 7.07e-5 7.84e-5 8.65e-5 6.4e-5];
ratetype=3;
angle=30; % Off-axis angle if angle == 90
totime=time(ratetype,5);
else
totime=time(ratetype,angle/15);
end
%totime=time(ratetype,2);
Nstep=450; % Total time steps
dtime=totime/Nstep; % Associated time increment rad=angle/180*pi;
%---
saf=[1/Ef1 -nuf12/Ef1 0; -nuf12/Ef1 1/Ef2 0; 0 0 1/Gf12]; % fiber compliance matrix eSm=[1/Em -num/Em 0; -num/Em 1/Em 0; 0 0 1/Gm]; % matrix compliance matrix (elastic)
v1=sqrt(volf); % volume fraction of AF with respective to region A v2=1-v1; % volume fraction of AM with respective to region A
% Transformation matrix to translate stress
Tsig=[cos(rad)^2 sin(rad)^2 2*cos(rad)*sin(rad);...
sin(rad)^2 cos(rad)^2 -2*cos(rad)*sin(rad);...
-cos(rad)*sin(rad) cos(rad)*sin(rad) cos(rad)^2-sin(rad)^2];
% Transformation matrix to translate strain
Teps=[cos(rad)^2 sin(rad)^2 -cos(rad)*sin(rad);...
sin(rad)^2 cos(rad)^2 cos(rad)*sin(rad);...
2*cos(rad)*sin(rad) -2*cos(rad)*sin(rad) cos(rad)^2-sin(rad)^2];
% ---matrix pre-setup--- sigam=zeros(3,1);
sigb=zeros(3,1);
Tstress=zeros(Nstep+1,1);
Tstrain=zeros(Nstep+1,1);
depsam=zeros(3,1);
depsaf=zeros(3,1);
dpepsam=zeros(3,1);
dpepsb=zeros(3,1);
depsb=zeros(3,1);
dsigaf=zeros(3,1);
dsigam=zeros(3,1);
%---Main Loop--- for n=1:Nstep;
t=n*dtime; % current time sigrate=strainrate(t,ratetype,angle/15);
dsig1=[sigrate; 0; 0]; % x-y stress rate Tstress(n+1,1)=Tstress(n,1)+dsig1(1,1)*dtime; % Overall stress dsig=Tsig*dsig1; % 1-2 stress rate
efeps=2/3*sqrt(3*(dpepsam(1)^2+dpepsam(2)^2+dpepsam(1)*dpepsam(2)+dpepsam(3)^2/
4));
if n = = 1 A=0;
else
A=xi*(efeps)^pm;
end
s1=1/3*(2*sigam(1,1)-sigam(2,1));
s2=1/3*(-sigam(1,1)+2*sigam(2,1));
s3=2*sigam(3,1);
efsig=sqrt((sigam(1,1)+sigam(2,1))^2-3*(sigam(1,1)*sigam(2,1)-sigam(3,1)^2));
st=[2*(sigam(1,1)+sigam(2,1))-3*sigam(2,1) 2*(sigam(1,1)+sigam(2,1))-3*sigam(1,1) 6*sigam(3,1)];
pSm=3/4*A*pn*efsig^(pn-3)*[s1*st; s2*st; s3*st];
sam=eSm+pSm;
% sm matrix
b1=(v1*Ef1*sam(1,1)+v2)^-1;
b2=-v1*(Ef1*sam(1,2)+nuf12)*(v1*Ef1*sam(1,1)+v2)^-1;
b3=-v1*Ef1*sam(1,3)*(v1*Ef1*sam(1,1)+v2)^-1;
a1=(1-v2*b1)/v1;
a2=-v2*b2/v1;
a3=-v2*b3/v1;
sa(1,1)=a1/Ef1;
sa(1,2)=a2/Ef1-nuf12/Ef1;
sa(1,3)=a3/Ef1;
sa(2,1)=sa(1,2);
sa(2,2)=v1*(1/Ef2-nuf12*a2/Ef1)+v2*(b2*sam(2,1)+sam(2,2));
sa(2,3)=-v1*a3*nuf12/Ef1+v2*(b3*sam(2,1)+sam(2,3));
sa(3,1)=sa(1,3);
sa(3,2)=sa(2,3);
sa(3,3)=v1/Gf12+v2*(b3*sam(3,1)+sam(3,3));
ca=inv(sa);
% CB matrix
efeps=2/3*sqrt(3*(dpepsb(1)^2+dpepsb(2)^2+dpepsb(1)*dpepsb(2)+dpepsb(3)^2/4));
if n = = 1 A=0;
else
A=xi*(efeps)^pm;
end
s1=1/3*(2*sigb(1,1)-sigb(2,1));
s2=1/3*(-sigb(1,1)+2*sigb(2,1));
s3=2*sigb(3,1);
efsig=sqrt((sigb(1,1)+sigb(2,1))^2-3*(sigb(1,1)*sigb(2,1)-sigb(3,1)^2));
st=[2*(sigb(1,1)+sigb(2,1))-3*sigb(2,1) 2*(sigb(1,1)+sigb(2,1))-3*sigb(1,1) 6*sigb(3,1)];
pSm=3/4*A*pn*efsig^(pn-3)*[s1*st; s2*st; s3*st];
sb=eSm+pSm;
cb=inv(sb);
c=v1*ca+v2*cb;
s=inv(c);
deps=s*dsig; % total strain rate
depsa=deps; % total strain rate of region A depsb=deps; % total strain rate of region B dsiga=ca*depsa; % stress rate of region A dsigb=cb*depsb; % stress rate of region B dsigaf(2:3,1)=dsiga(2:3,1);
dsigam(2:3,1)=dsiga(2:3,1);
depsaf(1,1)=depsa(1,1);
dsigaf(1,1)=depsaf(1,1)*Ef1+nuf12*dsigaf(2,1);
depsaf(2,1)=saf(2,:)*dsigaf;
depsaf(3,1)=saf(3,:)*dsigaf;
depsam(1,1)=depsa(1,1);
depsam(2,1)=(depsa(2,1)-v1*depsaf(2,1))/v2;
depsam(3,1)=(depsa(3,1)-v1*depsaf(3,1))/v2;
dsigam(1,1)=(dsiga(1,1)-v1*dsigaf(1,1))/v2;
dpepsam=depsam-Sm*dsigam;
dpepsb=depsb-Sm*dsigb;
deps1=Teps*deps;
Tstrain(n+1,1)=Tstrain(n,1)+deps1(1,1)*dtime;
sigam=sigam+dsigam*dtime;
sigb=sigb+dsigb*dtime;
end
B. A MATLAB Code for Generalized Method of Cells
Since some symbols in the GMC code are a repetition of those appearing in the SFM code, only differences are listed. This code is addressed on four regions RVE with fiber phase at β=γ=1.
Input Symbol List h, l – The height and the width of the RVE.
Nb, Nr - N , β N γ
Nfiber – Denote how many subcells are occupied by the fiber phase.
regionf – A matrix to store what cells are occupied by fibers. The fibrous subcell )
(βγ must be translated to a re-defined fiber number by using the equation γ
+
− β
=( 1)Nγ
N and N was arranged in the regionf matrix.
Matrix and Scalar Symbol List bb, hh – The length of subcell
1 1
= γ
= β and
2 1
= γ
=
β in the x2 direction.
hb, hr – Two matrix to store the geometry information bb and hh.
Sm – Compliance matrix including plastic part.
AG, J, AVPM, AVP, K, AtVP, Bvp -
[ ]
A , G[ ]
J ,[ ]
AVPM ,[ ]
AVP ,[ ]
K ,[ ]
A~VP , B*VPsubstress – A matrix to store current stress state of all subcells.
sigma – A transition matrix to get current stress state for specified subcell from substress matrix.
Cs – Stiffness matrix of all subcells.
Sf, Cf – Compliance matrix and stiffness matrix of the fiber.
eSm, eCm – Elastic compliance matrix and elastic stiffness matrix of the matrix.
index, index2 – Indexes to point out which subcell should be under calculation.
betahat, gamahat - βˆ and γˆ .
subteps, subsigrate – Strain rate and stress rate of specified subcell.
subeeps, subeps – Elastic strain rate and plastic strain rate of specified subcell.
Code
clear all
% Material Properties (transverse isotropic fiber + isotropic matrix) Ef1=235000; % Young’s modulus 11 of the fiber Ef2=18000; % Young’s modulus 22 of the fiber Ef3=Ef2; % Young’s modulus 33 of the fiber nuf12=0.2; % Poisson’s ratio 12 of the fiber nuf13=0.2; % Poisson’s ratio 13 of the fiber nuf23=0.25; % Poisson’s ratio 23 of the fiber Gf12=35000; % Shear modulus 12 of the fiber Gf13=Gf12; % Shear modulus 13 of the fiber Gf23=Ef2/(2*(1+nuf23)); % Shear modulus 23 of the fiber Em=3400; % Young’s modulus of the matrix num=0.373; % Poisson’s ratio of the matrix Gm=Em/(2*(1+num)); % Shear modulus of the matrix
% Coefficients of 3 parameter model
pn=5.62; % Coefficient n in the 3 parameters model xi=1.23e-14; % Coefficient xi in the 3 parameters model pm=-0.168; % Coefficient m in the 3 parameters model
% region dimension & cell situation
h=1; % the length of the RVE
l=1; % the width of the RVE
volf=0.69; % fiber volume fraction Nb=2;
Nr=2;
bb=sqrt(volf);
hh=1-bb;
hb=[bb hh];
hr=hb;
Nfiber=1;
regionf=[1];
% Loop & angle
time=[208 694 695 675 380;
1.83 7.25 10.9 6.9 4.54;
0.0323 0.0725 0.0981 0.0801 0.05;
5.1e-5 7.07e-5 7.84e-5 8.65e-5 6.4e-5];
ratetype=3;
angle=60; % Off-axis angle
if angle == 90
totime=time(ratetype,5);
else
totime=time(ratetype,angle/15);
end
Nstep=500; % Total time steps
dtime=totime/Nstep; % Time increment
rad=angle/180*pi;
% matrix pre-setup Sm=zeros(6,6);
AG=zeros(2*(Nb+Nr)+Nb*Nr+1,6*Nb*Nr);
J=zeros(2*(Nb+Nr)+Nb*Nr+1,6);
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);
Cs=zeros(6*Nb*Nr,6);
Tstrain=zeros(Nstep+1,1);
Tstress=zeros(Nstep+1,1);
AVP=zeros(6*Nb*Nr,6);
% compliance matrix of the fiber and the 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);
eSm=[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];
eCm=inv(eSm);
% initial Cs for beta=1:Nb for gama=1:Nr
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
multi=1;
region=regionf-index;
for j=1:Nfiber
multi=multi*region(1,j);
end
if multi ~= 0 Cs(p1:p2,:)=eCm;
else
Cs(p1:p2,:)=Cf;
end end
end
% 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 for beta=1:Nb
index=(beta-1)*Nr+gama;
AG(Nb*Nr+gama,(index-1)*6+2)=hb(1,beta);
J(Nb*Nr+gama,2)=h;
AG(Nb*Nr+Nr+beta,(index-1)*6+3)=hr(1,gama);
J(Nb*Nr+Nr+beta,3)=l;
AG(Nb*Nr+Nr+Nb+1,(index-1)*6+4)=hb(1,beta)*hr(1,gama)/2;
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;
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;
J(Nb*Nr+Nr+2*Nb+1+gama,6)=h/2;
end end
% Tsig & Teps (Transformation Matrix)
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];
% Main Loop for i=1:Nstep t=i*dtime;
sigrate=strainrate(t,ratetype,angle/15);
dsig1=[sigrate; 0; 0; 0; 0; 0];
dsig=Tsig*dsig1;
Tstress(i+1,1)=Tstress(i,1)+dsig1(1,1)*dtime;
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);
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);
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);
AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end end
for gama=1:Nr-1 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);
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);
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);
AVPM(count,6*index2+j)=-Cs(6*index2+4,j);
end end
AtVP=[AVPM; AG];
AVP=inv(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=inv(Bvp)*dsig;
deps1=Teps*deps;
Tstrain(i+1,1)=Tstrain(i,1)+deps1(1,1)*dtime;
for beta=1:Nb for gama=1:Nr
index=(beta-1)*Nr+gama;
p1=6*(index-1)+1;
p2=6*(index-1)+6;
multi=1;
region=regionf-index;
for j=1:Nfiber
multi=multi*region(1,j);
end
if multi ~= 0
subteps=AVP(p1:p2,:)*deps;
subsigrate=Cs(p1:p2,:)*subteps;
subeeps=eSm*subsigrate;
subeps=subteps-subeeps;
efeps=2/3*sqrt(1/2*((subeps(1,1)-subeps(2,1))^2+(subeps(2,1)-subeps(3,1 ))^2+(subeps(3,1)-subeps(1,1))^2)+3/4*(subeps(4,1)^2+subeps(5,1)^2+sube ps(6,1)^2));
substress(p1:p2,1)=substress(p1:p2,1)+subsigrate*dtime;
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));
A=xi*(efeps)^pm;
pSm=9/4*A*pn*efsig^(pn-3)*[s1; s2; s3; s4; s5; s6]*[s1; s2; s3; s4; s5;
s6]';
Sm=eSm+pSm;
Cs(p1:p2,:)=inv(Sm);
else
Cs(p1:p2,:)=Cf;
end end
end end
Table 1. Material properties used in the micromechanical analysis where the matrix properties were obtained from experiments and the fiber properties were determined to fit the linear elastic of experimental data at all off-axis angles.
(GPa)
Table 2. Material properties employed in the finite element analysis.
)
10 mm
=12 mm
φ
175 mm 75 mm
4.3 mm 20 mm
R=60 mm
10 mm 10 mm
=12 mm
φ
175 mm 75 mm
4.3 mm 20 mm
R=60 mm
10 mm
(a) (b)
Fig. 2.1 Dimensions of tensile and compression specimens. (a) Cylindrical compression specimen. (b) Coupon tensile specimen.
Specimen Self-Adjusting Device
Fig. 2.2 Experimental setup for compression tests.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0
20 40 60 80 100 120 140 160
1/s 0.01/s 0.0001/s
Fig. 2.3 Compression test results of polymer at 10-4, 10-2 and 1/s strain rates.
Axial Strain
LateralStrain
0 0.005 0.01 0.015 0.02
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Tensile Test
ε11 ε11 ε 22(%)
Fig. 2.4 Tensile test result to determine Poisson’s ratio of the polymer.
ε
aε
tε
aε
tFig. 2.5 Schematic for a strain gage subjected to a biaxial strain field.
∆E
V
Reference material Test
material
R1
R3 R4
R2
∆E
V
Reference material Test
material
R1
R3 R4
R2
Fig. 2.6 A half-bridge circuit for measuring the coefficient of thermal expansion.
Temperature (oC) Strain(10-3 )
0 20 40 60 80 100 120 140
-4 -3 -2 -1 0 1 2 3 4 5
entire process
curve-fitting (heating) curve-fitting (cooling)
y=0.058x-1.71
y=0.061x-3.43
Fig. 2.7 Thermal response for the polymer.
Effective Plastic Strain
EffectiveStress(MPa)
0 0.005 0.01 0.015 0.02
0 20 40 60 80 100 120 140 160
1/s 0.01/s 0.0001/s 1/s curve-fitting 0.01/s curve-fitting 0.0001/s curve-fitting
A=6.70E-14 (MPa)-n A=2.95E-14 (MPa)-n A=1.42E-14 (MPa)-n
Fig. 2.8 Effective stress – effective plastic strain curves for polymer at three different strain rates.
Time (s)
EffectivePlasticStrain(%)
0 100 200 300 400 500 600 700 0
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
0.0001/s
Fig. 2.9 Effective plastic strain versus time curve for epoxy at strain rate of 10-4/s.
Log(Effective Plastic Strain Rate)
Log(A)
10-5 10-4 10-3 10-2 10-1 100 1E-14
1E-13
0.0001/s 0.01/s 1/s
curve-fitting
Fig. 2.10 Log-log plot for determining the parameters in the viscoplasticity model.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04 0.05 0.06 0
30 60 90 120 150 180
650/s (SHPB)
Fig. 2.11 The stress – strain curve of the polymer from SHPB results.
Time
Stress(MPa)
0 10 20 30 40 50 60 70 80 90 0
30 60 90 120 150 180
SHPB Results
µs) (
Fig. 2.12 The stress – time curve from SHPB test.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04 0.05 0.06 0
30 60 90 120 150 180
650/s (Exp) 1/s (Exp) 10-2/s (Exp) 10-4/s (Exp) 650/s (Prediction) 1/s (Prediction) 10-2/s (Prediction) 10-4/s (Prediction)
Fig. 2.13 Prediction results of polymer at different strain rates by using three parameters model.
B
AF AM
x
2x
1x
3h
1h
2h
3h
4B
AF AM
x
2x
1x
3B
AF AM B
AF AM
x
2x
1x
3h
1h
2h
3h
4(a) (b) Fig. 3.1 Demonstration of RVEs with various fiber arrangements. (a) Square edge
packing array. (b) Square diagonal packing array.
Fig. 3.2 Geometry of square fiber model [12].
B
Fig. 3.3 Square edge packing array for modified square fiber model.
A
Fig. 3.4 Fiber distribution of square diagonal packing array (SDP) based on the fiber volume fraction. (a) Less than 39.3 %. (b) Equal to 39.3 %. (c) Greater than 39.3
%. (d) Attain to maximum fiber volume fraction 78.5 %.
AF1 AM1
Fig. 3.5 Square diagonal packing array for modified square fiber model.
BFL
nBM
nBFR
nFig. 3.6 Two fiber phases can be treated as a whole one if the constant stress or strain assumptions in eqn (3.2.21) were applied.
=1
Fig. 4.1 The coordinate system and geometry information of the generalized method of cells [15].
Fig. 4.2 Local coordinate systems of the generalized method of cells [15].
( β , γ )
( ) β ˆ , γ
( β , γ ˆ )
X
2X
3( ) β ˆ , γ ˆ
( 0 , 1 , 0 ) ( 0 , 0 , 1 ) ( β , γ )
( ) β ˆ , γ
( β , γ ˆ )
X
2X
3( ) β ˆ , γ ˆ
( 0 , 1 , 0 ) ( 0 , 0 , 1 )
Fig. 4.3 Normal vectors at the interfaces of subcells.
= 1 β
= 2 β
= 1
γ γ = 2
h
1h
2h
1h
2x 2
x 3
= 1 β
= 2 β
= 1
γ γ = 2
h
1h
2h
1h
2x 2
x 3
Fig. 4.4 A four regions RVE employed in the GMC, in which β=γ=1 represents the fiber phase.
X
1X
2X
3σ
11σ
22σ
12a
a a
X
1X
2X
3σ
11σ
22σ
12a
a a
Fig. 5.1 3-D square diagonal packing array employed in ANSYS.
Fig. 5.2 (a) A finite element mesh generated by ANSYS.
Fig. 5.2 (b) A full view of finite element mesh.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04 0.05 0
20 40 60 80 100 120 140
Assumed matrix
Fig. 5.3 An assumed stress–strain curve of the matrix.
x
z
y
Fiber Matrix
δ l
x y
Fiber
Matrix EmαmAm
f f
f A
E α
(a) (b)
x
z
y
Fiber Matrix
x
z
y x
z
y
Fiber Matrix
δ l
x y
Fiber
Matrix EmαmAm
f f
f A
E α
(a) (b)
Fig. 6.1 (a) Simplified model for unidirectional fiber composites. (b) Evaluation of thermal residual stress based on the displacement continuity in the x direction.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04
0 50 100 150 200 250
SFM (Thermal) SFM
Fig. 6.2(a) Thermal stress effect on the stress and strain curve of 300 fiber composite obtained from the square fiber model.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04 0.05 0.06 0
60 120 180 240
SFM (Thermal) SFM
Fig. 6.2(b) Thermal stress effect on the stress and strain curve of 900 fiber composite obtained from the square fiber model.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04
0 50 100 150 200 250
GMC (Thermal) GMC
Fig. 6.3(a) Thermal stress effect on the stress and strain curve of 300 fiber composite obtained from the generalized method of cells.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04
0 50 100 150 200 250 300
GMC (Thermal) GMC
Fig. 6.3(b) Thermal stress effect on the stress and strain curve of 900 fiber composite obtained from the generalized method of cells.
Fig. 6.4 The RVE with 26×26 subcells employed in the calculation of generalized method of cells (square edge packing).
B
A B
A
Fig. 6.5 The RVE with 50 subcells in fibrous region employed in the modified square fiber model (square edge packing).
Strain
Stress(MPa)
0 0.002 0.004 0.006 0.008 0.01 0
50 100 150 200 250 300
Square Fiber (GMC) Round Fiber (GMC) Square Fiber (SFM) Round Fiber (SFM)
Fig. 6.6(a) Fiber shape effects on the stress and strain curves of 150 fiber composites using the generalized method of cells (GMC) and the square fiber model (SFM).
Strain
Stress(MPa)
0 0.005 0.01 0.015 0.02
0 50 100 150 200
Square Fiber (GMC) Round Fiber (GMC) Square Fiber (SFM) Round Fiber (SFM)
Fig. 6.6(b) Fiber shape effects on the stress and strain curves of 300 fiber composites using the generalized method of cells (GMC) and the square fiber model (SFM).
Strain
Stress(MPa)
0 0.01 0.02 0.03
0 30 60 90 120 150 180
Square Fiber (GMC) Round Fiber (GMC) Square Fiber (SFM) Round Fiber (SFM)
Fig. 6.6(c) Fiber shape effects on the stress and strain curves of 450 fiber composites using the generalized method of cells (GMC) and the square fiber model (SFM).
Strain
Stress(MPa)
0 0.01 0.02 0.03
0 30 60 90 120 150 180
Square Fiber (GMC) Round Fiber (GMC) Square Fiber (SFM) Round Fiber (SFM)
Fig. 6.6(d) Fiber shape effects on the stress and strain curves of 600 fiber composites using the generalized method of cells (GMC) and the square fiber model (SFM).
Fig. 6.7 The RVE with 20×20 subcells employed in the calculation of generalized method of cells (square diagonal packing).
A
A B A
A B A
A B
Fig. 6.8 The RVE employed in the modified square fiber model (square diagonal packing).
Strain
Stress(MPa)
0 0.002 0.004 0.006 0.008
0 50 100 150 200 250 300
SEP (GMC) SDP (GMC) SEP (SFM) SDP (SFM)
Fig. 6.9(a) The effect of fiber arrangements on the stress and strain curves of 150 fiber composites obtained from the SFM and GMC.
Strain
Stress(MPa)
0 0.005 0.01 0.015 0.02 0.025 0
50 100 150 200 250
SEP (GMC) SDP (GMC) SEP (SFM) SDP (SFM)
Fig. 6.9(b) The effect of fiber arrangements on the stress and strain curves of 300 fiber composites obtained from the SFM and GMC.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04
0 40 80 120 160 200 240
SEP (GMC) SDP (GMC) SEP (SFM) SDP (SFM)
Fig. 6.9(c) The effect of fiber arrangements on the stress and strain curves of 450 fiber composites obtained from the SFM and GMC.
Strain
Stress(MPa)
0 0.01 0.02 0.03 0.04
0 50 100 150 200
SEP (GMC) SDP (GMC) SEP (SFM) SDP (SFM)
Fig. 6.9(d) The effect of fiber arrangements on the stress and strain curves of 600 fiber composites obtained from the SFM and GMC.
Strain
Stress(MPa)
0 0.002 0.004 0.006 0.008
0 50 100 150 200 250 300
SEP (FEM) SDP (FEM)
Fig. 6.10(a) The effect of fiber arrangements on the stress and strain curves of 150
Fig. 6.10(a) The effect of fiber arrangements on the stress and strain curves of 150