• Introduction and Definition of General Matrices
• Section 1: Monoclinic Material (m)
• Section 2: Orthotropic Material (o)
• Section 3: Transversely Isotropic Material (t)
• Section 4: Isotropic Material (i)
• Section 5: Constitutive Relations for an Isotropic Material. Lame Coefficients and Engineering Constants
• Section 6: Conclusion
• Section 7: References

Constitutive Equations for Linear Elastic Materials

by Tancredo Westphal Jr., Dr.-Ing.

Institute of Solid Mechanics, Karlsruhe University, Germany

NOTE: This worksheet demonstrates the use of Maple for calculating constitutive coefficients of linear elastic materials (Hooke's law in the three-dimensional theory of elasticity) by means of coordinate transformations (rotation of coordinates). We start from an anisotropic material embedded in a right-handed orthogonal Cartesian coordinate system  ( )  adopting Voigt's notation. In this case the starting object is a symmetric 6 by 6 matrix characterized by 21 independent coefficients.

Key words: linear elastic, elasticity, constitutive equations, Voigt notation

Introduction and Definition of General Matrices

Starting from an anisotropic material (21 independent coefficients) we obtain the following constitutive equations:
1. Monotropic Material, symmetry with respect to the plane 1-2 (13 independent coefficients);
2. Orthotropic Material, additional symmetry with respect to the plane 2-3 (9 independent coefficients);
3. Transversely Isotropic Material, additional isotropy with respect to the planes 1-2 (5 independent coefficients);
4. Isotropic Material, additional isotropy with respect to the plane 2-3 (2 independent coefficients).

At the end we determine the constitutive equations for an isotropic material in terms of the engineering constants (E, ) and Lame's coefficients ( , ).

The transformation for rotations from old coordinates  ( )  to new coordinates  ( )  is given by   , where   . We consider Einstein's notation, where repeated indices imply summing over the range of the index, 3 in this case.

In place of the conventional tensor notation     we define the matrix relationship Sigma33:=C3333.Epsilon33, where the numbers identify the order of the corresponding matrices.

In terms of Voigt's notation  11->1,  22->2,  33->3,  23-> 4,  13->5,  12->6  these matrix relationships are given by Sigma6:=C66.Epsilon6. According to Einstein's notation repeated indices imply summing from 1 to 6 in this case.

For the four special cases above we will append to the respective matrices the labels _m, _o, _t and _i. For the better understanding of the calculation steps see e.g. Ref. [2].

>	restart; with(LinearAlgebra): alias(C=cos(theta)): alias(S=sin(theta)):

>	Sigma33:=Matrix(3,3,shape=symmetric,[[sigma11],[sigma12,sigma22],[sigma13,sigma23,sigma33]]); Sigma6:=Vector([sigma11,sigma22,sigma33,sigma23,sigma13,sigma12]); C66:=Matrix(6,6,shape=symmetric,[ [c11], [c12,c22], [c13,c23,c33], [c14,c24,c34,c44], [c15,c25,c35,c45,c55], [c16,c26,c36,c46,c56,c66] ]);

>	Num_aniso:=21: printf("\nAn anisotropic material is characterized by %a independent elastic coefficients, given by a symmetric matrix of the form C66!",Num_aniso);

An anisotropic material is characterized by 21 independent elastic coefficients, given by a symmetric matrix of the form C66!

Section 1: Monoclinic Material (m)

The material presents a symmetry with respect to the plane 1-2. Instead to consider a mirrow transformation in respect to this plane it is possible to use a corresponding right-handed Cartesian coordinate system determined by a rotation of 180 degrees with respect to the axis 3:

>	alpha:=Pi: A33_m:=Matrix(3,3,[ [cos(alpha),sin(alpha),0], [-sin(alpha),cos(alpha),0], [0,0,1] ]);

According to the transformation rule, a tensor of second order     transforms according to   ,  i.e.,   ,  what corresponds to the matrix transformation   :

>	Sigma33_m:=A33_m.Sigma33.Transpose(A33_m);

Now consider the above matrix in terms of Voigt's notation, where the variable 'Voigt' gives the corresponding relationship between the indices:

>	Sigma6_m:=Vector(6): for i from 1 to 3 do   for j from i to 3 do     Voigt:=(-14*i+10*j+3*i*j)/(3*j-4);     Sigma6_m[Voigt]:=Sigma33_m[i,j];   od:j:='j': od:i:='i': print(Transpose(Sigma6_m));

We obtain the matrix A66_m according to the relation Sigma6_m:=A66_m.Sigma6

>	A66_m:=Matrix(6,6): for k from 1 to 6 do   for i from 1 to 6 do     x:=Sigma6_m[i];     if(has(x,Sigma6[k])) then       A66_m[i,k]:=A66_m[i,k]+select(has,x,Sigma6[k])/Sigma6[k]*sign(x);     fi:   od:i:='i': od:k='k': print(A66_m);

The following transformation is the same as before, used to calculate Sigma33_m, now for matrices 6 by 6:

>	C66_m:=A66_m.C66.Transpose(A66_m);

The constitutive equations of both systems are the same (for this transformation the result can be directly evaluated as follows):

>	C66_m:=(C66+C66_m)/2;

Transforming the above matrix to a set facilitates to identify the number of different coefficients in this matrix:

>	Num_mono:=nops(convert(C66_m,set)): printf("\nA monoclinic material is characterized by %a independent elastic coefficients, given by a symmetric matrix of the form C66_m!",Num_mono);

A monoclinic material is characterized by 13 independent elastic coefficients, given by a symmetric matrix of the form C66_m!

Section 2: Orthotropic Material (o)

The material presents an additional symmetry with respect to the plane 2-3. This corresponds to a rotation of 180 degrees with respect to the axis 1. The other steps are similar to those in the preceding section.

>	A33_o:=Matrix(3,3,[ [1,0,0], [0,cos(alpha),sin(alpha)], [0,-sin(alpha),cos(alpha)] ]);

>	Sigma33_o:=A33_o.Sigma33.Transpose(A33_o);

Sigma6_o:=A66_o.Sigma6

>	Sigma6_o:=Vector(6): for i from 1 to 3 do   for j from i to 3 do     Voigt:=(-14*i+10*j+3*i*j)/(3*j-4);     Sigma6_o[Voigt]:=Sigma33_o[i,j];   od:j:='j': od:i:='i': print(Transpose(Sigma6_o));

>	A66_o:=Matrix(6,6): for k from 1 to 6 do   for i from 1 to 6 do     x:=Sigma6_o[i];     if(has(x,Sigma6[k])) then       A66_o[i,k]:=A66_o[i,k]+select(has,x,Sigma6[k])/Sigma6[k]*sign(x);     fi:   od:i:='i': od:k='k': print(A66_o);

>	C66_o:=A66_o.C66_m.Transpose(A66_o);

The constitutive equations of both systems are the same:

>	C66_o:=(C66_m+C66_o)/2;

>	Num_ortho:=nops(convert(C66_o,set)): printf("\nAn orthotropic material is characterized by %a independent elastic coefficients, given by a symmetric matrix of the form C66_o!",Num_ortho);

An orthotropic material is characterized by 9 independent elastic coefficients, given by a symmetric matrix of the form C66_o!

Section 3: Transversely Isotropic Material (t)

We consider additionally that the plane 1-2 is a plane of isotropy:

>	A33_t:=Matrix(3,3,[ [cos(theta),sin(theta),0], [-sin(theta),cos(theta),0], [0,0,1] ]);

>	Sigma33_t:=A33_t.Sigma33.Transpose(A33_t);

Sigma6_t:=A66_t.Sigma6

>	Sigma6_t:=Vector(6): for i from 1 to 3 do   for j from i to 3 do     Voigt:=(-14*i+10*j+3*i*j)/(3*j-4);     Sigma6_t[Voigt]:=expand(Sigma33_t[i,j]);   od:j:='j': od:i:='i': print(Sigma6_t);

>	A66_t:=Matrix(6,6): for k from 1 to 6 do   for i from 1 to 6 do     x:=Sigma6_t[i];     if(has(x,Sigma6[k])) then       A66_t[i,k]:=expand(factor(A66_t[i,k]+select(has,x,Sigma6[k])/Sigma6[k]));     fi:   od:i:='i': od:k='k': print(A66_t);

>	C66_t:=A66_t.C66_o.Transpose(A66_t): C66_t:=map(simplify,C66_t):

The equality of the systems _o and _t leads to 'Sys_t' non-trivial equations:

>	Sys_t:=0: for i from 1 to 6 do   for j from 1 to 6 do     if(C66_t[i,j]-C66_o[i,j]<>0) then       Sys_t:=Sys_t+1:     fi:   od:j:='j': od:i:='i': Sys_t;

Setting equality of the matrices for both systems leads to a system of linear equations:

>	EqSys_t:=Vector(Sys_t): k:=1: for i from 1 to 6 do   for j from 1 to 6 do     if(C66_t[i,j]-C66_o[i,j]<>0) then #      print(i,j,C66_t[i,j]-C66_o[i,j]);       EqSys_t[k]:=C66_t[i,j]-C66_o[i,j];       k:=k+1:     fi:   od:j:='j': od:i:='i':

Now we identify the coefficients that are involved in the equation system. As can be seen, from the 9 remaining independent coefficients c33 is the only coefficient not involved in these relations:

>

total_t:=0: for i from 1 to 6 do   for j from 1 to 6 do     yes:=0;     for a from 1 to Sys_t while yes = 0 do             if(evalb(select(has,expand(EqSys_t[a]),c||i||j)<>0))then         yes:=1;         total_t:=total_t+1;         coef||total_t:=c||i||j;       fi;     od:a:='a':   od:j:='j': od:i:='i': C_t:=Vector(total_t): for i from 1 to total_t do   C_t[i]:=coef||i; od:i:='i': Transpose(C_t);

We obtain now the corresponding coefficient matrix K_t

>	EqSys_t:=map(expand,EqSys_t): K_t:=Matrix(Sys_t,total_t): for k from 1 to total_t do   for i from 1 to Sys_t do     x:=EqSys_t[i];     if(has(x,C_t[k])) then       K_t[i,k]:=expand(select(has,x,C_t[k])/C_t[k]);     fi:   od:i:='i': od:k='k': evalm(K_t);

The rank of the coefficient matrix gives the number of independent equations:

>	Rank_t:=Rank(K_t);

The corresponding vectors giving these independent equations are calculated now:

>	CS_t:=map(simplify,ColumnSpace(Transpose(K_t)));

And here are these equations:

>	for i from 1 to Rank_t do   g||i:=CS_t[i].C_t;   print(g||i); od:i:='i':

Solving the above 4 equations for the 4 coefficients c22,c23,c55,c66 gives these as function of the remaing ones in C_t:

>	sol_t:=solve({seq(g||i,i=1..4)},{c22,c23,c55,c66}); assign(sol_t);

After assigning the solution of the system all the equations must be satisfied (equal to zero):

>	Ok_t:=evalm(K_t.C_t): map(simplify,Ok_t);

It is only necessary to perform a simplification:

>	for i from 1 to 6 do   for j from 1 to 6 do     C66_t[i,j]:=factor(C66_t[i,j]);   od:j:='j': od:i:='i': print(C66_t);

>	Num_trans:=Num_ortho-Rank_t: printf("\nA transversely isotropic material is characterized by %a independent elastic coefficients, given by a symmetric matrix of the form C66_t!",Num_trans);

A transversely isotropic material is characterized by 5 independent elastic coefficients, given by a symmetric matrix of the form C66_t!

Section 4: Isotropic Material (i)

The plane 2-3 is an additional plane of isotropy. The necessary calculations are similar to those in the preceding section:

>	A33_i:=Matrix(3,3,[ [1,0,0], [0,cos(theta),sin(theta)], [0,-sin(theta),cos(theta)] ]);

>	Sigma33_i:=A33_i.Sigma33.Transpose(A33_i);

Sigma6_i=A66_i.Sigma6

>	Sigma6_i:=Vector(6): for i from 1 to 3 do   for j from i to 3 do     Voigt:=(-14*i+10*j+3*i*j)/(3*j-4);     Sigma6_i[Voigt]:=expand(Sigma33_i[i,j]);   od:j:='j': od:i:='i': print(Sigma6_i);

>	A66_i:=Matrix(6,6): for k from 1 to 6 do   for i from 1 to 6 do     x:=Sigma6_i[i];     if(has(x,Sigma6[k])) then       A66_i[i,k]:=expand(factor(A66_i[i,k]+select(has,x,Sigma6[k])/Sigma6[k]));     fi:   od:i:='i': od:k='k': print(A66_i);

>	C66_i:=A66_i.C66_t.Transpose(A66_i): C66_i:=map(simplify,C66_i):

The equality of the systems _t and _i  leads to 'Sys_i' non-trivial equations:

>	Sys_i:=0: for i from 1 to 6 do   for j from 1 to 6 do     if(simplify(C66_i[i,j]-C66_t[i,j])<>0) then       Sys_i:=Sys_i+1:     fi:   od:j:='j': od:i:='i': Sys_i;

Seting equality of the matrices for both systems:

>	EqSys_i:=Vector(Sys_i): k:=1: for i from 1 to 6 do   for j from 1 to 6 do     if(simplify(C66_i[i,j]-C66_t[i,j])<>0) then #      print(i,j,simplify(C66_i[i,j]-C66_t[i,j]));       EqSys_i[k]:=simplify(C66_i[i,j]-C66_t[i,j]);       k:=k+1:     fi:   od:j:='j': od:i:='i':

Identification of the involved coefficients:

>

total_i:=0: for i from 1 to 6 do   for j from 1 to 6 do     yes:=0;     for a from 1 to Sys_i while yes=0 do       if(evalb(select(has,expand(EqSys_i[a]),c||i||j)<>0))then         yes:=1;         total_i:=total_i+1;         coef||total_i:=c||i||j;       fi:     od:a:='a':   od:j:='j': od:i:='i': total_i;

>	C_i:=Vector(total_i): for i from 1 to total_i do   C_i[i]:=coef||i od:i:='i': Transpose(C_i);

>	C_i:=convert(convert(convert(C_i,set),list),Vector): Transpose(C_i);

>	EqSys_i:=map(expand,EqSys_i): K_i:=Matrix(Sys_i,5):

>	for k from 1 to 5 do   for i from 1 to Sys_i do     x:=EqSys_i[i];     if(has(x,C_i[k])) then       K_i[i,k]:=expand(select(has,x,C_i[k])/C_i[k]);     fi:   od:i:='i': od:k='k': evalm(K_i);

>	Rank_i:=Rank(K_i);

>	CS_i:=map(simplify,ColumnSpace(Transpose(K_i)));

>	for i from 1 to Rank_i do   g||i:=CS_i[i].C_i;   print(g||i); od:i:='i':

Solving the above 3 equations for the 3 coefficients c13,c33,c44 gives these as function of the remaing ones in C_i:

>	sol:=solve({seq(g||i,i=1..Rank_i)},{c13,c33,c44}); assign(sol);

Are all the equations satisfied?

>	Ok_i:=evalm(K_i.C_i): map(simplify,Ok_i);

>	for i from 1 to 6 do   for j from 1 to 6 do     C66_i[i,j]:=factor(C66_i[i,j]);   od:j:='j': od:i:='i': print(C66_i);

>	Num_iso:=Num_trans-Rank_i: printf("\nAn isotropic material is characterized by %a independent elastic coefficients, given by a symmetric matrix of the form C66_i!",Num_iso);

An isotropic material is characterized by 2 independent elastic coefficients, given by a symmetric matrix of the form C66_i!

Section 5: Constitutive Relations for an Isotropic Material. Lame Coefficients and Engineering Constants

S: compliance matrix

>	S_i:=1/C66_i;

Engi : Engineering Constants
Lame : Lame Coefficients

>	Engi:=Vector(2,[E,nu]); Lame:=Vector(2,[lambda,mu]);

The engineering constants have physical meaning:     and     (T: Transversal, L: Longitudinal):

>	f1:=S_i[1,1]=1/Engi[1]; f2:=S_i[1,2]=-Engi[2]/Engi[1];

Constitutive equations for isotropic materials in terms of engineering constants:

>	solEngi:=solve({f1,f2},{c11,c12}): SEngi_i:=map(factor,subs(solEngi,S_i)); CEngi_i:=map(factor,1/SEngi_i);

Lame's coefficients in terms of engineering constants:

>	L1:=mu=E/(2*(1+nu)); L2:=lambda=nu*E/((1+nu)*(1-2*nu));

Constitutive equations for isotropic materials in terms of Lame's coefficients. Observe that CLame_i is a very neat matrix:

>	solLame:=solve({L1,L2},{E,nu}); SLame_i:=map(factor,subs(solLame,SEngi_i)); CLame_i:=1/SLame_i;

>

Section 6: Conclusion

We showed how Maple can easily handle the tensor transformations that are necessary in order to treat some important material symmetries of mechanical materials, giving the corresponding constitutive equations in a very elegant way. This worksheet can be of great help for students learning the fundamentals of the theory of elasticity, to whom it was originally developed.

Section 7: References

There exists a lot of books on elasticity theory and continuum mechanics treating constitutive equations for linear elastic materials. In writing this program I used specially:
[1] Boresi, A.P. and Chong, K.P.: Elasticity in engineering mechanics, 2nd Ed. New York : John Wiley, 2000.
[2] Chung, T.J.: Continuum Mechanics. Englewood Cliffs, NJ : Prentice-Hall, 1988.
[3] Jones, R.M.: Mechanics of Composite Materials. New York : McGraw-Hill, 1975.
[4] Soutas-Little, R.W.: Elasticity. Mineola, N.Y. : Dover, 1999.

Disclaimer:  While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.