## Linear Elastic Materials: Matrix of Material Properties of Linear Elastic Materials

#### Learning Outcomes

that in general, after considering that the stress matrix, the strain matrix, and the coefficients matrix are symmetric, there are 21 constants to describe the relationship between the stress matrix and strain matrix.__Identify__why the coefficients matrix needs to be symmetric.__Explain__that for orthotropic materials, the coefficients matrix has 9 constants and the relationship is particular to a specific coordinate system.__Identify__the meaning of the 9 physical constants for orthotropic materials.__Explain__that for transversely isotropic materials, the coefficients matrix has 5 constants and that the relationship is particular to a specific coordinate system.__Identify__**Explain**the meaning of the 5 physical constants for transversely isotropic materials.that for isotropic materials, the coefficients matrix has 2 constants and that the relationship is valid in any coordinate system. Explain the meaning of the 2 physical constants for isotropic materials.__Identify__- Given 6 stress and/or strain components,
the remaining components.__calculate__ the physical meaning of and__Describe__between Young’s modulus, Poisson’s ratio, the shear modulus, and the bulk modulus for isotropic materials. Given any two of the constants above,__differentiate__the remaining ones.__compute__

### Introduction

An elastic material is defined as a material whose Cauchy stress tensor is a function of the current state of deformation independent of the history or path, I.e., . For a non-homogeneous material, this function would depend on the position as well but this dependence is dropped at the moment for brevity. The principle of material frame-indifference restricts the possible forms for to:

Equivalently, the constitutive law for the second Piola-Kirchhoff stress is restricted to have the form:

The equivalence can be shown as follows, given , therefore:

and

Therefore:

The opposite statement can be shown similarly. In other words, for an elastic material, the second Piola-Kirchhoff stress tensor can be viewed as a function of the right stretch tensor independent of the rotation while the Cauchy stress matrix, as it is a spacial measure of stress depends on both and such that . For example, if a material is rotating in space, then, the second Piola-Kirchhoff stress tensor is expected to be constant, while the Cauchy stress matrix will change as it is tied to the spacial coordinates.

#### Small Deformations

Small deformations can be characterized by a deformed configuration function given by:

with being a small parameter and is a vector field independent of the parameter such that:

Following the notation adopted in Taylor’s Theorem, the expression indicates an expression (scalar or tensor) that is directly proportional to .

##### Measures of Deformation

In the case of small deformations, the deformation gradient can be written as:

As a reminder, the displacement gradient can be decomposed into the small strain tensor and the infinitesimal rotation tensor resulting in:

The right Cauchy-Green deformation tensor adopts the following form:

The right stretch tensor can be approximated as:

With the small strain tensor appearing naturally in the approximation. The rotation can be assumed to have the form:

such that:

It is important to note that the small deformations as shown by the above approximations imply both small strains and small rotations. In this case, the small strain and the Green Strain matrices can be assumed equal as they differ by a small term:

The volumetric strain can be evaluated as follows:

which is based on the following approximation (utilizing the relationship between the determinant and the trace as shown in the invariants section):

##### Stress Tensors

Under the assumption of small deformations, the stress tensors are all equal by ignoring the terms :

Similarly

##### Constitutive Relationship

As shown above, for elastic materials, the second Piola-Kirchhoff stress tensor is a function of . Utilizing the assumption of small deformation and using a Taylor expansion around the identity matrix we get:

Assuming that the reference configuration is stress free, i.e., and that the Cauchy and Piola-Kirchhoff stress tensors are equal:

where is a matrix of 81 constants.

### Matrix of Material Properties of Linear Elastic Materials

A linear elastic material is a material that exhibits a linear relationship between the components of the stress tensor and the components of the strain tensor. A linear elastic material constitutive law, under the assumption of small deformation, is fully represented by a linear map between the stress matrix and the infinitesimal strain matrix such that . In general, the stress matrix has components and the strain matrix has components. A linear relationship between these components would be represented by a fourth-order tensor with components. These components are unique for a particular linear elastic material, and thus, they are considered to be material constants. By adopting the Einstein summation convention the linear map between the stress and the strain tensors can be written in component form as

The components of the fourth-order tensor can be visualized in the following form:

