## Constitutive Laws: Linear Elastic Materials

#### 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 compnents. 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 to state 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

### Restrictions on the Elastic Moduli:

The strain energy density per unit volume for a linear elastic material under a general state of stress is equal to:

As shown above, since each of the stress and strain matrices are symmetric, they can be represented as vectors in . The dot product can then be used to find a simplified expression for the strain energy density function:

We can use two physical restrictions on the materials to impose mathematical restrictions on the entries of the matrix . These two restrictions are that the material is stable and deformable. Stability means that whenever a state of strain is applied, energy is stored in the material. Stability is also equivalent to the statement that the material deformation increases in the direction of the increasing stress. The equivalent mathematical restriction on is:

Since is symmetric, then has to be positive definite. This implies that all its eigenvalues have to be positive. This argument will be used to find restrictions on the values of the elastic moduli described above.

#### Restrictions on the Elastic Moduli of Linear Elastic Isotropic Materials:

The distinct eigenvalues of the matrix in Equation 6 are:

For these to be positive, the following inequalities have to be satisfied:

The following mathematica code presents the calculations above:

View Mathematica Code

Cc = {{1/Ee, -Nu/Ee, -Nu/Ee, 0, 0, 0}, {-Nu/Ee, 1/Ee, -Nu/Ee, 0, 0,0}, {-Nu/Ee, -Nu/Ee, 1/Ee, 0, 0, 0}, {0, 0, 0, 1/G, 0, 0}, {0, 0,0, 0, 1/G, 0}, {0, 0, 0, 0, 0, 1/G}};

Cc // MatrixForm

a = Eigenvalues[Cc] Reduce[{a[[1]]>0, a[[2]]>0}, {Ee, Nu}]

Materials possessing a negative Poisson’s ratio are common, see for example the article here or the Wikipedia article on Auxetics (materials possessing negative Poisson’s ratio).

A more intuitive method to find the same restrictions is by assuming different states of stresses and ensuring that the energy stored per unit volume under these states of stresses is positive. The first state is a uniaxial state of stress where is the only non-zero stress component. Using Equation 6, the stresses and the strains have the following forms:

The restriction that the energy stored has to be greater than zero leads to as follows:

The second state is a shear state of stress where are the only non-zero stress components. Using Equation 6, the stresses and the strains have the following forms:

The restriction that the energy stored has to be greater than zero leads to as follows:

The third state is a state of hydrostatic stress where are the only non-zero stress components. Using Equation 6, the stresses and the strains have the following forms:

The restriction that the energy stored has to be greater than zero leads to as follows:

#### Restrictions on the Elastic Moduli of Linear Elastic Orthotropic Materials:

For orthotropic linear elastic materials, the positive definiteness of imply the following inequalities:

The above equations show that values for the different orthotropic Poisson’s ratios are not bounded by –1 and 0.5. In fact, orthotropic materials (typically, fibrous biological materials like ligaments) can exhibit values beyond those bounds. For example, ligaments can exhibit a Poisson’s ratio that is higher than 0.5.

##### Important Note:

