7
$\begingroup$

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}]

Manipulate[model["Wireframe"["MeshElement" -> "MeshElements", "MeshElementStyle" -> {Directive[FaceForm[Green]], Directive[FaceForm[Red]]}, PlotRange -> {All, {i, 10.1}, All}]], {{i, 5.0}, 0, 10}]

enter image description here

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

Y = If[ElementMarker == 0, 70, 120];
vi = If[ElementMarker == 0, 0.25, 0.2];
op3D := -(vi Y /((1 + vi) (1 - 2 vi)) + Y/(2 (1 + vi)) ) \!\(\*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 + vi))  \!\(\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\({u[x, y, z],v[x, y, z], w[x, y, z]}\)\);

enter image description here

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:

Show[{Graphics[{Circle[{5, 5}, 4]}], 
VectorPlot[{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]]],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]]]}, {x, 0.1, 10}, {z, 0, 10},RegionFunction ->Function[{x, z, vx, vz, n}, 14 < (x - 5)^2 + (z - 5)^2 < 22]]}]

enter image description here

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];
DensityPlot[(epsilon11[x, 5, z])*1000000, {x, 0, 10}, {z, 0, 10}, Frame -> True, ColorFunction -> "TemperatureMap", PlotPoints -> 100, PlotLegends -> BarLegend[Automatic, LegendLabel -> "strain_xx [\[Mu]strain]"]]

enter image description here

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? I read that for time derivative it works so I courious if it can work for spacial derivative too or not.

$\endgroup$
  • $\begingroup$ When I run your code I get an error messages: NDSolveValue::dvnoarg: The function If[ElementMarker==0,0.25,0.2] appears with no arguments. Also, could you include the code for the plots? $\endgroup$ – user21 Dec 28 '17 at 8:41
  • $\begingroup$ I uploaded the code for plots too. Yes there was a mistake because Poisson constant had the same mark then the displacement function v[x,y,z]. I corrected it. $\endgroup$ – Gergi Jan 8 '18 at 13:04
  • $\begingroup$ Have you tried to make an Anti-periodic boundary condition on u? $\endgroup$ – user21 Jan 11 '18 at 10:20
  • $\begingroup$ You mean like this: PeriodicBoundaryCondition[-u[x,y, z], x == 10, TranslationTransform[{-10, 0, 0}]] ? It could not work because the u on the two opposite boundary plane dont have mirror symetry. Actually it should be used like this: PeriodicBoundaryCondition[u[x,y, z]+c, x == 10, TranslationTransform[{-10, 0, 0}]] which would give correct results but the constant c have to be calculated analytically. In my case that is quite hard task because elastic properties of the phases are not even isotropic (this example just for representation). $\endgroup$ – Gergi Jan 13 '18 at 19:43
  • $\begingroup$ Unfortunately, I do not understand what you are trying to do. Perhaps you can simplify the problem further? $\endgroup$ – user21 Jan 15 '18 at 6:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.