# Stress analysis in a disk with circular holes

Writing the following code:

Rext = 16;
Rint = 3;
thickness = 1/10;

nholes = 8;
distcenters = 9;
Rholes = 2;

density[x_, y_, z_] := (x^2 + y^2 + z^2) 8 10^-3;

domain = "-thickness/2\[LessEqual]z\[LessEqual]thickness/2" <>
"&&(x-0)^2+(y-0)^2\[LessEqual]Rext^2" <>
"&&(x-0)^2+(y-0)^2\[GreaterEqual]Rint^2";
Table[domain = domain <> "&&(x-distcenters Cos[2(" <> ToString[i] <>
"-1)Pi/nholes])^2+(y-distcenters Sin[2(" <> ToString[i] <>
"-1)Pi/nholes])^2\[GreaterEqual]Rholes^2", {i, 1, nholes}];

disk = ImplicitRegion[ToExpression[domain], {x, y, z}];
min = NMinValue[{density[x, y, z], {x, y, z} \[Element] disk}, {x, y, z}];
max = NMaxValue[{density[x, y, z], {x, y, z} \[Element] disk}, {x, y, z}];

Magnify[RegionPlot3D[disk,
ColorFunction -> Function[{x, y, z},
ColorData["SandyTerrain"][density[x, y, z]]],
ColorFunctionScaling -> False,
PlotLegends -> BarLegend[{"SandyTerrain", {min, max}}],
PlotPoints -> 75], 1.5]

Iz = NIntegrate[(x^2 + y^2) density[x, y, z], {x, y, z} \[Element] disk]


I get: 13474.5

which is wonderfully as desired.

Now, assuming the elastic theory is valid and assuming that this disk is subject to a perpendicular load, I would like to calculate the stress state to which it is subject by plotting as above. For this purpose I tried the following approach:

NDSolveValue[{
D[w[x, y, z], {x, 4}] + 2 D[w[x, y, z], {x, 2}, {y, 2}] + D[w[x, y, z], {y, 4}] == 1,
DirichletCondition[w[x, y, z] == 0, x^2 + y^2 == Rint^2],
DirichletCondition[D[w[x, y, z], {x, 2}] == 0, x^2 + y^2 == Rint^2],
DirichletCondition[D[w[x, y, z], {y, 2}] == 0, x^2 + y^2 == Rint^2]},
w, {x, y, z} \[Element] disk]


but unfortunately I get the following error:

NDSolveValue::femcmsd: The spatial derivative order of the PDE may not exceed two.

Of course before opening this topic for many hours I read other topics and the most similar is this. Could someone help me adapt it to my case? Thank you!

• Currently, the FEM can not handle higher than second derivatives. Having said that, you could try to rewrite your equation as a system of two second order equations. – user21 Dec 28 '17 at 9:16