The above restrictions ensure that, physically, a material is stable (possessing a rising stress strain curve) and that any state of stress will cause nonzero strain values inside the material. It is important to realize that it cannot be argued that and/or that are physically impossible. They are simply unacceptable constants for linear elastic materials with a positive definite constitutive model. Materials that become unstable with loading (e.g., decrease in stress accompanying an increase in strain) or incompressible materials (materials exhibiting no strain for a stress state may exist but cannot be modelled using the linear elastic constitutive model

### Plane Isotropic Linear Elastic Materials Constitutive Laws:

In many solid mechanics problems, it is possible, due to loading and geometry symmetries, to simplify particular problems to be solved in rather than . In such cases, the matrix of elastic moduli shown in Equations 6 and 7 can be reduced from a matrix to a or matrix that accounts for such symmetries. In the following, the material is assumed to be isotropic linear elastic and the measure of strain is the infinitesimal strain tensor.

#### Plane Stress:

The plane stress assumption is used when the non-zero stress components are , , and . Usually, this assumption is used for situations when a material is allowed to freely contract and/or expand in one direction (called the thickness direction) while the loads are applied only in the plane perpendicular to the thickness direction. In this situation, the stress matrix has the following form:

Using Equation 6, and the strain matrix has the following form:

The relationship between the stress and strain components shown in Equation 6 can now be reduced to only account for the stresses that are not necessarily equal to zero:

(9)

The inverse of this relationship has the following form:

(10)

In the plane stress assumption, the material is allowed to expand and contract freely in the third direction. The corresponding value of the axial strain in the thickness direction can be obtained using the equation:

#### Plane Strain:

The plane strain assumption is used when the non-zero strain components are , , and . Usually, this assumption is used for situations when a material is restrained from contracting or expanding in one direction (called the thickness direction). The loads are applied in the plane perpendicular to the thickness direction and a restraining stress is applied in the thickness direction. In this situation, the strain matrix has the following form:

Using Equation 7, and the stress matrix has the following form:

The relationship between the stress and strain components shown in Equation 7 can now be reduced to only account for the strains that are not necessarily equal to zero:

(11)

The inverse of this relationship has the following form:

(12)

In the plane strain assumption, the material is restrained from expanding or contracting in the third direction. The restraining axial stress in the thickness direction can be obtained using the equation:

#### Axisymmetry:

Axisymmetry is used for modelling cylindrical shapes when the loads and boundary conditions are also cylindrical. By considering an cross section where is aligned with the axis of symmetry, the shear stresses and can be assumed to be zero. In this situation, the stress matrix has the following form:

Using Equation 6, and the strain matrix has the following form:

The axisymmetry assumption leads to the following relationship between the displacement , the radius and the strain component as follows:

The relationship between the stress and strain components shown in Equation 6 can now be reduced to only account for the stresses that are not necessarily equal to zero:

The inverse of this relationship has the following form:

### Examples and Problems:

#### Example 1:

The state of strain inside a continuum is represented by the small strain matrix:

Assume that the material is steel with Young’s modulus of 210 GPa and Poisson’s ratio of 0.3. Assume that the material follows the von Mises yield criterion with MPa. Find the principal strains and the principal strains directions, the principal stresses and their directions, and comment on whether this state of strain is possible before yielding of the material.

##### Solution:

The principal strains (eigenvalues) are

and the corresponding principal directions (eigenvectors) of the strains are:

To apply the linear elastic constitutive relations, the strain matrix will be converted into a vector of strains, as follows:

Using Equation 7, the corresponding stress vector is:

The corresponding stress matrix is:

The principal stresses (eigenvalues) of the stress matrix are:

the principal directions (the corresponding eigenvalues of the stress matrix) are:

Since the material is an isotropic linear elastic material, the eigenvectors of the stress and the strain matrix coincide. The maximum principal stress is 714.53 MPa; however, the von Mises stress is equal to 1031.19 MPa, which is higher than MPa. This indicates that the material cannot sustain the stress matrix above, and in fact, the strain state above is impossible to occur before yielding or failure of the steel metal.

View Mathematica Codea=Eigensystem[str];

a//MatrixForm

str//MatrixForm

Ee=210000;

Nu=0.3;

G=Ee/2/(1+Nu);

Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};

Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};

Dd=FullSimplify[Inverse[Cc]];

stressvector=Dd.Strainvector

Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};

Stressmatrix//MatrixForm

b=Eigensystem[Stressmatrix];

b//MatrixForm

VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-

sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];

VonMises[Stressmatrix]

#### Example 2:

In a confined compression test of a linear elastic isotropic material, the confining stress was -500MPa while the vertical stress was -100 MPa. The axial strains measured along the directions of the vectors , , and were -0.003, 0.001, and -0.003, respectively. Find Young’s modulus and Poisson’s ratio.