For example, the elements of the matrix on the top left of the above equation can be viewed as the matrix of coefficients that are used to calculate the stress component as follows:

#### Reductions Due to Symmetries of the Stress and the Strain

Because the stress matrix is symmetric, i.e., , then this leads to the following relationship among the components of

This can be further explained using an example. Assume a strain state defined by a strain matrix with zero values except for the strain . Then, and . Since , then, . This applies to all the coefficients in the form of and .

Similarly, because the strain matrix is symmetric, i.e., , then this leads to the following relationship among the components of

This can be further explained using an example. Assume a strain state defined by a strain matrix with zero values except for the strain components . Then, . Since we have: , i.e., the coefficients and can be chosen to be equal. This applies to all the coefficients in the form of and .

Therefore, because of the symmetries of the stress and the strain matrices, the number of independent coefficients in the matrix can be reduced to independent coefficients. This simplifies the relationship which can now be represented in matrix form by assuming that the stress and the strain matrices are vectors in with the following forms:

Notice that traditionally, the engineering shear strains , , and are used in the constitutive relationships and these are equal to double the corresponding shear strains. The linear map between the stress and the strain can be represented by a matrix . The relationship between the vector representation of the stress and the vector representation of the strain has the following form:

The above relationship is described by 36 independent coefficients. So far, there are no restrictions imposed on for it to be symmetric. However, in the next section we will show that the restriction imposed by the preservation of energy leads to the symmetry of .

#### Reductions Due to Preservation of Energy

For a linear elastic material, the energy stored during loading and/or the energy released during unloading have to be independent of the path of loading and/or unloading. Thus, for such a material, a potential energy function has to exist such that the stress is derived from this potential energy function as follows:

See chapter 6 in the book “Elasticity: Theory, Applications, and Numerics” by Martin H. Sadd for a detailed description of this concept.

The stress component can be obtained using the following relationship:

The equality of the second partial derivatives of leads to the symmetry of as follows:

In matrix form, the relationship between the vector representation of the stress and the vector representation of the strain has the following form with 21 independent coefficients:

The symmetry of the coefficient matrix is crucial so that the material model would not predict energy creation or dissipation during loading cycles. This can be illustrated with an example in which the relationship between the in-plane strain and in-plane stress of a material is given by the following non-symmetric relationship:

Assume that this material is loaded such that is increased gradually from 0 to 100MPa. Then, is kept constant and is increased gradually from 0 to 100MPa. Then, assume that both stresses are decreased simultaneously from 100MPa to 0. Since this is an elastic material, no energy should be created or dissipated through the entire loading process since the final state of stress is identical to the initial state of stress even though the paths of loading and unloading are different. The energy in a path from state a to state b can be calculated using the equation:

To simplify the integrals, we will use the following relationships:

Therefore, the integrals have the following form:

The loading can be divided into three paths. In the first path, increases from 0 to 100MPa while stays constant at zero. Therefore, the energy can be calculated as follows:

In the second path, stays constant at 100MPa while increases from 0 to 100MPa. Therefore, the energy in the second path can be calculated as follows:

In the third path, , and both change from 100MPa to 0. Therefore, the energy in the third path can be calculated as follows:

Therefore, the total energy after the three paths is equal to:

which is contrary to the expectation that the material model preserves energy. The above loading cycle shows a dissipation of energy. The amount of energy dissipated is related to the difference between the off-diagonal components in the matrix of material coefficients. Had the matrix of coefficients been symmetric, the total energy stored or generated during the loading cycle would have been zero. Can you identify a different loading cycle in which the material would generate energy?

#### Reductions Due to Material Symmetry

The number of constants can be further reduced, depending on the symmetries of the material. If a material possesses a certain symmetry, this implies that the components of do not change upon coordinate transformations within the specified symmetry. Since is a fourth-order tensor, upon changing the coordinate system using an orthogonal matrix , the new material constants have the following form:

If a material’s behaviour is independent of a certain change of the coordinate system described by a particular orthogonal matrix , then the components of should be equal to those of . Such orthogonal matrix is termed a material symmetry. For example, if a material response is symmetric around a plane, then a reflection matrix across this plane would produce the same components for the coefficients matrix. If a material response is the same within any rotation inside a plane, then a rotation within that plane would produce the same components for as those for . The number of the independent constants of can be further reduced from 21, depending on the level of anisotropy of the material. In the following section, the different levels of anisotropy are described, and the physical interpretation of some of the material constants is given.

