(* Modules for implementation of Classical Lamination Theory *)
(* by Zafer Gurdal *)
(* Copyright 1993  *)
(* Department of Engineering Science and Mechanics *)
(* Virginia Polytechnic Institute and State University *)
(* Blacksburg, VA 24061-0219 *)
(* 703-231-5905 *)
(* e-mail: gurdal@gurdal.esm.vt.edu (NeXT Mail Accepted) *)
 
BeginPackage["CompositeLaminate`"]

CLT::usage = "The following modules are available in CLT
		MaterialProperties
		Qmatrix
		Nu21
		InvariantU
		Tmatrix
		Qbarmatrix
		Amatrix
		Bmatrix
		Dmatrix
		AddSubLaminae
		"

Qmatrix::usage = "Qmatrix[{E1, E2, G12, nu12}] 
	Calculates the reduced stiffness matrix of a
	layer with fiber direction modulus E1, transverse
	modulus E2, shear modulus G12, and Poisson's
	ratio nu12 input as a list. If the ply properties 
	are not provided numerically, then it returns the 
	relations for the elements of the reduced stiffness 
	matrix in terms of the layer properties." 

Nu21::usage = "Nu21[{E1, E2, G12, nu12}] 
	Calculates the minor Poisson's ratio in terms of
	the fiber direction modulus E1, transverse
	modulus E2, and Poisson's ratio nu12 input as a list. 
	If the ply properties are not provided numerically, 
	then it returns the relation for the minor Poisson's 
	ration in terms of the other properties." 

InvariantU::usage = "InvariantU[{E1, E2, G12, nu12}] 
	Calculates the material vector U from the properties 
	of a unidirectional layer and returns it in a list
	{U1, U2, U3, U4, U5}. If the ply properties are not 
	provided numerically, then the returned list is in
	terms of the elements of the reduced stiffness matrix 
	Qmatrix." 

Tmatrix::usage = "Tmatrix[theta] 
	Calculates the transformation matrix for a specified
	fiber orientation angle \"theta\" . If the variable
	\"theta\" is a list (ie., theta = {th1, th2, th3, etc.}),
	then the Tmatrix will return a set of matrices 
	corresponding to each of the angles (not operational yet)" 

Qbarmatrix::usage = "Qbarmatrix[theta, uvect] 
	Given the material property vector \"uvect\" (uvect must 
	be a list), Qbarmatrix calculates the stiffness matrix of 
	a layer rotated by an angle \"theta\" from the reference 
	x-y axes.  If the variable \"theta\" is a list 
	(ie., theta = {th1, th2, th3, etc.}), then the Qbarmatrix 
	will return a set of matrices corresponding to each of 
	the angles (not operational yet)" 

Amatrix::usage = "Amatrix[theta, uvect, thick] 
	Calculates the A matrix for a laminate whose 
	stacking sequence is a list of angles \"theta\"
	(for example theta = {45, -45, -45, 45} ).  
	Individual layer thickness is specified in 
	\"thick\" , and the material property vector is in 
	list \"uvect\" ." 
  
Bmatrix::usage = "Bmatrix[theta, uvect, thick] 
	Calculates the B matrix for a laminate whose 
	stacking sequence is a list of angles \"theta\"
	(for example theta = {45, -45, -45, 45}).  
	Individual layer thickness is specified in 
	\"thick\" , and the material property vector is in 
	list \"uvect\" ." 

Dmatrix::usage = "Dmatrix[theta, uvect, thick] 
	Calculates the D matrix for a laminate whose 
	stacking sequence is a list of angles \"theta\"
	(for example theta = {45, -45, -45, 45}).  
	Individual layer thickness is specified in 
	\"thick\" , and the material property vector is in 
	list \"uvect\" ." 

MaterialProperties::usage = "MaterialProperties[\"material\", option] 
	Returns typical elastic properties of the material system
	\"material\"  in a list {E1, E2, G12, nu12}. The default units 
	are the customary U.S. units (psi).  For metric units use the 
	option, MetricUnits -> True, in calling MaterialProperties. 
	Properties for the following list of material systems are
	currently available, 
	\"T300/5208\", \"AS4/3501\", \"Br/Ep\", \"Gl/Ep\", \"Ar/Ep\".
	The first two are commonly used Graphite/Epoxy systems.  The
	last three are typical Boron, Glass, and Aramid fibers in
	Epoxy resin."

MetricUnits::usage = "MetricUnits is an option of MaterialProperties 
	which specifies whether to return the material properties in 
	U.S. customary units or metric units.  The default is 
	MetricUnits -> False." 

AddSubLaminae::usage = "AddSubLaminae[{Ain,Bin,Din,th}] Calculates 
	the A,B,D matrices for a laminate which is formed by joining 
	two or more sublaminates. Input requires the list form of the
	A,B,D matrices of each sublaminate in the same order as the
	stacking sequence (i.e., Ain = {Asub1,Asub2,..}). th is also 
	a list of the sublaminate thicknesses.
	Output will be in the
	form of a list {A,B,D}, where the letters A,B,D correspond to 
	the stiffness matrices for the overall laminate."