##### Solution:

Using the constitutive relationships for linear elastic isotropic materials:

Solving the two equations:

#### Example 3:

Consider a transversely isotropic linear elastic piece of material with the isotropy plane being the plane of and . The values of Young’s modulus and Poisson’s ratio in the plane of isotropy are equal to 100 MPa and 0.2, respectively. In a uniaxial tension test, when the material was stretched in the direction of by a stress of a value 1 MPa the strain measured in the direction of the vector was recorded to be –0.0005. In another uniaxial tension test, when the material was stretched in the direction of by a stress of a value 1 MPa, the strain measured in the direction of was 0.002. Find an estimate for the material constants , , and .

##### Solution:

The Young’s modulus in the direction of can be obtained from the second experiment:

From the first experiment, the value of Poisson’s ratio can be obtained as follows:

The value of Poisson’s ratio can be obtained from the symmetry of the coefficients matrix:

#### Example 4:

Consider an orthotropic material whose major axes of orthotropy lie in the directions of the vectors:

The material properties in the coordinate systems defined by the orthotropy axes are: , , , , , , , , and . Consider the coordinate system defined by the basis set where

If the material is in a state of uniaxial stress such that the only nonzero stress component is , find the engineering strain matrix in the coordinate systems defined by the basis sets and .

##### Solution:

Since the material is orthotropic, the coefficients matrix depends on the coordinate system involved. The values of the material coefficients are given in the coordinate system defined by . Therefore, it is best to apply the constitutive relationship in the coordinate system of .

The stress matrix in the coordinate system of is given by:

The rotation matrix required for the coordinate transformation is:

The stress matrix in the coordinate system defined by is given by:

The relationship between the vector representations of the stress and the strain in the coordinate system of can be used to find the strains in the coordinate system of (See Equation 2) as follows:

Therefore, the strain matrix in the coordinate system of is given by:

On the other hand, the strain matrix in the coordinate system of is given by:

It is very important to notice that, while a uniaxial state of stress was applied in the coordinate system defined by , shear strain was developed in the same coordinate system. This is because the uniaxial state of stress was applied in a coordinate system different from the orthotropy axes coordinate system!

View Mathematica Code

nu31=0.2;nu32=0.1;nu21=0.05;

nu13=nu31/E33*E11;nu12=nu21/E22*E11;nu23=nu32/E33*E22;

Ss={{1/E11,-nu21/E22,-nu31/E33,0,0,0} ,{-nu12/E11,1/E22,-nu32/E33,0,0,0} ,{-nu13/E11,-nu23/E22,1/E22,0,0,0},{0,0,0,1/G12,0,0},{0,0,0,0,1/G13,0},

{0,0,0,0,0,1/G23}};

s={{1,0,0},{0,0,0},{0,0,0}};

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

sdash=Q.s.Transpose[Q];

sdash//MatrixForm

mat=sdash

sigmavector={mat[[1,1]],mat[[2,2]],mat[[3,3]],mat[[1,2]],mat[[1,3]],mat[[2,3]]}

sigmavector//MatrixForm

strvector=Ss.sigmavector;

strvector//MatrixForm

strdash={{strvector[[1]],strvector[[4]]/2,strvector[[5]]/2},{strvector[[4]]/2,strvector[[2]],strvector[[6]]/2},{strvector[[5]]/2,strvector[[6]]/2,strvector[[3]]}} ;

strdash//MatrixForm

str=Transpose[Q].strdash.Q;

str//MatrixForm

#### Example 5:

The in-plane strain components , , and of a specimen are 0.001, -0.002, and 0.003, respectively. If Young’s modulus is 10,000MPa and Poisson’s ratio is 0.2, find the stress and strain matrices if the material is in the state of plane strain and if the material is in the state of plane stress.

##### Solution:

For the state of plane strain, the strain matrix has the following form:

Equation 11 can be used to find the stress components , , and as follows:

