2
$\begingroup$

I wish to solve following system of equations with Dirchlet and Neumann boundary conditions on an piecewise polynomial (cubic spline) shaped domain as showed here. enter image description here

Emod = 2*10^6; \[Nu] = 0.3; \[Rho] = 7860; g = 10; F = 100000; Gmod = 769231;
System = {Emod/((1 + \[Nu]) (1 - 
     2 \[Nu])) ((1 - \[Nu]) (D[r*U[r, z], r, r])/r - \[Nu]*
   D[r* V[r, z], r, z]/r) +(Emod/(2 (1 + \[Nu]))) (D[U[r, z], z, z] + D[V[r, z], r, z]),(Emod/(2 (1 + \[Nu]))) (D[r* U[r, z], r, z]/ r + 
D[r* V[r, z], r, r]/r) + (Emod/((1 + \[Nu]) (1 - 2 \[Nu]))) ((1 - \[Nu]) D[V[r, z], z, z] - \[Nu]*D[U[r, z], r, z])};
{uif, vif} = NDSolveValue[{System == {0, NeumannValue[F, z == 4]}, DirichletCondition[{U[r, z] == 0., V[r, z] == 0.}, z == 0]}, {U, V}, {r, 0, 1}, {z, 0, 4}];

I know how to solve and plot it in a rectangular domain: {r, 0, 1}, {z, 0, 4}. I'd also like to solve this on the region, where {z, 0, 4} and r(z) is defined as a cubic spline function: $r_{i}=a_{i}+b_{i}(z-z_{i})+c_{i}(z-z_{i})^2+d_{i}(z-z_{i})^3$ for an arbitrary number of intervals.

This cubic spline is interpolation of an exponential function: $r(z)=r_0 e^\frac{\rho r_0^2 gz\pi}{2F}$, where $r_0=0.5$. The function coefficients $(a_{i},..d_{i})$ depend on the number of intervals. The coefficients are (for 10 intervals):

a={0.5, 0.571811, 0.653935, 0.747854, 0.855263, 0.978097, 1.11857, 1.27922, 1.46295};

b={0.15105, 0.172638, 0.197461, 0.225813, 0.258247, 0.295335, 0.337756, 0.386248, 0.441788};

c={0.0224672, 0.0261063, 0.0297453, 0.0340466, 0.0389299, 0.0445182, 0.0509301, 0.0581749, 0.0667915};

d={0.00272929, 0.00272929, 0.00322598, 0.00366244, 0.0041912, 0.00480898, 0.00543356, 0.0064625, 0.0064625}.

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You did not define $a_i,b_i,c_i,d_i,z_i$. $\endgroup$ – Alex Trounev May 24 at 14:17
  • $\begingroup$ They depend on number of intervals, and are obtained by interpolation of exponential function., r(z). Should I also write that? $\endgroup$ – LejlaS May 24 at 14:27
  • $\begingroup$ Yes, add the equation r=r[z] $\endgroup$ – Alex Trounev May 24 at 14:29
3
$\begingroup$

At the beginning, we solve equations in the domain of a bounded curve r=r0*Exp[k*z]

<< NDSolve`FEM`
r0 = .5; Emod = 

2*10^6; \[Nu] = 0.3; \[Rho] = 7860; g = 10; u = 0.1; Gmod = 769231; \
g = 10; F = 100000; k = \[Rho]*g*.5^2*Pi/2/F;
reg = ImplicitRegion[0 <= r <= r0*Exp[k*z] && 0 <= z <= 4, {r, z}];
mesh = ToElementMesh[reg];
mesh["Wireframe"]

System = {Emod/((1 + \[Nu]) (1 - 
          2 \[Nu])) ((1 - \[Nu]) (D[r*U[r, z], r, r])/r - \[Nu]*
        D[r*V[r, z], r, z]/r) + (Emod/(2 (1 + \[Nu]))) (D[U[r, z], z, 
        z] + D[V[r, z], r, 
        z]), (Emod/(2 (1 + \[Nu]))) (D[r*U[r, z], r, z]/r + 
       D[r*V[r, z], r, r]/
        r) + (Emod/((1 + \[Nu]) (1 - 2 \[Nu]))) ((1 - \[Nu]) D[
         V[r, z], z, z] - \[Nu]*D[U[r, z], r, z])};
{uif, vif} = 
  NDSolveValue[{System == {0, 
      NeumannValue[10000., r > 0]}, 
    DirichletCondition[{U[r, z] == 0., V[r, z] == 0.}, z == 4]}, {U, 
    V}, {r, z} \[Element] mesh];
mesh = uif["ElementMesh"];
    {Show[{
   mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
   ElementMeshDeformation[mesh, {uif, vif}][
    "Wireframe"[
     "ElementMeshDirective" -> 
      Directive[EdgeForm[Red], FaceForm[]]]]}], 
 ContourPlot[uif[r, z], {r, z} \[Element] mesh, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "u", Contours -> 20, 
  PlotLegends -> Automatic], 
 ContourPlot[vif[r, z], {r, z} \[Element] mesh, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "v", Contours -> 20, 
  PlotLegends -> Automatic]}

