Determination of the contact pressure and plot of the equilibrium shapes of an elastic ring with points of contact

Contents of this notebook :
Part I. Definitions and input data
Part II. Determination of the contact values of the parameters σ, q and η
Part III. Plot of the equilibrium shape
Part IV. Confirmation of the contact conditions -- formulas (45)
r((Subscript[s, 1])) = r(T - (Subscript[s, 1]))	t((Subscript[s, 1])) = - t(T - (Subscript[s, 1]))

Part I. Definitions and input data

Definitions from Section 5. Analytic expressions for the ring curvatures and slope angles

μ=(1/4)((α^2)+(β^2)+(γ^2)+α β+α γ+β γ);

γ=-(( α+β)/2)+I η ; δ=-(( α+β)/2)-I η;

λ=(1/4)(Sqrt[A B]); k=(1/2)-((4 (η^2)+(3 α+β) (α+3 β))/(2 A B));

A=(Sqrt[4 (η^2)+(((3 α+β))^2)]); B=(Sqrt[4 (η^2)+(((α+3 β))^2)]);

T=(4/λ) EllipticK[k];

Left-hand side of the closure condition -- formula (42)

CC=(((A+B) (α-β))/(2 λ (A-B))) EllipticPi[-((((A-B))^2)/(4 A B)),k] +  ((A β-B α)/(λ (A-B))) EllipticK[k];

New parametrization -- formula (44)

α=((4 σ)/((q^2)+(η^2)))-q; β=((4 σ)/((q^2)+(η^2)))+q;

Input data: choose the number of the symmetry axes

n=3;

Part II. Determination of the contact values of the parameters σ, q and η

Contour plots of the three transcendental equations (42), (43) and (60)

F1=ContourPlot3D[CC==(π/(2 n)),{q,2,8},{η,2,8},{σ,10,30},ContourStyle->Blue,Lighting->"Neutral"];

F2=ContourPlot3D[(T/4)==(π/(2 n)),{q,2,8},{η,2,8},{σ,10,30},ContourStyle->Yellow,Lighting->"Neutral"];

φ[s_]= ((A β-B α)/(A-B)) s+(((A+B) (α-β))/(2 λ (A-B))) EllipticPi[-((((A-B))^2)/(4 A B)),JacobiAmplitude[s λ,k],k]+((α-β)/(2 λ (Sqrt[k+((((A-B))^2)/(4 A B))])))  ArcTan[(Sqrt[k+((((A-B))^2)/(4 A B))])((JacobiSN[s λ,k])/(JacobiDN[s λ,k]))];

s1=(1/λ) EllipticF[ArcCos[((A β+B α+ (Sqrt[2 μ]) (A+B))/(A β-B α+ (Sqrt[2 μ]) (A-B)))],k];

F3=ContourPlot3D[φ[s1]==(( π)/n)-(( π)/2),{q,2,8},{η,2,8},{σ,10,30},ContourStyle->Red,Lighting->"Neutral"];

Show[F1, F2, F3]

Identification of the contact values of the parameters σ, q and η

qguess=5;  ηguess=5; σguess = 15;

Comment: The foregoing guess values qguess, ηguess and σguess for q, η and σ, respectively, are checked for n≤10. For other values of n these guess values may not work and should be increased accordingly. New guess values may be identified  converting the Style of the three ContourPlot3D and Show commands above from Text to Input.

FR=FindRoot[{(T/4)==(π/(2 n)),CC==(π/(2 n)), φ[s1]==(( π)/n)-(( π)/2)},{q,qguess},{η,ηguess}, {σ,σguess},AccuracyGoal->60,WorkingPrecision->150,MaxIterations->3000];

FR1=Extract[FR,1];

FR2=Extract[FR,2];

FR3=Extract[FR,3];

q=Abs[q/.FR1];

η=Abs[η/.FR2];

σ=Abs[σ/.FR3];

Part III. Plot of the equilibrium shape

Expressions for the curvature (38) and its derivative

κ[s_]=(((A β+B α)-(A β-B α) JacobiCN[s λ,k])/((A+B)-(A-B) JacobiCN[s λ,k]));

κ1[s_]=(∂_s) κ[s];

Parametric equations (25)

X[s_]=(1/σ) κ1[s] Cos[φ[s]]+(( ((κ[s])^2)-2 μ)/(2 σ)) Sin[φ[s]];

Y[s_]=(1/σ) κ1[s] Sin[φ[s]]-(( ((κ[s])^2)-2 μ)/(2 σ)) Cos[φ[s]];

Plot

ParametricPlot[{X[s],Y[s]},{s,0, n T},AspectRatio->Automatic, Axes -> False,PlotStyle->{Black, Thick}]

Part IV. Confirmation of the contact conditions -- formulas (45)

s2=T-s1;

N[X[s1]-X[s2]]

N[Y[s1]-Y[s2]]

X1[s_]=(∂_s) X[s]; Y1[s_]=(∂_s) Y[s];

N[X1[s1]+X1[s2]]

N[Y1[s1]+Y1[s2]]

N[κ[s1]+(Sqrt[2 μ])]

N[σ,10]