(*Equilibrium Equations*)
eqn1 = {D[\[Sigma]rD[σr, r] + (\[Sigma]rσr - \[Sigma]\[Theta]σθ)/r +
D[\[Tau]D[τ, z], D[\[Sigma]zD[σz, z] + D[\[Tau]D[τ, r] + \[Tau]τ/r,
1/r D[\[Sigma]\[Theta]D[σθ, r]};
(*Stress Strain*)
eqn2 = {\[Sigma]rσr ->
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]rϵr - \[Nu]ν (\[Epsilon]\[Theta]ϵθ + \[Epsilon]zϵz)), \
\[Sigma]\[Theta]σθ ->
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]\[Theta]ϵθ - \[Nu]ν (\[Epsilon]rϵr + \[Epsilon]zϵz)), \
\[Sigma]zσz ->
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]zϵz - \[Nu]ν (\[Epsilon]rϵr + \[Epsilon]\[Theta]ϵθ)), \
\[Tau]τ -> Y/(2 (1 + \[Nu]ν)) \[Gamma]rzγrz};
(*Strain Displacement*)
eqn3 = {\[Epsilon]rϵr -> D[u[r, z], r], \[Epsilon]zϵz ->
D[w[r, z], z], \[Epsilon]\[Theta]ϵθ -> u[r, z]/r, \[Gamma]rzγrz ->
D[u[r, z], z] + D[w[r, z], r]};
eq = {\[Sigma]rσr =
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]rϵr - \[Nu]ν (\[Epsilon]\[Theta]ϵθ + \[Epsilon]zϵz)), \
\[Sigma]\[Theta]σθ =
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]\[Theta]ϵθ - \[Nu]ν (\[Epsilon]rϵr + \[Epsilon]zϵz)), \
\[Sigma]zσz =
Y/((\[Nu]ν + 1) (2 \[Nu]ν - 1)) ((\[Nu]ν -
1) \[Epsilon]zϵz - \[Nu]ν (\[Epsilon]rϵr + \[Epsilon]\[Theta]ϵθ)), \
\[Tau]τ = Y/(2 (1 + \[Nu]ν)) \[Gamma]rzγrz} /. eqn3;
sys = {D[eq[[1]], r] + (eq[[1]] - eq[[2]])/r + D[eq[[4]], z],
D[eq[[3]], z] + D[eq[[4]], r] + eq[[4]]/r, 1/r D[eq[[2]], r]};
To derive the system of equations we use the code
(*Equilibrium Equations*)
eqn1 = {D[\[Sigma]r, r] + (\[Sigma]r - \[Sigma]\[Theta])/r +
D[\[Tau], z], D[\[Sigma]z, z] + D[\[Tau], r] + \[Tau]/r,
1/r D[\[Sigma]\[Theta], r]};
(*Stress Strain*)
eqn2 = {\[Sigma]r ->
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]r - \[Nu] (\[Epsilon]\[Theta] + \[Epsilon]z)), \
\[Sigma]\[Theta] ->
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]\[Theta] - \[Nu] (\[Epsilon]r + \[Epsilon]z)), \
\[Sigma]z ->
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]z - \[Nu] (\[Epsilon]r + \[Epsilon]\[Theta])), \
\[Tau] -> Y/(2 (1 + \[Nu])) \[Gamma]rz};
(*Strain Displacement*)
eqn3 = {\[Epsilon]r -> D[u[r, z], r], \[Epsilon]z ->
D[w[r, z], z], \[Epsilon]\[Theta] -> u[r, z]/r, \[Gamma]rz ->
D[u[r, z], z] + D[w[r, z], r]};
eq = {\[Sigma]r =
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]r - \[Nu] (\[Epsilon]\[Theta] + \[Epsilon]z)), \
\[Sigma]\[Theta] =
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]\[Theta] - \[Nu] (\[Epsilon]r + \[Epsilon]z)), \
\[Sigma]z =
Y/((\[Nu] + 1) (2 \[Nu] - 1)) ((\[Nu] -
1) \[Epsilon]z - \[Nu] (\[Epsilon]r + \[Epsilon]\[Theta])), \
\[Tau] = Y/(2 (1 + \[Nu])) \[Gamma]rz} /. eqn3;
sys = {D[eq[[1]], r] + (eq[[1]] - eq[[2]])/r + D[eq[[4]], z],
D[eq[[3]], z] + D[eq[[4]], r] + eq[[4]]/r, 1/r D[eq[[2]], r]};
Here we get three equations for two unknowns. The third equation can be integrated independently of the first two. It gives expression to $\sigma _\theta$. Then we can build a solution using FEM. To solve a specific problem, one needs to know the boundary conditions. Here is an example of deformation under compression:
r0 = 1;
r1 = 6;
r2 = 8;
z1 = 4;
z2 = 5; Y = 10^3; \[Nu] = 1/3;
<< NDSolve`FEM`
mesh = ToElementMesh[
RegionUnion[Rectangle[{r0, 0}, {r2, z2}],
Rectangle[{r0, z2}, {r1, z1 + z2}]]];
mesh["Wireframe"]
{ufun, wfun} =
NDSolveValue[{sys[[1]] == NeumannValue[-10, z == z1 + z2],
sys[[2]] == 0,
DirichletCondition[{u[r, z] == 0, w[r, z] == 0}, z == 0]}, {u,
w}, {r, z} \[Element] mesh]
mesh = ufun["ElementMesh"];
Show[{
mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[mesh, {ufun, wfun}][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