Begin["`Private`"]

Options[MaterialProperties] =
	{MetricUnits -> False}

MaterialProperties[material_, opts___Rule] :=
   Module[ { unitm},
	unitm = MetricUnits /. {opts} /. Options[MaterialProperties];
	If[ !unitm , 
(* 			properties in U.S. conventional units 	 *)
		If[ material == "T300/5208" , 
			{26.25 10^6, 1.49 10^6, 1.04 10^6, 0.28}, 
		If[ material == "AS4/3501"  , 
			{20.01 10^6, 1.30 10^6, 1.03 10^6, 0.30},
		If[ material == "Br/Ep"     ,
			{29.59 10^6, 2.68 10^6, 0.811 10^6, 0.23},
		If[ material == "Gl/Ep"     ,
			{5.598 10^6, 1.20 10^6, 0.600 10^6, 0.26},
		If[ material == "Ar/Ep"     ,
			{11.02 10^6, 0.798 10^6,0.334 10^6, 0.34},  
		  	"material must be one of \"T300/5208\", 
			\"AS4/3501\", \"Br/Ep\", \"Gl/Ep\", \"Ar/Ep\" "
		   ] ] ] ] ] , 
(* 			properties in metric units 		 *)
		If[ material == "T300/5208" , 	
			{181.0 10^9, 10.3 10^9, 7.17 10^9, 0.28}, 
		If[ material == "AS4/3501"  ,
			{138.0 10^9, 8.96 10^9, 7.10 10^9, 0.30},
		If[ material == "Br/Ep"     ,
			{204.0 10^9, 18.5 10^9, 5.59 10^9, 0.23},
		If[ material == "Gl/Ep"     ,
			{38.6  10^9, 8.27 10^9, 4.14 10^9, 0.26},
		If[ material == "Ar/Ep"     ,
			{76.0  10^9, 5.50 10^9, 2.30 10^9, 0.34},
		        "material must be one of \"T300/5208\", 
			\"AS4/3501\", \"Br/Ep\", \"Gl/Ep\", \"Ar/Ep\" "
		  ] ] ] ] ]  
	  ]
	]

Qmatrix[prop_List] :=  
   Module[ { div1, qmatr },
	If[And @@ Map[NumberQ, prop], 
		nu21 =(prop[[2]]/prop[[1]])*prop[[4]], 
		If[ nucheck == 1011, nu21 = nu21, nu21= "Nu21"] ];
	div1 = 1 - prop[[4]]*nu21;
	qmatr = { { prop[[1]] / div1, prop[[4]]*prop[[2]] / div1, 0},
		  { prop[[4]]*prop[[2]] / div1, prop[[2]] / div1, 0},
		  { 0, 0, prop[[3]] } }
	]

Nu21[prop_List] := 
    Module[ {  },
	nucheck = 1011;
	nu21 = (prop[[2]] / prop[[1]] )*prop[[4]]
	 ]

InvariantU[prop_List] :=  
   Module[ { div1, q11, q22, q12, q66, empt, u},
	If[ And @@ Map[NumberQ, prop],
            	{{q11,q12,empt},{q12,q22,empt},{empt,empt,q66}} = 
		  Qmatrix[prop],
		{ q11 ="Q11", q22 = "Q22", q12 = "Q12", q66 = "Q66"} ];
	u = { (1/ 8)*( 3*q11 + 3*q22 + 2*q12 + 4*q66), 
	      (1/ 2)*(q11 - q22), 
	      (1/ 8)*( q11 + q22 - 2*q12 - 4*q66), 
	      (1/ 8)*( q11 + q22 + 6*q12 - 4*q66), 
	      (1/ 8)*( q11 + q22 - 2*q12 + 4*q66)  }	
         ]

Tmatrix[theta_] :=  
   Module[ { m, n, rthet, t},
	If[ NumberQ[theta], rthet = N[theta*Degree], rthet = theta ];
	m = Cos[ rthet];
	n = Sin[ rthet];
	t = {   { m^2 , n^2 , 2 m n },
		{ n^2 , m^2 , -2 m n}, 
		{-m n , m n , m^2-n^2} }	
         ]

Qbarmatrix[theta_, u_List] :=  
   Module[{v1, v2, v3, v4, rthet, q11 , q12 , q16, q22 , q26, q66, qb},
	If[ NumberQ[theta], rthet = N[theta*Degree], rthet = theta ];
	v1 = Cos[2*rthet];
	v2 = Sin[2*rthet];
	v3 = Cos[4*rthet];
	v4 = Sin[4*rthet];
	q11 = u[[1]]+u[[2]]*v1+u[[3]]*v3;
	q12 = u[[4]]-u[[3]]*v3;
	q16 =        u[[2]]*v2/2+u[[3]]*v4;
	q22 = u[[1]]-u[[2]]*v1+u[[3]]*v3;
	q26 =        u[[2]]*v2/2-u[[3]]*v4;
	q66 = u[[5]]-u[[3]]*v3;
	qb = {{ q11 , q12 , q16 },
		{ q12 , q22 , q26 }, 
		{ q16 , q26 , q66 }}	
         ]

