I have a structural mechanical problem, where i try to calculate the strain field in the modelled material. I simplified the problem in the following example: There is a material with two phases, where the modelled area is periodical in all direction.

`Needs["NDSolve`FEM`"] spher[{x0_,y0_,z0_},r_] := (x - x0)^2 + (y - y0)^2 + (z - z0)^2 >= r^2 Drm = BoundaryDiscretizeRegion[ImplicitRegion[spher[{5, 5, 5}, 4], {x, y, z}], {{0, 10}, {0, 10}, {0, 10}},MaxCellMeasure -> 9.9] boundaryMarkerFunction = Compile[{{boundaryElementCoords, _Real, 3},{pointMarkres, _Integer,2}}, Module[{pt1 = #[[1]], pt2 = #[[2]], pt3 = #[[3]]}, Which[Abs[-5 + pt1[[1]]] < 4.1 && Abs[-5 + pt1[[2]]] < 4.1 && Abs[-5 + pt1[[3]]] < 4.1, 1, True, 2]] & /@ boundaryElementCoords]; bmesh = ToBoundaryMesh[Drm, "BoundaryMarkerFunction" -> boundaryMarkerFunction, MaxCellMeasure -> 1]; model = ToElementMesh[Drm, MaxCellMeasure -> .015, "RegionHoles" -> None, "BoundaryMarkerFunction" -> boundaryMarkerFunction, "MeshOrder" -> 1, "RegionMarker" -> {{5, 5, 5}, 1}] `

The mechanical properties are described by the Navier-Cauchy operator with the proper constants:

`Y = If[ElementMarker == 0, 70, 120]; v = If[ElementMarker == 0, 0.25, 0.2]; op3D := -(v Y /((1 + v) (1 - 2 v)) + Y/(2 (1 + v)) ) \!\(\*SubscriptBox[\(\[Del]\), \({x, y, z}\)]\(\*SubscriptBox[\(\[Del]\), \({x, y, z}\)] . {u[x, y, z], v[x, y, z],w[x, y, z]}\)\) - Y/(2 (1 + v)) \!\(\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\({u[x, y, z],v[x, y, z], w[x, y, z]}\)\); `

In the next step there are given the Neumann conditions for the interface of the two phases:

`Subscript[gamma, NX1] := NeumannValue[Sign[5 - x]*Sin[ArcTan[(z - 5)/(Sign[x - 5]*Sqrt[((x - 5))^2 + (y - 5)^2] +0.001)] -54/180*Pi*ArcTan[Tan[54/180*Pi]/Sqrt[1 + ((y - 5)/(Abs[(x - 5)] + 0.001))^2]]],ElementMarker == 1]; Subscript[gamma, NZ1] := NeumannValue[Sign[x - 5]*Sin[ArcTan[(z - 5)/(Sign[x - 5]*Sqrt[((x - 5))^2 + (y - 5)^2] +0.001)] -54/180*Pi*ArcTan[Tan[54/180*Pi]/Sqrt[1 + ((y - 5)/(Abs[(x - 5)] + 0.001))^2]]],ElementMarker == 1]; `

These are plane forces acting to the boundary of the two phases. Just for illustration I draw it with vectorplot:

As it is visible forces don’t have the same symmetry than the material. Therefore i can’t apply periodic boundary condition for the displacement field (u[x,y,z]). Due to the periodicity of material the strain must be the same on the opposite sides, which means I should define periodic boundary conditions on the derivation of displacement field {u,v,w}. I tried it with DirichletCondition and also with the following way:

`Subscript[gamma, DX2] = PeriodicBoundaryCondition[Derivative[1, 0, 0][u][x,y, z], x == 10, TranslationTransform[{-10, 0, 0}]] `

None of them works because mathematica requires linear equations instead of derivations. Without boundary conditions the result looks the following way:

`{ufun, vfun, wfun} = NDSolveValue[{op3D == {Subscript[gamma, NX1], 0, Subscript[gamma, NZ1]}}, {u, v, w}, {x, y, z} \[Element] model]; epsilon11 = Derivative[1, 0, 0][ufun]; `

This is obviously not realistic for a periodic material because the strain is different in opposite sides. Can someone help me how to handle these periodic conditions?