In a plane strain condition, the movement in the direction of is restrained, and thus, a restraining stress component develops according to the relationship:

Thus, the stress matrix has the following form:

View Mathematica Code for the State of Plane Strain:

Cc=Ee/(1+Nu)*{{(1-Nu)/(1-2Nu),Nu/(1-2Nu),0} ,{Nu/(1-2Nu),(1-Nu)/(1-2Nu),0},{0,0,1/2}};

Cc//MatrixForm

strvector={0.001,-0.002,0.003};

stressvector=Cc.strvector;

stressvector//MatrixForm

s=Table[0,{i,1,3},{j,1,3}];

s[[1,1]]=stressvector[[1]];

s[[2,2]]=stressvector[[2]];

s[[1,2]]=stressvector[[3]];

s[[2,1]]=s[[1,2]];

s[[3,3]]=Nu (s[[1,1]]+s[[2,2]]);

s//MatrixForm

In a plane stress condition, a nonzero strain component can develop, and thus, the strain matrix has the form:

The relationship between the stresses and the strains in a plane stress condition (Equation 10) can be utilized to find the values of the non-zero stress components, as follows:

Thus, the stress matrix has the following form:

In a plane stress condition, develops according to the relationship:

View Mathematica Code for the State of Plane Stress:

Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};

Cc//MatrixForm

strvector={0.001,-0.002,0.003};

stressvector=Cc.strvector;

stressvector//MatrixForm

s=Table[0,{i,1,3},{j,1,3}];

s[[1,1]]=stressvector[[1]];

s[[2,2]]=stressvector[[2]];

s[[1,2]]=stressvector[[3]];

s[[2,1]]=s[[1,2]];

s[[3,3]]=0;

s//MatrixForm

str=Table[0,{i,1,3} ,{j,1,3}];

str[[1,1]]=strvector[[1]];

str[[2,2]]=strvector[[2]];

str[[1,2]]=str[[2,1]]=strvector[[3]]/2;

str[[3,3]]=-Nu/Ee*(s[[1,1]]+s[[2,2]]);

str//MatrixForm

#### Example 6:

The position function in dimensions of m of a plate adheres to the following form:

If the plate has dimensions 2m in the direction of , and 1m in the direction of , with the bottom left corner coinciding with the origin of the coordinate system, then:

- Show that this is a state of plane strain.
- Find the engineering small strain matrix, and the Green strain matrix and show that the engineering small strain is adequate to describe the motion.
- Draw the vector plot of the displacement function.
- If Young’s modulus and Poisson’s ratio are equal to 210 GPa and 0.3, respectively, then find the equilibrium body forces to which this plate is subjected.
- Draw the contour plot of the von Mises stress and the normal out-of-plane stress component.

##### Solution:

Unlike the previous examples, this example shows a nonlinear displacement function resulting in a strain field that varies according to position. I.e., every material point will have a different local strain matrix. The displacement function for this plate has the following form:

The following figures shows the vector plot of the displacement function as drawn by Mathematica:

The gradient of the displacement tensor is:

The engineering and the Green strain matrices are:

The difference between the engineering and the Green strain matrices is very small and can be ignored. Thus, the engineering strain can be used to adequately describe the motion (i.e., the displacement is not associated with large rotations). In addition, the displacement component and the other displacement components are independent of the coordinate . Therefore, the strain components , i.e., a state of plane strain.

In order to calculate the stresses, the vector representation of the strain will be used:

Using the constitutive equation for the linear elastic isotropic material (Equation 7), the vector representation of the stress is:

Since the displacement is small, the coordinates can be used instead of the coordinates to find the equilibrium body forces:

The von Mises stress:

The contour plot of the von Mises stress shows that the areas with the maximum von Mises stresses are the right and left bottom corners of the plate. Compare with the vector plot of the displacement (shown above) and note that the maximum displacements are not necessarily correlated with the maximum strain or stress. Since the plate is under a state of plane strain, the restraining stress component is non-zero. The contour plot of is shown below as well (Units are in MPa).

