# Formulating equations for 3D stress in the finite element method

I would like to know how you formulate equations for the finite element method for stress calculations. We know the answer because user21 has put it here. It involves usage of Inactive in stopping the equations from simplifying too early. I am almost there with formulating the equations but need a straightforward method because I need to look at variants.

The general equations for the finite element method must be put into the form See the second on Systems of Partial Differential Equations

here u, v and w are the displacements in the three directions. The $$c_{i j}$$ are 3 by 3 matrices which we have to find. The other coefficients are vectors which we don't need for a stress formulation.

From a standard engineering textbook we find that the strain is related to the stress by

strainStress = {
ϵx[x, y, z] ==
1/Y (σx[x, y,
z] - ν (σy[x, y, z] + σz[x, y, z])),
γxy[x, y, z] == 1/G τxy[x, y, z],
γxz[x, y, z] == 1/G τxz[x, y, z],

γyx[x, y, z] == 1/G τyx[x, y, z],
ϵy[x, y, z] ==
1/Y (σy[x, y,
z] - ν (σx[x, y, z] + σz[x, y, z])),
γyz[x, y, z] == 1/G τyz[x, y, z],

γzx[x, y, z] == 1/G τzx[x, y, z],
γzy[x, y, z] == 1/G τzy[x, y, z],
ϵz[x, y, z] ==
1/Y (σz[x, y,
z] - ν (σx[x, y, z] + σy[x, y, z]))
} /. G -> Y/(2 (1 + ν));


Here σx[x, y, z], σy[x, y, z] and σz[x, y, z] are the normal stresses and τxy[x, y, z], τxz[x, y, z] and τyx[x, y, z] are shear stresses. I have also put in shear stresses τyx[x, y, z], τzx[x, y, z] and τxy[x, y, z] which I might have been able to avoid since there is symmetry. Y is Young's modulus and G is the shear modulus. ν is Poisson's ratio. I have substituted for the shear modulus because it is related to Young's modulus.

The relationship between strain and stress is given by

strainDisp={Derivative[1, 0, 0][u][x, y, z], Derivative[0, 1, 0][u][x, y, z] +
Derivative[1, 0, 0][v][x, y, z], Derivative[0, 0, 1][u][x, y, z] +
Derivative[1, 0, 0][w][x, y, z], Derivative[0, 1, 0][u][x, y, z] +
Derivative[1, 0, 0][v][x, y, z], Derivative[0, 1, 0][v][x, y, z],
Derivative[0, 0, 1][v][x, y, z] + Derivative[0, 1, 0][w][x, y, z],
Derivative[0, 0, 1][u][x, y, z] + Derivative[1, 0, 0][w][x, y, z],
Derivative[0, 0, 1][v][x, y, z] + Derivative[0, 1, 0][w][x, y, z],
Derivative[0, 0, 1][w][x, y, z]};


Thus the job we have before us is to find a matrix that relates stress and displacements and then put it in the form above. We can do this as follows: First find a matrix that relates stress and strains. Then find a second matrix that relates strains and displacements and then multiply the two together. I do this as follows:

1. Make a list of the stress symbols and strain symbols

2. Solve for the stresses and then use CoefficientArrays to get a matrix (mat) for the stresses in terms of the strains.

3. Multiply the mat by the strains to get equations for the stresses in terms of the displacement derivatives.

4. Make a list of displacement derivatives.

5. Use CoefficientArrays again to get a matrix (mat2) for the stresses in terms of the displacement derivatives.

6. Partition mat2 to get the matrices that are the $$c_{i j}$$ we need.

 stress = {
σx[x, y, z], τxy[x, y, z], τxz[x, y, z],
τyx[x, y, z], σy[x, y, z], τyz[x, y, z],
τzx[x, y, z], τzy[x, y, z], σz[x, y, z]
};
strain = strainStress[[All, 1]];
sol = stress /. Solve[strainStress, stress];
{vec, mat} =
CoefficientArrays[Thread[stress == sol[]], strain] // Normal;
dd = Join[Grad[u[x, y, z], {x, y, z}], Grad[v[x, y, z], {x, y, z}],
Grad[w[x, y, z], {x, y, z}]];
eq2 = -mat.strainDisp;
{vec2, mat2} =
CoefficientArrays[Thread[stress == eq2], dd] // Normal // Simplify;
mat3 = Partition[mat2 // Simplify // Factor, {3, 3}];


Is there a better way to go from textbook equations to here? Let me know if there is.

Now we have to put in lots of Inactive to get the equation we need for the finite element formulation in Mathematica. We need an Inactive[Div] and an Inactive[Grad] in lots of places. This is where I get lost and need help. My idea was to formulate this as in the expression above using symbols for the $$c_{i j}$$. Having got this formulation, in Inactive form I could then substitute in my 3 by 3 matrices for each $$c_{i j}$$.

Thus next I

1. Make a list of coefficients for the $$c_{i j}$$, and a list of Inactive Grads.

2. Make an expression for the $$c_{i j}$$ and the Grads. This is half of what I need.

For the next step I am lost. I need to premultiply each $$c_{i j}$$ by an Inactive Div. Then I would have my expression. How do I do this?

cc = Array[c, {3, 3}];
coffs = Thread[Flatten[cc, 1] -> Flatten[mat3, 1]];
idd = {Inactive[Grad][u[x, y, z], {x, y, z}],
Inactive[Grad][v[x, y, z], {x, y, z}],
Inactive[Grad][w[x, y, z], {x, y, z}]}
e1 = cc.idd // ExpandAll • I derived these with pen and paper; I have a very rough notes notebook, that I can send you if you ping me at ruebenko AT wolfram.com – user21 May 31 '19 at 12:45