fig1

Now we take nint points on the curve and construct an interpolation by third order polynomials

    nint = 20; point = Table[{z, r0*Exp[k*z]}, {z, 0, 4, 4/nint}]; f = 
 Interpolation[point];

Build a new mesh and solve the problem.

reg1 = ImplicitRegion[0 <= r <= f[z] && 0 <= z <= 4, {r, z}];
mesh1 = ToElementMesh[reg1];
mesh1["Wireframe"]
{uif1, vif1} = 
  NDSolveValue[{System == {0, NeumannValue[10000., r > 0]}, 
    DirichletCondition[{U[r, z] == 0., V[r, z] == 0.}, z == 4]}, {U, 
    V}, {r, z} \[Element] mesh1];

mesh1 = uif1["ElementMesh"];
{Show[{
   mesh1["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
   ElementMeshDeformation[mesh1, {uif1, vif1}][
    "Wireframe"[
     "ElementMeshDirective" -> 
      Directive[EdgeForm[Red], FaceForm[]]]]}], 
 ContourPlot[uif1[r, z], {r, z} \[Element] mesh1, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "u", Contours -> 20, 
  PlotLegends -> Automatic], 
 ContourPlot[vif1[r, z], {r, z} \[Element] mesh1, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "v", Contours -> 20, 
  PlotLegends -> Automatic]} 

Fugure 2

Then we can determine the difference of the two solutions when r = 0 as

{Row[{"nint = ", nint}], 
 Plot[{uif[0, z] - uif1[0, z]}, {z, 0, 4}, 
  AxesLabel -> {"z", "\[Delta]u"}, PlotLegends -> Automatic], 
 Plot[{vif[0, z] - vif1[0, z]}, {z, 0, 4}, 
  AxesLabel -> {"z", "\[Delta]v"}, PlotLegends -> Automatic]}

Figure 3

If coefficients a, b, c, d are given, then the solution to the problem is as follows

<< NDSolve`FEM`
a = {0.5, 0.571811, 0.653935, 0.747854, 0.855263, 0.978097, 1.11857, 
   1.27922, 1.46295};
b = {0.15105, 0.172638, 0.197461, 0.225813, 0.258247, 0.295335, 
   0.337756, 0.386248, 0.441788};
c = {0.0224672, 0.0261063, 0.0297453, 0.0340466, 0.0389299, 0.0445182,
    0.0509301, 0.0581749, 0.0667915};
d = {0.00272929, 0.00272929, 0.00322598, 0.00366244, 0.0041912, 
   0.00480898, 0.00543356, 0.0064625, 0.0064625};
zz = Range[0, 4, 4/9];
f[z_] = Piecewise[
   Table[{a[[i]] + b[[i]] (z - zz[[i]]) + c[[i]] (z - zz[[i]])^2 + 
      d[[i]] (z - zz[[i]])^3, zz[[i]] <= z < zz[[i + 1]]}, {i, 
     Length[a]}]];

reg = ImplicitRegion[0 <= r <= f[z] && 0 <= z <= 4, {r, z}];
mesh = DiscretizeRegion[reg, AccuracyGoal -> 10, 
  MaxCellMeasure -> .001]

r0 = .5; Emod = 
 2*10^6; \[Nu] = 0.3; \[Rho] = 7860; g = 10; u = 0.1; Gmod = 769231; \
g = 10; F = 100000; k = \[Rho]*g*.5^2*
  Pi/2/F; System = {Emod/((1 + \[Nu]) (1 - 
         2 \[Nu])) ((1 - \[Nu]) (D[r*U[r, z], r, r])/r - \[Nu]*
       D[r*V[r, z], r, z]/r) + (Emod/(2 (1 + \[Nu]))) (D[U[r, z], z, 
       z] + D[V[r, z], r, 
       z]), (Emod/(2 (1 + \[Nu]))) (D[r*U[r, z], r, z]/r + 
      D[r*V[r, z], r, r]/
       r) + (Emod/((1 + \[Nu]) (1 - 2 \[Nu]))) ((1 - \[Nu]) D[V[r, z],
         z, z] - \[Nu]*D[U[r, z], r, z])};
{uif1, vif1} = 
 NDSolveValue[{System == {0, NeumannValue[10000., r > 0]}, 
   DirichletCondition[{U[r, z] == 0., V[r, z] == 0.}, z == 4]}, {U, 
   V}, {r, z} \[Element] mesh]



meshu = uif1["ElementMesh"];
{Show[{meshu["Wireframe"["MeshElement" -> "BoundaryElements"]], 
   ElementMeshDeformation[meshu, {uif1, vif1}][
    "Wireframe"[
     "ElementMeshDirective" -> 
      Directive[EdgeForm[Red], FaceForm[]]]]}], 
 ContourPlot[uif1[r, z], {r, z} \[Element] meshu, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "u", Contours -> 20, 
  PlotLegends -> Automatic], 
 ContourPlot[vif1[r, z], {r, z} \[Element] meshu, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "v", Contours -> 20, 
  PlotLegends -> Automatic]}

Figure 4

Code for version 10.3

<< NDSolve`FEM`
a = {0.5, 0.571811, 0.653935, 0.747854, 0.855263, 0.978097, 1.11857, 
   1.27922, 1.46295};
b = {0.15105, 0.172638, 0.197461, 0.225813, 0.258247, 0.295335, 
   0.337756, 0.386248, 0.441788};
c = {0.0224672, 0.0261063, 0.0297453, 0.0340466, 0.0389299, 0.0445182,
    0.0509301, 0.0581749, 0.0667915};
d = {0.00272929, 0.00272929, 0.00322598, 0.00366244, 0.0041912, 
   0.00480898, 0.00543356, 0.0064625, 0.0064625};
zz = Range[0, 4, 4/9];
f[z_] := Piecewise[
   Table[{a[[i]] + b[[i]] (z - zz[[i]]) + c[[i]] (z - zz[[i]])^2 + 
      d[[i]] (z - zz[[i]])^3, zz[[i]] <= z <= zz[[i + 1]]}, {i, 
     Length[a]}]];

reg = ImplicitRegion[0 <= r <= f[z] && 0 <= z <= 4, {r, z}];
mesh = DiscretizeRegion[reg, {{0, f[4]}, {0, 4}}, 
  MaxCellMeasure -> .001, Frame -> True, Method -> "Continuation"]
r0 = .5; Emod = 
 2*10^6; \[Nu] = 0.3; \[Rho] = 7860; g = 10; u = 0.1; Gmod = 769231; \
g = 10; F = 100000; k = \[Rho]*g*.5^2*
  Pi/2/F; System = {Emod/((1 + \[Nu]) (1 - 
         2 \[Nu])) ((1 - \[Nu]) (D[r*U[r, z], r, r])/r - \[Nu]*
       D[r*V[r, z], r, z]/r) + (Emod/(2 (1 + \[Nu]))) (D[U[r, z], z, 
       z] + D[V[r, z], r, 
       z]), (Emod/(2 (1 + \[Nu]))) (D[r*U[r, z], r, z]/r + 
      D[r*V[r, z], r, r]/
       r) + (Emod/((1 + \[Nu]) (1 - 2 \[Nu]))) ((1 - \[Nu]) D[V[r, z],
         z, z] - \[Nu]*D[U[r, z], r, z])};
{uif1, vif1} = 
 NDSolveValue[{System == {0, NeumannValue[10000., r > 0]}, 
   DirichletCondition[{U[r, z] == 0., V[r, z] == 0.}, z == 4]}, {U, 
   V}, {r, z} \[Element] mesh]



meshu = uif1["ElementMesh"];
{Show[{meshu["Wireframe"["MeshElement" -> "BoundaryElements"]], 
   ElementMeshDeformation[meshu, {uif1, vif1}][
    "Wireframe"[
     "ElementMeshDirective" -> 
      Directive[EdgeForm[Red], FaceForm[]]]]}], 
 ContourPlot[uif1[r, z], {r, z} \[Element] meshu, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "u", Contours -> 20, 
  PlotLegends -> Automatic], 
 ContourPlot[vif1[r, z], {r, z} \[Element] meshu, PlotRange -> All, 
  AspectRatio -> Automatic, ColorFunction -> Hue, 
  FrameLabel -> {"r", "z"}, PlotLabel -> "v", Contours -> 20, 
  PlotLegends -> Automatic]}

Figure 5

$\endgroup$
  • $\begingroup$ This looks good, i just find one mistake in your code, mesh1 = ToElementMesh[reg1]; instead of [reg}, then i get different results than with the exponential curve, but i think that this is normal because the discrete points are different. Mathematica then shows the error: InterpolatingFunction::dmval: Input value {4.00028} lies outside the range of data in the interpolating function. Extrapolation will be used. $\endgroup$ – LejlaS May 29 at 12:31
  • $\begingroup$ @LejlaS You're right, there is a difference between the two solutions. I corrected the code and added Figure 3. $\endgroup$ – Alex Trounev May 29 at 19:38
  • $\begingroup$ Thank you for your answer. How do we know its third order polynomial interpolation? I don't see it in your code. its normaly with [InterpolationOrder -> 3] $\endgroup$ – LejlaS Jun 6 at 9:13
  • $\begingroup$ @LejlaS This is an automatic option InterpolationOrder -> 3 for Interpolation[]. $\endgroup$ – Alex Trounev Jun 6 at 17:58
  • $\begingroup$ How to get values of displacements at specific nodes? $\endgroup$ – LejlaS Jun 24 at 11:05

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.