#### Orthotropic Linear Elastic Materials

In an orthotropic linear elastic material, the 21 material constants are reduced to 9 material constants fully describing the relationship between the stress and the small strain. A proof of this argument is illustrated using the mathematica code found here. Assuming that the basis vectors , , and are the axis of orthotropy, then the relationship between the stress and the strain is given by:

(1)

Notice that symmetry around each plane leads to the normal stresses being independent of the shear strains and vice versa. In other words, if an orthotropic material is pulled along a direction perpendicular to a plane of symmetry, then the associated shear strains are zero. For orthotropic materials, the relationship between the small strains and the stresses is traditionally given as function of nine physical material constants , , , , , , , , and where: : The ratio of the uniaxial stress to uniaxial strain in the direction of the basis vector . : Negative the ratio between the axial strain in the direction of the basis vector to the axial strain in the direction of the basis vector when the material is uniaxially stressed in the direction of the basis vector . : The ratio between the shear stress to the engineering shear strain in the plane defined by the directions of and .

The relationship can be written in terms of the above physical material constants as follows:

The relationship has the following matrix form:

(2)

Notice that symmetry of the matrix implies that

#### Transversely Isotropic Linear Elastic Materials

In a transversely isotropic linear elastic material, five material constants describe the relationship between the stress and the small strain. A proof of this argument is illustrated using the mathematica code found here. Assuming full isotropy in the plane formed by the basis vectors and . i.e., the material response is identical given any rotation around . Also, assuming that the plane formed by the vectors and is a plane of symmetry, then, the relationship between the stress and the strain is given by five material constants:

(3)

For transversely isotropic materials, the relationship between the small strains and the stresses are traditionally given as function of five physical material constants , , , , and where:

: The ratio of the uniaxial stress to uniaxial strain in any direction in the isotropy plane.

: Negative the ratio between the transverse strain to the axial strain when the material is uniaxially stressed in any direction in the isotropy plane.

: The ratio of the uniaxial stress to uniaxial strain in the direction perpendicular to the isotropy plane.

: Negative the ratio between the axial strain to the axial strain when the material is stressed uniaxially in the direction of

: The ratio between the shear stress to the engineering shear strain in the plane defined by the directions of and .

In this case, the relationship between the strains and the stresses can be written as:

The relationship has the following matrix form:

(4)

The symmetry of the matrix implies that

#### Isotropic Linear Elastic Materials

In an isotropic linear elastic material, the response of the material to stress is independent of the coordinate system chosen. In other words, the relationship between the stress and the strain is independent of the coordinate system. This assumption reduces the number of material constants that describe the behaviour of the material in any direction to only two. A proof of this argument is illustrated using the mathematica code found here. The relationship between the stress and the strain components has the following form independent of the chosen coordinate system:

(5)

The two constants that are generally used are the Young’s modulus and Poisson’s ratio , with the following physical definitions: Young’s modulus : The slope of the stress-strain curve during uniaxial tensile tests of a sample of a linear elastic isotropic material. Poisson’s ratio : Negative the ratio of the transverse strain (perpendicular to the applied load) to the axial strain (in the direction of the applied load) during uniaxial tensile tests of a sample of a linear elastic isotropic material.

For such materials, the strains are given by the following equations

The relationship has the following matrix form:

(6)

Where:

The inverse of the above relationship has the following form:

(7)

Note that equations 1, 2, 3, and 4 are only valid in particular coordinate systems while equations 5, 6, and 7 can be used in any coordinate system.

##### Alternative Material Constants

The stress-strain relationships (Equations 6 and 7) rely on the definitions of and . However, other material constants can be used. The wikipedia page “Elastic Modulus” offers an extensive list of such constants. In particular, the shear modulus , the bulk modulus , and the Lamé’s constants and are sometimes used. The following are their respective definitions:

###### Shear Modulus

The shear modulus is the slope of the shear stress versus the engineering shear strain of a linear elastic isotropic material. It is related to and as follows:

The above relationship can be deduced for an isotropic material as follows. We first use Equation 6 but assume that , , and are independent of each other. Assume a state of pure shear stress with being the only non-zero stress component. Using the isotropic linear elastic constitutive law equation (Equation 6), then the only non-zero strain components are . Therefore, the stress and strain matrices are:

We now apply a coordinate transformation in the plane described by and using the rotation matrix :

In the new coordinate system, and have the following forms:

Since the material is isotropic, the relationship shown in Equation 6 applies in the new coordinate system as well. Therefore:

However, we already know from applying the constitutive equation in the original coordinate system that:

However, since we get the required relationship: .

###### Bulk Modulus

The bulk modulus is the slope of the applied hydrostatic stress versus the volumetric strain of a linear elastic isotropic material. It is related to and as follows:

The relationship for the bulk modulus can be obtained as follows. First we assume a hydrostatic state of stress where is the identity matrix, and is a real number representing the applied pressure. In this case, when adopting Equation 6, the developed strains are . Therefore,

Therefore, the slope of the relationship between and the volumetric strain is equal to .

###### LAMÉ’S CONSTANTS:

These constants are defined as follows:

These two constants provide a simplified way of writing 7 as follows.

(8)

or in simple component form

### Videos

### Derivation of the Reduced Number of Independent Matrix Coefficients due to Material Symmetries

The following Mathematica code is used to reduce the number of independent coefficients in the fourth-order tensor. In the first part, the symmetry of the stress tensor and the symmetry of the strain tensor are utilized to reduce the number of independent coefficients from 81 to 36. Then, the preservation of energy concept is utilized to reduce the number of independent coefficients from 36 to 21. Next, reflections along the directions of , , and are applied to find the number of coefficients required for an orthotropic material, with the planes of symmetry being the planes perpendicular to the three basis vectors. In this case, the number of independent constants is shown to be 9. Then, the number of independent constants in a transversely isotropic material is found by utilizing the orthotropic material fourth-order tensor and assuming that the behaviour of the material is independent of a rotation in the plane of and . The resulting number of independent material constants is shown to be 5. Finally, another rotation around is used to show that the number of independent material constants for an isotropic material is just 2. Copy and paste the code into your Mathematica browser, and follow it step by step. You should also show that the results for the transversely isotropic and isotropic materials are independent of the angle of rotation used in the derivation, as long as the angle is not a multiple of .

View Mathematica Code

