# How to model wooden joints with Mathematica's FEM?

This is a dovetail joint:

and I'd like to see the stresses and deformation on the joint.I haven't seen any modeling of disconnected regions with FEM, only connected regions, so I'm curious if you can do it. I guess you "need" disconnected regions since one would like to study how the fit of the joint determines its strength.

I don't need a full answer, just to know if it's possible to couple the boundary conditions between these two meshes, and how to. The idea being fixing one side on one beam and applying a force on the other beam, when the join is together.

code for the regions:

L = 1
Subscript[\[CapitalGamma], D] =
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, y == -L];
r = RegionUnion[Rectangle[{0, -L}, {10, 1}],
Triangle[{{9, 0}, {11, .75}, {11, -.75}}]];
r2 = RegionDifference[Rectangle[{10, -5}, {12, 5}],
Triangle[{{9, 0}, {11, .75}, {11, -.75}}]];
mesh =
ToElementMesh[r, "MaxCellMeasure" -> 0.05,
MeshQualityGoal -> "Maximal"];
mesh2 =
ToElementMesh[r2, "MaxCellMeasure" -> 0.05,
MeshQualityGoal -> "Maximal"];


Update:

the code of user21 is not exactly what I look for since the mesh regions are continuous. What's the problem with this? Wood joints only work on compression, not on tension, so the results we obtain with this mesh are totally unrealistic.

Note how the lower part of the green beam is pulled by the "endgrain" of the red beam, this cannot happen in reality. This is why I'm asking for discontinuous (disconnected) regions that I haven't seen on the documentation. Maybe another solution is to define lines where only compression happens?

To further elaborate on my comment, one could add a third region "air" around the joints, and give it a very low young modulus. This however will make the meshing far too small when the joints are precisely "cut".

• I could imagine that making a sandwich of 3 materials, the middle one with a very low young modulus could do the trick, but I rather have two distinct entities. Apr 6, 2016 at 5:44
• Very interesting problem. Did you find a solution for this joint-problem with one sided contact? Feb 25 at 10:20
• This is a tricky problem that requires much work. You need to iterate and look for nodes that are in tension not compression. Then you have to re-mesh with a gap at that node and calculate again. You probably need many more nodes than you have here. Similarly, if you have inter-penetration between surfaces you need to re-mesh with connections at that point. There can also be sliding at interfaces so you need to add tangential friction at the surfaces. All this is standard for research in the bolted joint community. Without doubt this is a huge challenge to do properly. Can you approximate?
– Hugh
Mar 8 at 15:29

We can use method described in my answer here. First we define mesh with a small space between regions as follows

L = 1; eps = 1/10;

r1 = RegionUnion[Rectangle[{0, -L}, {10 - eps, 1}],
Triangle[{{9, 0}, {11, .75}, {11, -.75}}]];
r2 = RegionDifference[Rectangle[{10, -5}, {12, 5}],
Triangle[{{9, 0}, {11, .75}, {11, -.75}}]];
mesh1 = ToElementMesh[r, "MaxCellMeasure" -> 0.05,
MeshQualityGoal -> "Maximal"];
mesh2 = ToElementMesh[r2, "MaxCellMeasure" -> 0.05,
MeshQualityGoal -> "Maximal"];
r = RegionUnion[r1, r2]; mesh =
ToElementMesh[r, "MaxCellMeasure" -> 0.05,
MeshQualityGoal -> "Maximal"];

Show[mesh1["Wireframe"["MeshElementStyle" -> FaceForm[Red]]],
mesh2["Wireframe"["MeshElementStyle" -> FaceForm[Green]]]]


Using mesh we can compute plane stress with force applied to the left end and reaction force in small depth

 planeStress = {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}} . Inactive[Grad][v[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}} .
Inactive[Grad][u[y, z], {y, z}], {y, z}],
Inactive[
Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}} . Inactive[Grad][u[y, z], {y, z}], {y, z}] +
Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
0}, {0, -(Y/(1 - \[Nu]^2))}} .
Inactive[Grad][v[y, z], {y, z}], {y, z}]} /. {Y ->
10^3, \[Nu] -> 33/100};

sol = ParametricNDSolveValue[{planeStress == {NeumannValue[-g,
y == 10 - eps && z > 0.5] +
NeumannValue[g, y == 10 && z > 0.5],
NeumannValue[1, z == -L && y < 1]},
DirichletCondition[u[y, z] == 0, z == -5 || z == 5],
DirichletCondition[v[y, z] == 0, z == -5 || z == 5]}, {u,
v}, {y, z} \[Element] mesh, {g}]


Visualization

dmesh = ElementMeshDeformation[mesh, sol[10], "ScalingFactor" -> 1];
dmesh["Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Red], FaceForm[Green]]]]