Predicting Young's Modulus of Banana-Epoxy Composite Using MATLAB

[barbarahowser]

Hi,
I am trying to analyze a composite part which is made of carbon/epoxy and has [+45/-45]2s as the layup. I got the following material properties for carbon/epoxy from the textbook "Engineering Mechanics of Composite Materials", by I.M. Daniel & O. Ishai,
Longitudinal Modulus, E1=147Gpa
Transverse In-plane, E2=10.3Gpa
Transverse out of plane, E3=10.3Gpa
In plane shear modulus, G12=7Gpa
Out of plane shear modulus, G23=3.7Gpa
Out of plane shear modulus, G13=7Gpa
Major in-plane Poisson's ratio, niu12=0.27
Major in-plane Poisson's ratio, niu23=0.54
Major in-plane Poisson's ratio, niu13=0.27

The coordinate system of the part is shown in figure attached with the fiber laid up in the YZ plane. Can someone please tell me how I can transform the material properties given above into the coordinate system shown in (composite part.jpeg)

Thanks,
Barbara

 

[PerennialII]

Hi Barbara,

you've a reference of lamina & laminate analysis nearby (that textbook might do it although haven't got it myself so can't check)? 1st you derive the properties for different orthotropic lamina (you've the basic elements there which you've collected), like the stiffness or compliance matrix (depending on what you're after), and then you "stack them" in the laminate part (laminate=several lamina) of the analysis (which can then contain for example a stress analysis of some "element" using the laminate in question). The properties are transformed between coordinate systems in the process, typically using so called transformation matrices.

It's somewhat lengthy to present, so if you'd have a book on the subject available with a consistent run through of it all it would help a lot.

 

[barbarahowser]

Hi,
Thanks for the reply. Here is what I did,
1. Calculate the compliance matrix using the available material properties.
2. Calculate transformation matrices for +45 and -45. This will be transformation about X-axis (i.e rotation about X-axis) since X-axis is the perpendicular out-of-plane axis.
3. Calculate ABD matrix using lamination theory
4. Calculate inverse of ABD matrix
5. Calculate material properties using the inverse of ABD matrix

The problem is that using this code, I get the out-of plane stiffness EX higher than the in-plane stiffness EY and EZ which I think is wrong since out-of-plane stiffness should always be lower.

Can you tell me whether the procedure I am using is right or not? I would really appreciate any help on this!

Below is the MATLAB code used according to procedure given above,

format long
%Material Properties
E1=147e9
E2=10.3e9
E3=10.3e9
G12=7e9
G23=3.7e9
G13=7e9
niu12=0.27
niu23=0.54
niu13=0.27
tply=6.35e-4
t=0.00508

S=[(1/E1) (-niu12/E1) (-niu13/E1) 0 0 0;(-niu12/E1) 1/E2 (-niu23/E2) 0 0 0;(-niu13/E1) (-niu23/E2) 1/E3 0 0 0;0 0 0 1/G23 0 0;0 0 0 0 1/G13 0;0 0 0 0 0 1/G12]
Q = inv(S)