Amatrix[theta_, u_List, thick_] :=  
   Module[{nn, v1, v2, v3, v4, rthet, a, a11, a12, a16, a22, a26, a66},
	nn=Length[theta];
	If[VectorQ[theta,NumberQ],rthet=N[theta*Degree],rthet=theta];
	v0 = nn*thick;
	v1 = thick*Apply[Plus, Map[Cos,2*rthet]];
	v2 = thick*Apply[Plus, Map[Sin,2*rthet]];
	v3 = thick*Apply[Plus, Map[Cos,4*rthet]];
	v4 = thick*Apply[Plus, Map[Sin,4*rthet]];
	a11 = u[[1]]*v0+u[[2]]*v1+u[[3]]*v3;
	a12 = u[[4]]*v0-u[[3]]*v3;
	a16 = u[[2]]*v2/2+u[[3]]*v4;
	a22 = u[[1]]*v0-u[[2]]*v1+u[[3]]*v3;
	a26 = u[[2]]*v2/2-u[[3]]*v4;
	a66 = u[[5]]v0-u[[3]]*v3;
	a = {{a11, a12, a16},{a12, a22, a26},{a16, a26, a66}}//N         		
	 ]
   
Bmatrix[theta_, u_List, thick_] :=  
   Module[{nn, v1, v2, v3, v4, rthet, b, b11, b12, b16, b22, b26, b66,
		ho2, zpos2},    
	nn=Length[theta];
	If[VectorQ[theta,NumberQ],rthet=N[theta*Degree],rthet=theta];
	ho2 = -nn*thick/2.0;
	zpos2 = Array[(ho2+#*thick)^2 - (ho2+(#-1)*thick)^2 &, nn];
	v1 = 1/2*Apply[Plus,MapThread[Times,{zpos2,Map[Cos,2*rthet]}]];
	v2 = 1/2*Apply[Plus,MapThread[Times,{zpos2,Map[Sin,2*rthet]}]];
	v3 = 1/2*Apply[Plus,MapThread[Times,{zpos2,Map[Cos,4*rthet]}]];
	v4 = 1/2*Apply[Plus,MapThread[Times,{zpos2,Map[Sin,4*rthet]}]];
	b11 =  u[[2]]*v1+u[[3]]*v3;
	b12 = -u[[3]]*v3;
	b16 =  u[[2]]/2*v2+u[[3]]*v4 ;
	b22 = -u[[2]]*v1+u[[3]]*v3;
	b26 =  u[[2]]/2*v2-u[[3]]*v4 ;
	b66 = -u[[3]]*v3;
	b = {{b11, b12, b16},{b12, b22, b26},{b16, b26, b66}}//N        		
	 ]
   
Dmatrix[theta_, u_List, thick_] :=  
   Module[{nn, v1, v2, v3, v4, rthet, d, d11, d12, d16, d22, d26, d66,
		ho2, zpos3},    
	nn=Length[theta];
	If[VectorQ[theta,NumberQ],rthet=N[theta*Degree],rthet=theta];
	ho2 = -nn*thick/2.0;
	zpos3 = Array[(ho2+#*thick)^3 - (ho2+(#-1)*thick)^3 &, nn];
	v0 = (nn*thick)^3 / 12;
	v1 = 1/3*Apply[Plus,MapThread[Times,{zpos3,Map[Cos,2*rthet]}]];
	v2 = 1/3*Apply[Plus,MapThread[Times,{zpos3,Map[Sin,2*rthet]}]];
	v3 = 1/3*Apply[Plus,MapThread[Times,{zpos3,Map[Cos,4*rthet]}]];
	v4 = 1/3*Apply[Plus,MapThread[Times,{zpos3,Map[Sin,4*rthet]}]];
	d11 = u[[1]]*v0+u[[2]]*v1+u[[3]]*v3;
	d12 = u[[4]]*v0-u[[3]]*v3;
	d16 = u[[2]]*v2/2+u[[3]]*v4;
	d22 = u[[1]]*v0-u[[2]]*v1+u[[3]]*v3;
	d26 = u[[2]]*v2/2-u[[3]]*v4;
	d66 = u[[5]]v0-u[[3]]*v3;
	d = {{d11, d12, d16},{d12, d22, d26},{d16, d26, d66}}//N
         ]
   
AddSubLaminae[A_List,B_List,D_List,th_List] :=
		
	Module[{total,abd,nn,zb},
		nn=Length[th];
		total=Apply[Plus,th]; 
		zb = FoldList[ Plus, th[[1]]/2-total/2,
			Array[ (th[[#]] + th[[#+1]])/2 &, nn-1] ];	
		abd={Apply[Plus,A],Apply[Plus,(B+(A*zb))],
			Apply[Plus,D+zb^2*A+2*zb*B]}
	       ]

nucheck = 0   
End[]
EndPackage[]