View the Mathematica codex1=X1+0.001*X2+0.0002*X1^2;

x2=1.001*X2-0.0002*X2^2;

x3=X3;

X={X1,X2,X3};

x={x1,x2,x3};

u=x-X;

u//MatrixForm

uplane=Table[u[[i]],{i,1,2}]

VectorPlot[uplane,{X1,0,2},{X2,0,1},AspectRatio->Automatic,PlotLabel->”u”,VectorStyle->Black]

gradu=Table[D[u[[i]],X[[j]]],{i,1,3},{j,1,3}];

gradu//MatrixForm

strGreen=FullSimplify[1/2*(gradu+Transpose[gradu]+Transpose[gradu].gradu)];

Expand[strGreen]//MatrixForm

str=FullSimplify[1/2*(gradu+Transpose[gradu])];

str//MatrixForm

Strainvector={str[[1,1]],str[[2,2]],str[[3,3]],2*str[[1,2]],2*str[[1,3]],2*str[[2,3]]};

Strainvector//MatrixForm

Ee=210000;

Nu=0.3;

G=Ee/2/(1+Nu);

Cc={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};

Dd=FullSimplify[Inverse[Cc]];

stressvector=FullSimplify[Chop[Dd.Strainvector]];

stressvector//MatrixForm

Stressmatrix={{stressvector[[1]],stressvector[[4]],stressvector[[5]]},{stressvector[[4]],stressvector[[2]],stressvector[[6]]},{stressvector[[5]],stressvector[[6]],stressvector[[3]]}};

Stressmatrix=FullSimplify[Chop[Stressmatrix]];

Stressmatrix//MatrixForm

X={X1,X2,X3};

s=Stressmatrix;

-Sum[D[s[[i,1]],X[[i]]],{i,1,3}]

-Sum[D[s[[i,2]],X[[i]]],{i,1,3}]

-Sum[D[s[[i,3]],X[[i]]],{i,1,3}]

VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-

sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];

a=FullSimplify[Chop[VonMises[Stressmatrix]]]