Ctensor=Table[Row[{c,i,j,k,l}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“A general C matrix relating 9 stress tensor components to 9 strain tensor components:”]

Ctensor//MatrixForm

(*Applying stress symmetry*)

RuleStressSymmetry={};

Do[AppendTo[RuleStressSymmetry,Ctensor[[i,j,k,l]]->Ctensor[[j,i,k,l]]],{i,2,3},{j,1,i-1},{k,1,3},{l,1,3}];

Print[“Stress tensor symmetry leads to the following relations:”,RuleStressSymmetry]

CtensorSt=Ctensor/.RuleStressSymmetry;

Print[“The resulting C matrix from stress tensor symmetry is:”]

CtensorSt//MatrixForm

(*Applying Strain symmetry*)

RuleStrainSymmetry={};

Do[AppendTo[RuleStrainSymmetry,CtensorSt[[i,j,k,l]]->CtensorSt[[i,j,l,k]]],{i,1,3},{j,1,3},{k,2,3},{l,1,k-1}];

Print[“Strain tensor symmetry leads to the following relations:”,RuleStrainSymmetry]

CtensorStStr=CtensorSt/.RuleStrainSymmetry;

Print[“The resulting C matrix from both stress and strain tensors symmetry is:”]

CtensorStStr//MatrixForm

(*Using symmetry due to Energy Preservation*)

RuleEnergySymmetry={};

Do[AppendTo[RuleEnergySymmetry,CtensorStStr[[i,j,k,l]]->CtensorStStr[[k,l,i,j]]],{i,1,3},{j,i,3},{k,1,i},{l,k,If[k

*Print[“Symmetry due to the existence of an energy function leads to the following relations:”,RuleEnergySymmetry]*

CtensorStStrEn=CtensorStStr/.RuleEnergySymmetry;

Print[“The resulting C matrix after using stress, strain and energy symmetries is:”]

CtensorStStrEn//MatrixForm

Print[“The original C tensor and the fully symmetric Ctensor:”]

{Ctensor//MatrixForm,CtensorStStrEn//MatrixForm}

Print[“The difference (minus) between C tensor and the fully symmetric Ctensor:”]

Ctensor-CtensorStStrEn//MatrixForm

(*Adopting the new Ctensor after symmetries*)

Ctensor=CtensorStStrEn;

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

Cmatrix=Table[ToExpression[ToString[Row[{“Ctensor”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

Print[“The Following is the resulting C matrix after adopting all the symmetries:”]

Cmatrix//MatrixForm

(*Orthotropic material*)

Print[“Deriving C tensor and C matrix for Orthotropic Materials”]

(*The following are three reflections around the three planes of symmetry perpendicular to e1=ex and e2=ey and e3=ez*)

Qx=ReflectionMatrix[{1,0,0}];

Qy=ReflectionMatrix[{0,1,0}];

Qz=ReflectionMatrix[{0,0,1}];

(*Using plane perpendicular to x as a symmetry plane*)

Cdashtensor=Table[Sum[Qx[[i,m]]*Qx[[j,n]]*Qx[[k,o]]*Qx[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“The following is the C tensor after reflection around the plane perpendicular to e1=ex:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e1=ex and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e1=ex:”]

ax=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensorx=Ctensor/.ax[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensorx//MatrixForm

(*Using plane perpendicular to y as a symmetry plane*)

Cdashtensor=Table[Sum[Qy[[i,m]]*Qy[[j,n]]*Qy[[k,o]]*Qy[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the C tensor after reflection around the plane perpendicular to e2=ey:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e2=ey and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e2=ey:”]

ay=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensory=Ctensor/.ay[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensory//MatrixForm

(*Using plane perpendicular to z as a symmetry plane*)

Cdashtensor=Table[Sum[Qz[[i,m]]*Qz[[j,n]]*Qz[[k,o]]*Qz[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the C tensor after reflection around the plane perpendicular to e3=ez:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e3=ez and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e3=ez:”]

az=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensorz=Ctensor/.az[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensorz//MatrixForm

CtensorOrthotropic=Ctensor/.Flatten[{az[[1]],ay[[1]],ax[[1]]}];

Print[“This is the resulting C tensor after adopting the orthotropic material symmetries:”]

CtensorOrthotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixOrthotropic=Table[ToExpression[ToString[Row[{“CtensorOrthotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

Print[“Following is the resulting C matrix after adopting the orthotropic material symmetries:”]

CmatrixOrthotropic//MatrixForm

(*Transverse Isotropy*)

Print[“Deriving C tensor of a Transversely Isotropic material from C matrix of Orthotropic Materials assumine that e1=ex and e2=ey form the plane of isotropy”]

(*Rotation around z with 45 Degrees*)

Q=RotationMatrix[45Degree,{0,0,1}];

Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorOrthotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the resulting C tensor after rotating the orthotropic material C matrix with some angle in the plane of isotropy:”]

Cdashtensor=FullSimplify[Cdashtensor];

Cdashtensor//MatrixForm

Print[“Following are the differences between the previous C tensor and the Orthotropic material C tensor:”]

(CtensorOrthotropic-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Vars={};

Do[If[CtensorOrthotropic[[i,j,k,l]]==0,,,AppendTo[Vars,CtensorOrthotropic[[i,j,k,l]]]],{i,3,1,-1},{j,3,1,-1},{k,3,1,-1},{l,3,1,-1}];

Vars;

Print[“Following are the relationships after equating C tensor after rotating in the plane of isotropy to the Ctensor of orthotropic material:”]

a=Solve[(CtensorOrthotropic-Cdashtensor)==zero,Vars] CtensorTransverseIsotropic=CtensorOrthotropic/.a[[1]];

Print[“Following is the C tensor for a Transversely Isotropic Material:”]

CtensorTransverseIsotropic=FullSimplify[CtensorTransverseIsotropic];

CtensorTransverseIsotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixTransverseIsotropic=Table[ToExpression[ToString[Row[{“CtensorTransverseIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

CmatrixTransverseIsotropic=FullSimplify[CmatrixTransverseIsotropic];

Print[“Following is the C matrix for a Transversely Isotropic Material:”]

CmatrixTransverseIsotropic//MatrixForm

(*Full Isotropy*)

Print[“Deriving C tensor of an Isotropic material from the previous Transversely Isotropic Materials C tensor”]

(*Rotation around x with 45 Degrees*)

Q=RotationMatrix[45Degree,{1,0,0}];

Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorTransverseIsotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the resulting C tensor after rotating the Transversely Isotropic material C matrix with some angle around x:”]

Cdashtensor=FullSimplify[Cdashtensor];

Cdashtensor//MatrixForm

Print[“Following are the differences between the previous C tensor and the Transversely Isotropic material C tensor:”]

(CtensorTransverseIsotropic-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following are the relationships after equating C tensor after rotating to the Ctensor of transversely isotropic material:”]

a=Solve[(CtensorTransverseIsotropic-Cdashtensor)==zero,Vars]

CtensorIsotropic=CtensorTransverseIsotropic/.a[[1]];

Print[“Following is the C tensor for an Isotropic Material:”]

CtensorIsotropic=FullSimplify[CtensorIsotropic];

CtensorIsotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixIsotropic=Table[ToExpression[ToString[Row[{“CtensorIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

CmatrixIsotropic=FullSimplify[CmatrixIsotropic];

Print[“Following is the C matrix for an Isotropic Material:”]

CmatrixIsotropic//MatrixForm

CtensorStStrEn=CtensorStStr/.RuleEnergySymmetry;

Print[“The resulting C matrix after using stress, strain and energy symmetries is:”]

CtensorStStrEn//MatrixForm

Print[“The original C tensor and the fully symmetric Ctensor:”]

{Ctensor//MatrixForm,CtensorStStrEn//MatrixForm}

Print[“The difference (minus) between C tensor and the fully symmetric Ctensor:”]

Ctensor-CtensorStStrEn//MatrixForm

(*Adopting the new Ctensor after symmetries*)

Ctensor=CtensorStStrEn;

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

Cmatrix=Table[ToExpression[ToString[Row[{“Ctensor”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

Print[“The Following is the resulting C matrix after adopting all the symmetries:”]

Cmatrix//MatrixForm

(*Orthotropic material*)

Print[“Deriving C tensor and C matrix for Orthotropic Materials”]

(*The following are three reflections around the three planes of symmetry perpendicular to e1=ex and e2=ey and e3=ez*)

Qx=ReflectionMatrix[{1,0,0}];

Qy=ReflectionMatrix[{0,1,0}];

Qz=ReflectionMatrix[{0,0,1}];

(*Using plane perpendicular to x as a symmetry plane*)

Cdashtensor=Table[Sum[Qx[[i,m]]*Qx[[j,n]]*Qx[[k,o]]*Qx[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“The following is the C tensor after reflection around the plane perpendicular to e1=ex:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e1=ex and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e1=ex:”]

ax=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensorx=Ctensor/.ax[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensorx//MatrixForm

(*Using plane perpendicular to y as a symmetry plane*)

Cdashtensor=Table[Sum[Qy[[i,m]]*Qy[[j,n]]*Qy[[k,o]]*Qy[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the C tensor after reflection around the plane perpendicular to e2=ey:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e2=ey and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e2=ey:”]

ay=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensory=Ctensor/.ay[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensory//MatrixForm

(*Using plane perpendicular to z as a symmetry plane*)

Cdashtensor=Table[Sum[Qz[[i,m]]*Qz[[j,n]]*Qz[[k,o]]*Qz[[l,p]]*Ctensor[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the C tensor after reflection around the plane perpendicular to e3=ez:”]

Cdashtensor//MatrixForm

Print[“Following is the difference between the C tensor after reflection around the plane perpendicular to e3=ez and the original symmetric C tensor:”]

(Ctensor-Cdashtensor)//MatrixForm

Print[“Following are the relationships after equating C tensor after reflection to the Ctensor before reflection around the plane perpendicular to e3=ez:”]

az=Solve[(Ctensor-Cdashtensor)==zero,Flatten[Ctensor]]

Ctensorz=Ctensor/.az[[1]];

Print[“Following is the resulting C tensor after adopting the previous relations:”]

Ctensorz//MatrixForm

CtensorOrthotropic=Ctensor/.Flatten[{az[[1]],ay[[1]],ax[[1]]}];

Print[“This is the resulting C tensor after adopting the orthotropic material symmetries:”]

CtensorOrthotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixOrthotropic=Table[ToExpression[ToString[Row[{“CtensorOrthotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

Print[“Following is the resulting C matrix after adopting the orthotropic material symmetries:”]

CmatrixOrthotropic//MatrixForm

(*Transverse Isotropy*)

Print[“Deriving C tensor of a Transversely Isotropic material from C matrix of Orthotropic Materials assumine that e1=ex and e2=ey form the plane of isotropy”]

(*Rotation around z with 45 Degrees*)

Q=RotationMatrix[45Degree,{0,0,1}];

Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorOrthotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the resulting C tensor after rotating the orthotropic material C matrix with some angle in the plane of isotropy:”]

Cdashtensor=FullSimplify[Cdashtensor];

Cdashtensor//MatrixForm

Print[“Following are the differences between the previous C tensor and the Orthotropic material C tensor:”]

(CtensorOrthotropic-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Vars={};

Do[If[CtensorOrthotropic[[i,j,k,l]]==0,,,AppendTo[Vars,CtensorOrthotropic[[i,j,k,l]]]],{i,3,1,-1},{j,3,1,-1},{k,3,1,-1},{l,3,1,-1}];

Vars;

Print[“Following are the relationships after equating C tensor after rotating in the plane of isotropy to the Ctensor of orthotropic material:”]

a=Solve[(CtensorOrthotropic-Cdashtensor)==zero,Vars] CtensorTransverseIsotropic=CtensorOrthotropic/.a[[1]];

Print[“Following is the C tensor for a Transversely Isotropic Material:”]

CtensorTransverseIsotropic=FullSimplify[CtensorTransverseIsotropic];

CtensorTransverseIsotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixTransverseIsotropic=Table[ToExpression[ToString[Row[{“CtensorTransverseIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

CmatrixTransverseIsotropic=FullSimplify[CmatrixTransverseIsotropic];

Print[“Following is the C matrix for a Transversely Isotropic Material:”]

CmatrixTransverseIsotropic//MatrixForm

(*Full Isotropy*)

Print[“Deriving C tensor of an Isotropic material from the previous Transversely Isotropic Materials C tensor”]

(*Rotation around x with 45 Degrees*)

Q=RotationMatrix[45Degree,{1,0,0}];

Cdashtensor=Table[Sum[Q[[i,m]]*Q[[j,n]]*Q[[k,o]]*Q[[l,p]]*CtensorTransverseIsotropic[[m,n,o,p]],{m,1,3},{n,1,3},{o,1,3},{p,1,3}],{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following is the resulting C tensor after rotating the Transversely Isotropic material C matrix with some angle around x:”]

Cdashtensor=FullSimplify[Cdashtensor];

Cdashtensor//MatrixForm

Print[“Following are the differences between the previous C tensor and the Transversely Isotropic material C tensor:”]

(CtensorTransverseIsotropic-Cdashtensor)//MatrixForm

zero=Table[0,{i,1,3},{j,1,3},{k,1,3},{l,1,3}];

Print[“Following are the relationships after equating C tensor after rotating to the Ctensor of transversely isotropic material:”]

a=Solve[(CtensorTransverseIsotropic-Cdashtensor)==zero,Vars]

CtensorIsotropic=CtensorTransverseIsotropic/.a[[1]];

Print[“Following is the C tensor for an Isotropic Material:”]

CtensorIsotropic=FullSimplify[CtensorIsotropic];

CtensorIsotropic//MatrixForm

(*Order of the stress vector*)

order={{1,1},{2,2},{3,3},{1,2},{1,3},{2,3}};

(*Converting Ctensor to Cmatrix*)

CmatrixIsotropic=Table[ToExpression[ToString[Row[{“CtensorIsotropic”,”[[“,Row[Table[If[OddQ[i],Flatten[{order[[j]],order[[k]]}][[(i+1)/2]],”,”],{i,1,7}]],”]]”}]]],{j,1,6},{k,1,6}];

CmatrixIsotropic=FullSimplify[CmatrixIsotropic];

Print[“Following is the C matrix for an Isotropic Material:”]

CmatrixIsotropic//MatrixForm