Solving a system of PDEs on a piecewise polynomial domain

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.

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

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

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

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

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

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

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

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

Code for version 10.3

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

• 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. – LejlaS May 29 at 12:31
• @LejlaS You're right, there is a difference between the two solutions. I corrected the code and added Figure 3. – Alex Trounev May 29 at 19:38
• 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] – LejlaS Jun 6 at 9:13
• @LejlaS This is an automatic option InterpolationOrder -> 3 for Interpolation[]. – Alex Trounev Jun 6 at 17:58
• How to get values of displacements at specific nodes? – LejlaS Jun 24 at 11:05