ContourPlot[a,{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”vonMises”,ColorFunction->”GrayTones”] ContourPlot[Stressmatrix[[3,3]],{X1,0,2},{X2,0,1},AspectRatio->Automatic,ContourLabels->All,PlotLabel->”s33″,ColorFunction->”GrayTones”]

#### Example 7:

The in plane displacements of a plate of dimensions of 5m and 1m in the directions of and , respectively, with the bottom left corner situated at the origin are according to the relationship:

If the plate is in the state of plane stress, draw the contour plot of the change of thickness of the plate, assuming that Young’s modulus is 1000 MPa, Poisson’s ratio is 0.4, and the thickness of the plate is 1mm.

##### Solution:

Since the plate is in a state of plane stress, while a strain compoenent develops according to the equation:

The in plane strain components can be obtained from the in plane gradient of the displacement function, as follows:

The stresses can be obtained using the plane stress constitutive relation (Equation 10):

Thus, the strain component is:

The change in thickness in mm can be obtained using the formula:

The required contour plot is shown below:

View Mathematica Code

u={0.02X1+0.02X1*X2,0.0015X2*X1^2+0.02X2};

VectorPlot[u,{X1,0,5},{X2,0,1},AspectRatio->Automatic]

gradu=Table[D[u[[i]],X[[j]]],{i,1,2},{j,1,2}];

gradu//MatrixForm

esmall=1/2*(gradu+Transpose[gradu]);

esmall//MatrixForm

strvector={esmall[[1,1]],esmall[[2,2]],2*esmall[[1,2]]};

strvector//MatrixForm

Ee=1000;Nu=0.4;Cc=Ee/(1+Nu)*{{1/(1-Nu),Nu/(1-Nu),0},{Nu/(1-Nu),1/(1-Nu),0},{0,0,1/2}};

stressvector=Cc.strvector;

stressvector//MatrixForm

strain=Table[0,{i,1,3},{j,1,3}];

strain[[1;;2,1;;2]]=esmall;

strain[[3,3]]=-Nu/Ee*(stressvector[[1]]+stressvector[[2]]);

strain//MatrixForm

thickness=1;

delta=FullSimplify[strain[[3,3]]*thickness]

ContourPlot[delta,{X1,0,5},{X2,0,1},ContourLabels->All,AspectRatio->Automatic,PlotLabel->”delta Thickness mm”]

#### Example 8:

A cylindrical pressure vessel is made out of a steel material. The material’s Young modulus and Poisson’s ratio are 210 GPa and 0.3 respectively. The pressure vessel is used to hold a gas with an internal pressure of 3 MPa. The vessel’s average diameter is 2 m, and its thickness is 10mm. Find the von Mises stress, the change in diameter, and the change in thickness of the cylinder in the following two conditions:

- Condition 1: Capped ends.
- Condition 2: Restrained ends.

##### Solution:

As shown in the sketch, the orthonormal basis vectors are chosen in this example such that is aligned with the longitudinal axis of the pipe and is vertical. The analyzed part of the pressure wall is taken to be

perpendicular to .

For condition 1: The normal stress component in the longitudinal direction is:

The normal stress component in the direction of the vector is:

All the remaining components of the stress matrix are zero. Thus, the stress matrix has the following form:

The von Mises stress is:

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

The strain component gives the relative change of diameter as follows:

Therefore the change in diameter is:

The component provides the relative change in the thickness of the pressure vessel. Therefore:

View Mathematica code for Condition 1:

s=Table[0,{i,1,3},{j,1,3}];

s[[1,1]]=P*r/2/t

s[[2,2]]=P*r/t

VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-

sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];

a=VonMises[s]

G=Ee/2/(1+Nu);

Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};

stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}

strainvector=Ss.stressvector;

strainvector//MatrixForm

deltadiameter=strainvector[[2]]*2*r

deltathickness=strainvector[[3]]*t

For condition 2, the normal stress in the longitudinal direction () is not known since the vessel is restrained in the longitudinal direction. The circumferential stress (normal stress in the direction of the vector ) is given by:

Thus, the stress matrix has the form:

Since for condition 2, the pipe is restrained from expanding or contracting, therefore, . Therefore:

The von Mises stress is thus given by:

The strain-stress relationship (Equation 6) can be used to find the vector representation of the strain:

The strain component gives the relative change of diameter as follows:

Therefore the change in diameter is:

The component provides the relative change in the thickness of the pressure vessel. Therefore:

View Mathematica code for Condition 2:

s=Table[0,{i,1,3},{j,1,3}];

s[[2,2]]=P*r/t

s[[1,1]]=Nu*s[[2,2]]

VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]-sigma[[2,2]])^2+(sigma[[2,2]]-sigma[[3,3]])^2+(sigma[[3,3]]-sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))];

a=VonMises[s]

G=Ee/2/(1+Nu);

Ss={{1/Ee,-Nu/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,1/Ee,-Nu/Ee,0,0,0},{-Nu/Ee,-Nu/Ee,1/Ee,0,0,0},{0,0,0,1/G,0,0},{0,0,0,0,1/G,0},{0,0,0,0,0,1/G}};

stressvector={s[[1,1]],s[[2,2]],s[[3,3]],s[[1,2]],s[[1,3]],s[[2,3]]}

strainvector=Ss.stressvector;

strainvector//MatrixForm

deltadiameter=strainvector[[2]]*2*r

deltathickness=strainvector[[3]]*t

#### Problems:

- What does a negative Poisson’s ratio mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with negative Poisson’s ratio exist (cite any references used).
- What does a Poisson’s ratio>0.5 mean? Conduct a literature search and write a small paragraph supporting or opposing that materials with Poisson’s ratio>0.5 exist (cite any references used).
- Assuming a linear elastic material such that the relationship between the stress and the strain is described by the following matrix:
Assuming that the material is isotropic (i.e., the relationship above applies in any coordinate system) and starting from a state of stress described by the matrix:

show that

(Hint: use a procedure similar to the one shown below the definition of the shear modulus).

- Show the following relationships for linear elastic isotropic materials:
where:

are the components of the deviatoric stress tensor

are the components of the deviatoric strain tensor

and are Lamé’s constants. - A state of stress of a linear isotropic material is represented by a positive definite symmetric stress matrix. Find the volumetric strain in terms of the principal stresses, Young’s modulus and Poisson’s ratio. Comment on whether the volume increases, decreases, stays constant or whether the information provided is not enough to answer this question.
- In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the confining pressure was equal to 50% of the applied vertical compressive stress . Find the relationship between and the volumetric strain in terms of Young’s modulus and Poisson’s ratio. Also, find the maximum shear stress, and the von Mises stress in terms of .
- In a confined compression test of a linear elastic isotropic material, the testing apparatus was such that the specimen was restrained from expanding (or contracting) laterally. Find the relationship between the vertical compressive stress and the vertical strain in terms of Young’s modulus and Poisson’s ratio.
- Find the value of Poisson’s ratio at which both the shear modulus and the bulk modulus are equal to each other in value.
- A cylindrical sample is subjected to a confined compression test in which the horizontal stress is kept constant at –200 MPa, while the magnitude of the compressive vertical stress is increased to reach the value of –600 MPa. Find the volumetric strain if the shear and bulk moduli are 22.72 GPa and 20.833 GPa, respectively. (Answer: –1.6%) .
- A linear elastic isotropic material with Young’s modulus of 210GPa and Poisson’s ratio of 0.3 is in a state of strain described by the infinitesimal strain matrix:
- Find the principal strains and their directions.
- Find the principal stresses and their directions.
- Are the principal stresses and the principal strains aligned?

- Consider a linear elastic orthotropic material in a state of strain described by the infinitesimal strain matrix as in the previous problem. Assume the following material properties:

, , .- Find the principal stresses and their directions.
- Are the principal stresses and the principal strains aligned?

- Using sound mathematical arguments show that the principal stresses and the principal strains are always aligned for linear elastic isotropic materials. (Hint: use the fact that Equations 6 and 7 are valid in any coordinate system)
- The in-plane displacement function of the two dimensional plate shown below is given by:

Assuming that the plate is of a linear elastic isotropic material with Young’s modulus and Poisson’s ratio =0.2, draw the contour plot of the stress components , , , , , and the hydrostatic stress for each of the following two cases:- the plate is in a state of plane stress.
- the plate is in a state of plane strain.

- A cylindrical pressure vessel with an average diameter of 1m and a wall thickiness of 10mm is made of a material with Young’s modulus of 2.1 GPa and Poisson’s ratio 0.3. A constant pressure of 1.3MPa is maintained inside the vessel
- Assuming that the length of the pressure vessel is 5m, find the increase in length. (Answer: 31mm)
- Find the increase in the diameter of the pressure vessel. (Answer: 26.3mm)
- Find the change in thickness of the pressure vessel. (Answer: –0.14mm)
- Find the volumetric strain throughout the pressure vessel wall. (Answer: 1.9%)

- An elastic cylindrical pressure vessel with average diameter , thickness , length , Young’s modulus , Poisson’s ratio is subjected to an internal pressure . Assume that the thickness is very small compared to the diameter. Find the circumferential stress, the longitudinal stress, the circumferential strain, and the longitudinal strain in the following two conditions:
- Condition 1: The pressure vessel has two caps on its ends.
- Condition 2: The pressure vessel is infinitely long and, thus, cannot expand or contract in the longitudinal direction.

### 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 CodeCtensor=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<i,3,j]}];

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