%Calculation of Qbar matrix for different fibre orientations
m=cos(pi/4)
n=sin(pi/4)
%Roattion about X axis
Tsigma=[1 0 0 0 0 0;0 m^2 n^2 2*m*n 0 0;0 n^2 m^2 -2*m*n 0 0;0 -m*n m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Tepsilon=[1 0 0 0 0 0;0 m^2 n^2 m*n 0 0;0 n^2 m^2 -m*n 0 0;0 -2*m*n 2*m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Qbar45=inv(Tsigma)*Q*Tepsilon

m=cos(-pi/4)
n=sin(-pi/4)
Rotation about X axis
Tsigma=[1 0 0 0 0 0;0 m^2 n^2 2*m*n 0 0;0 n^2 m^2 -2*m*n 0 0;0 -m*n m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Tepsilon=[1 0 0 0 0 0;0 m^2 n^2 m*n 0 0;0 n^2 m^2 -m*n 0 0;0 -2*m*n 2*m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]

Qbar_45=inv(Tsigma)*Q*Tepsilon

h0=-2.54e-3
h1=-1.905e-3
h2=-1.27e-3
h3=-6.35e-4
h4=0
h5=6.35e-4
h6=1.27e-3
h7=1.905e-3
h8=2.54e-3

%ABD matirx
A=Qbar45*(h1-h0)+Qbar_45*(h2-h1)+Qbar45*(h3-h2)+Qbar_45*(h4-h3)+Qbar_45*(h5-h4)+Qbar45*(h6-h5)+Qbar_45*(h7-h6)+Qbar45*(h8-h7)
B=(1/2)*(Qbar45*(h1^2-h0^2)+Qbar_45*(h2^2-h1^2)+Qbar45*(h3^2-h2^2)+Qbar_45*(h4^2-h3^2)+Qbar_45*(h5^2-h4^2)+Qbar45*(h6^2-h5^2)+Qbar_45*(h7^2-h6^2)+Qbar45*(h8^2-h7^2))
D=(1/3)*(Qbar45*(h1^3-h0^3)+Qbar_45*(h2^3-h1^3)+Qbar45*(h3^3-h2^3)+Qbar_45*(h4^3-h3^3)+Qbar_45*(h5^3-h4^3)+Qbar45*(h6^3-h5^3)+Qbar_45*(h7^3-h6^3)+Qbar45*(h8^3-h7^3))

ABD=[A B;B D]
abd=inv(ABD)

%Calculate material properties
EZ=1/((abd(1,1))*t)
Ey=1/((abd(2,2))*t)
EX=1/((abd(3,3))*t)

niuzy=-(abd(2,1)/abd(1,1))
niuxy=-(abd(3,2)/abd(2,2))
niuxz=-(abd(3,1)/abd(1,1))

Gzy=1/((abd(6,6))*t)
Gxz=1/((abd(5,5))*t)
Gxy=1/((abd(4,4))*t)



Thanks,
Barbara

 

[barbarahowser]

Im sorry there was a probelm in the last part of the code. Here is the correct one,

format long
%Material Properties
E1=147e9
E2=10.3e9
E3=10.3e9
G12=7e9
G23=3.7e9
G13=7e9
niu12=0.27
niu23=0.54
niu13=0.27
tply=6.35e-4
t=0.00508

S=[(1/E1) (-niu12/E1) (-niu13/E1) 0 0 0;(-niu12/E1) 1/E2 (-niu23/E2) 0 0 0;(-niu13/E1) (-niu23/E2) 1/E3 0 0 0;0 0 0 1/G23 0 0;0 0 0 0 1/G13 0;0 0 0 0 0 1/G12]
Q = inv(S)

%Calculation of Qbar matrix for different fibre orientations
m=cos(pi/4)
n=sin(pi/4)
%Roattion about X axis
Tsigma=[1 0 0 0 0 0;0 m^2 n^2 2*m*n 0 0;0 n^2 m^2 -2*m*n 0 0;0 -m*n m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Tepsilon=[1 0 0 0 0 0;0 m^2 n^2 m*n 0 0;0 n^2 m^2 -m*n 0 0;0 -2*m*n 2*m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Qbar45=inv(Tsigma)*Q*Tepsilon

m=cos(-pi/4)
n=sin(-pi/4)
%Rotation about X axis
Tsigma=[1 0 0 0 0 0;0 m^2 n^2 2*m*n 0 0;0 n^2 m^2 -2*m*n 0 0;0 -m*n m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]
Tepsilon=[1 0 0 0 0 0;0 m^2 n^2 m*n 0 0;0 n^2 m^2 -m*n 0 0;0 -2*m*n 2*m*n (m^2-n^2) 0 0;0 0 0 0 m -n;0 0 0 0 n m]

Qbar_45=inv(Tsigma)*Q*Tepsilon

h0=-2.54e-3
h1=-1.905e-3
h2=-1.27e-3
h3=-6.35e-4
h4=0
h5=6.35e-4
h6=1.27e-3
h7=1.905e-3
h8=2.54e-3

%ABD matirx
A=Qbar45*(h1-h0)+Qbar_45*(h2-h1)+Qbar45*(h3-h2)+Qbar_45*(h4-h3)+Qbar_45*(h5-h4)+Qbar45*(h6-h5)+Qbar_45*(h7-h6)+Qbar45*(h8-h7)
B=(1/2)*(Qbar45*(h1^2-h0^2)+Qbar_45*(h2^2-h1^2)+Qbar45*(h3^2-h2^2)+Qbar_45*(h4^2-h3^2)+Qbar_45*(h5^2-h4^2)+Qbar45*(h6^2-h5^2)+Qbar_45*(h7^2-h6^2)+Qbar45*(h8^2-h7^2))
D=(1/3)*(Qbar45*(h1^3-h0^3)+Qbar_45*(h2^3-h1^3)+Qbar45*(h3^3-h2^3)+Qbar_45*(h4^3-h3^3)+Qbar_45*(h5^3-h4^3)+Qbar45*(h6^3-h5^3)+Qbar_45*(h7^3-h6^3)+Qbar45*(h8^3-h7^3))

ABD=[A B;B D]
abd=inv(ABD)

%Calculate material properties
EX=1/((abd(1,1))*t)
Ey=1/((abd(2,2))*t)
EZ=1/((abd(3,3))*t)

niuxy=-(abd(2,1)/abd(1,1))
niuyz=-(abd(3,2)/abd(2,2))
niuxz=-(abd(3,1)/abd(1,1))

Gxy=1/((abd(6,6))*t)
Gyz=1/((abd(5,5))*t)
Gxz=1/((abd(4,4))*t)

And here are the results,
EX = 1.470000000000000e+011 Pa

Ey = 1.112375948628138e+010 Pa

EZ = 1.112375948628138e+010 Pa

niuxy = 0.27000000000000

niuyz = 0.50321074138938

niuxz = 0.27000000000000

Gxy = 7.000000000000000e+009 Pa

Gyz = 7.000000000000000e+009 Pa

Gxz = 3.344155844155844e+009 Pa


Thanks,
Barbara

 

[venkatcad]

Hi anyone help me to predict youngs modulus of banana-epoxy composite using matlab