Surface Envelopes for Engineering Applications

Stephan G Bush, P.E.

Principal Engineer, The Procter & Gamble Company

Import[FileNameJoin[{NotebookDirectory[],"gearset.png"}]]

Abstract

Applications of envelopes of manifolds (especially surfaces) are presented in the context of mechanical engineering.  The mathematical analog of cutting away material to create a useful shape, envelopes find myriad applications where analytically determined, physically derived shapes are required.  The generation and exchange of surfaces of envelopes takes on greater importance as computer numerically controlled (CNC) machine tools are able to create surfaces directly from computer aided drawing (CAD) data, and as many mathematically precise surfaces cannot be generated directly by CAD tools.

To enable a complete application, three processes must work in harmony:

Generation

A mathematically precise description of the surface, preferably derived from conceptually simple building blocks.

Conversion

Conversion from exact form to an inexact form, compatible with solid modeling tools.

Integration

Parsing to a form that enables exchange of data between applications.

Motivation

Envelopes

Animate[
Block[
{f,g,h,w,a=2,b=1,c=1,d=1,e=2,constraint},
f=Function[{u},{a Sin[u],b Cos[u]}];
g=Function[{v},{c+d v Sin[2 v],e (v^2)}];
h=Function[{u,v},TranslationTransform[g[v]][f[u]]];
constraint=Solve[
Det[D[h[u,v],{{u,v}}]]==0,u
][[{2,3}]];
w=Function[{v},
Evaluate[
ReplaceAll[
h[u,v],
constraint
]
]
];
Show[
ParametricPlot[g[v],{v,0,3},PlotStyle->Dashed],
ParametricPlot[w[v],{v,0,3}],
ParametricPlot[h[u,V],{u,0,2π}],
PlotRange->{{-5,5},{-1,20}},
ImageSize->200
]
],
{V,0,3}
]

Envelopes in Engineering

f=Function[{θ,z},{r Cos[θ],r Sin[θ],z}]; (* generating shape *)
g=Function[{u},{a u,a (u^2),a u}]; (* path *)
h=Function[{θ,z,u},TranslationTransform[g[u]][f[θ,z]]];(* family 
of surfaces *)

q=Function[{z,u},
Evaluate[
h[θ,z,u]/.
Part[
Simplify@Solve[
Det[D[h[θ,z,u],{{θ,z,u}}]]==0,θ],
{2,3}
]
]
];

Manipulate[
Block[{a=10,r=2},
Show[
ParametricPlot3D[g[u],{u,0,1}],
ParametricPlot3D[h[θ,z,U],{θ,0,2π},{z,0,4},PlotStyle->Directive[Orange,Specularity[White,40]]],
ParametricPlot3D[q[z,u],{u,0,1},{z,0,4}],
PlotRange->{{-5,12},{-5,12},{0,15}},
ImageSize->500
]
],
{{U,0.5},0,1}
]

Generation

Governing Equation

Planar Case

For a planar curve:

f(y,t)
((∂f(y,t))/(∂t))==0

For a parametrically defined curve:

F=={x(u,t),y(u,t)}
(Det[(({{((∂x(u,t))/(∂u)), ((∂x(u,t))/(∂t))}, {((∂y(u,t))/(∂u)), ((∂y(u,t))/(∂t))}}))])==0

General Case

F=={(Subscript[f, 1])[(Subscript[u, 1]),...,(Subscript[u, m])],...,(Subscript[f, n])[(Subscript[u, 1]),...,(Subscript[u, m])]}
Det[Ja[F]]==0

In the typical case, this looks like:

F(θ,z,ϕ)   ,    (Det[((∂(F(θ,z,ϕ)))/(∂{θ,z,ϕ}))])==0

Cam Surfaces

The generating function:

f=Function[{x,θ},{x,a Cos[θ],b Sin[θ]}];

Block[{a=2,b=1},
ParametricPlot3D[f[x,θ],{x,0,10},{θ,0,2π},Axes->False]
]

The transformation is what makes the envelope technique work.  We can describe motion of nearly any complexity in a declarative way through the use of TransformationFunction[].

Tf1=Function[{ϕ},
Evaluate[Composition[
RotationTransform[ϕ,{0,0,1}],
TranslationTransform[{0,0,c Cos[ϕ]}]
]
]
];

Block[
{a=2,b=1,c=3,expr},
expr=Cases[ParametricPlot3D[f[x,θ],{x,0,10},{θ,0,2π}],_GraphicsComplex,Infinity];
Animate[
Evaluate[
Graphics3D[{
Triad,
GeometricTransformation[Triad,Tf1[ϕ]],
GeometricTransformation[expr,Tf1[ϕ]]
},
PlotRange->{{-10,10},{-10,10},{-10,10}},
ImageSize->500
]
],
{ϕ,0,2π}
]
]

The family of surfaces:

g=Function[{x,θ,ϕ},
Evaluate[Tf1[ϕ][f[x,θ]]]
];

The envelope surface.  Note we solve the constraint explicitly, and substitute into the parametric expression for the family of surfaces.  There are multiple solutions, which we could explore in more detail through the use of Reduce[].

h=Function[{x,ϕ},Evaluate[
ReplaceAll[
g[x,θ,ϕ],
Flatten@Solve[
Det[D[g[x,θ,ϕ],{{θ,ϕ,x}}]]==0,
θ
]
]
]
];

Block[
{a=2,b=1,c=3,expr},
expr=Cases[ParametricPlot3D[f[x,θ],{x,0,10},{θ,0,2π}],_GraphicsComplex,Infinity];
Animate[
Evaluate[
Show[
ParametricPlot3D[h[x,ϕ],{x,5,10},{ϕ,0,2π},PlotStyle->Directive[Orange,Specularity[White,20]]],
Graphics3D[{
Triad,
GeometricTransformation[Triad,Tf1[ϕ]],
GeometricTransformation[expr,Tf1[ϕ]]
},
],
PlotRange->{{-10,10},{-10,10},{-10,10}},
ImageSize->500
]
],
{ϕ,0,2π}
]
]

Prescribed Motion

The generating shape:

f=Function[{θ,z},{a Cos[θ],b Sin[θ],z}];

Block[
{a=2,b=1,c=4,d=1,e=0.2},
ParametricPlot3D[f[θ,z],{θ,0,2π},{z,-1,1}]
]

In this case, we will be describing the coordinated motion of two bodies, so there will be two distinct transformations that ultimately will be composited.  In particular, body 1 will translate and body 2 will rotate.

The transformations:

Tf1=Function[{ϕ},
Composition[
TranslationTransform[{0,c,0}],
TranslationTransform[{d ϕ+e (ϕ^2),0,0}]
]
];
Tf2=Function[{ϕ},
Composition[
RotationTransform[ϕ,{1,0,0}]
]
];
T12=Function[{ϕ},
Composition[
Tf2[ϕ],
InverseFunction[Tf1[ϕ]]
]
];
T2f=Function[{ϕ},InverseFunction[Tf2[ϕ]]];
T1f=Function[{ϕ},InverseFunction[Tf1[ϕ]]];

Illustrating the generator and transformations from the body 2 perspective:

Block[
{a=2,b=1,c=4,d=1,e=0.2,surf},
surf=Cases[ParametricPlot3D[f[θ,z],{θ,0,2π},{z,-1,1}],_GraphicsComplex];
Animate[
Evaluate[
Graphics3D[
{
GeometricTransformation[surf,T12[ϕ]],
Triad
},
PlotRange->{{-20,2},{-6,6},{-6,6}},
ImageSize->700
]
],
{ϕ,0,2π}
]
]

Solving for the generated surface:

g=Function[{θ,z,ϕ},Evaluate[T12[ϕ][f[θ,z]]]];

constraint=Solve[Det[D[T12[ϕ][f[θ,z]],{{θ,ϕ,z}}]]==0,θ][[2]];

h=Function[{z,ϕ},g[θ,z,ϕ]/.constraint];

The end result, in the fixed frame of reference:

Block[{a=2,b=1,c=4,d=1,e=0.2,s1,s2},
s1=Cases[ParametricPlot3D[f[θ,z],{θ,0,2π},{z,-1,1}],_GraphicsComplex];
s2=Cases[ParametricPlot3D[h[z,ϕ],{ϕ,0,4π},{z,-2,2},PlotStyle->Directive[Orange,Specularity[White,20]]],_GraphicsComplex];
Animate[
Evaluate[
Show[
Graphics3D[GeometricTransformation[s1,T1f[u]]],
Graphics3D[GeometricTransformation[s2,T2f[u]]],
PlotRange->All,ImageSize->7.5*72
]
],
{{u,π},0,4π}
]
]

The above may represent a bottle being moved along a manufacturing line.  Such a conveyance is called a screw conveyor, and so has considerable application in the real world.  This example shows how powerful this technique is, where a very complex shape can be expressed in just a few lines of code.

Block[{a=2,b=1,c=4,d=1,e=0.2,s1,s2},
ConveyorFunction=Function[{z,ϕ},Evaluate[h[z,ϕ]]]
];

Gear Surfaces

For these examples, we will be using a simple gear, a spur gear, to generate the surfaces for a type of right angle drive called a face gear.  In the first example, we consider the simple case.  In the second case, we consider one example among more exotic variations, which allow us to operate a different angles, have stronger teeth, or a less noisy meshing.

Explicit

The generating surface will be an involute curve, extended in the z-direction.  The function is nested, with parameters (Subscript[r, b]) and ϵ as arguments to make the solution more generic.  This is the equation of the tooth surface of a spur gear.

f=Function[{rb,ϵ},
Function[{θ,z},
{rb (-θ Cos[ϵ+θ]+Sin[ϵ+θ]),rb (Cos[ϵ+θ]+θ Sin[ϵ+θ]),z}
]
];

Here, we need to use the geometry to define reasonable limits ((Subscript[θ, i]) and (Subscript[θ, f])) for the parameter θ.  The reason for this will become apparent when we plot the generating surface.

θi=0;
θf=Replace[θ,
Last[
Block[
{rb=1,ϵ=π/12,z=0},
NMinimize[
{
Abs[Sqrt[Dot[f[rb,ϵ][θ,z],f[rb,ϵ][θ,z]]]-1.2rb],
θ>0.2
},θ]
]
]
]

Block[
{rb=1,ϵ=π/12},
Show[
ParametricPlot3D[f[rb,ϵ][θ,z],{θ,θi,θf},{z,0,1}],
Graphics3D[Triad],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,0},{0,0,1}},rb]}],
PlotRange->All,
Axes->None,
ImageSize->500
]
]

When we generate a gear, we need to consider the motion of two bodies, the generating one and the generated one.  Here are the transforms needed for the motion of the two bodies:

T1f=Function[{ϕ,γ},
Evaluate[
Composition[
RotationTransform[-γ,{1,0,0}],RotationTransform[ϕ,{0,0,1}]
]
]
];
Tf1=Function[{ϕ,γ},Evaluate[Simplify[InverseFunction[T1f[ϕ,γ]]]]];
T2f=Function[{ϕ,m},
Evaluate[
Composition[
RotationTransform[-m ϕ,{0,0,1}]
]
]
];
Tf2=Function[{ϕ,m},Evaluate[Simplify[InverseFunction[T2f[ϕ,m]]]]];
T12=Function[{ϕ,γ,m},Evaluate[Simplify[Composition[Tf2[ϕ,m],T1f[ϕ,γ]]]]];
T21=Function[{ϕ,γ,m},Evaluate[Simplify[InverseFunction[T12[ϕ,γ,m]]]]];

The transformations introduce the additional parameter ϕ, so we estimate some reasonable bounds:

ϕi=-ArcCos[1/1.2];
ϕf=ArcCos[1/1.2];

Manipulate[
Block[{m=10/28,rb=1,ϵ=π/12,c1,c2},
c2={
Opacity[0.5],
Cylinder[{{0,0,-rb},{0,0,-rb-1}},rb/m]
};
c1={
Opacity[0.5],
Cylinder[{{0,0,-rb-1},{0,0,-rb-2}},rb],
Opacity[1],Cases[ParametricPlot3D[f[rb,ϵ][θ,z]+{0,0,-rb-2},{θ,θi,θf},{z,0,1}],_GraphicsComplex,Infinity],
};
Graphics3D[
{
GeometricTransformation[Triad,T1f[ϕ,γ]],
GeometricTransformation[Triad,T2f[ϕ,m]],
EdgeForm[Dashed],
GeometricTransformation[c1,T1f[ϕ,γ]],
GeometricTransformation[c2,T2f[ϕ,m]],
Triad
},
PlotRange->{{-3,3},{-4,4},{-3,3}},
ImageSize->600
]
],
{ϕ,0,2π},
{{γ,π/2},0,π/2}
]

Utilizing the generating surface and the composition of the transforms, we can produce the family of surfaces in the coordinate system of body 2:

g=Function[{rb,γ,ϵ,m},
Function[{θ,z,ϕ},
T12[ϕ,γ,m][f[rb,ϵ][θ,z]]
]
];

Finally, we produce the generated surface by solving the constraint and eliminating one parameter from the equation of the family of surfaces.  Note that in this case, once again, the constraint can be solved explicitly.

h=Function[{rb,γ,ϵ,m},
Evaluate@Function[{θ,ϕ},
Evaluate@Simplify@Block[
{constraint},
constraint=Flatten@Solve[
Simplify@Det[
D[g[rb,γ,ϵ,m][θ,z,ϕ],{{θ,ϕ,z}}]
]==0,
z
];
ReplaceAll[g[rb,γ,ϵ,m][θ,z,ϕ],constraint] 
]
]
];

{xi,xf}=Block[{values},
values=Block[{rb=1,m=12/28,ϵ=π/12,γ=π/2},Flatten@Table[h[rb,γ,ϵ,m][θ,ϕ][[2]],{θ,θi,θf,(θf-θi)/10},{ϕ,ϕi,ϕf,(ϕf-ϕi)/10}]];
{Min[values],Max[values]}
]

Animate[
Block[{m=12/28,γ=90Degree,s1,s2},
s1=Block[{rb=1,ϵ=π/12},Cases[ParametricPlot3D[f[rb,ϵ][θ,z],{θ,θi,θf},{z,xi,xf},MeshStyle->None],_GraphicsComplex]];
s2=Block[{rb=1,ϵ=π/12},Cases[ParametricPlot3D[h[rb,γ,ϵ,m][θ,ϕ],{θ,θi,θf},{ϕ,ϕi,ϕf},PlotStyle->Directive[Orange,Specularity[White,20]],MeshStyle->None],_GraphicsComplex]];
Graphics3D[
{
GeometricTransformation[Triad,T1f[ϕ,γ]],
GeometricTransformation[Triad,T2f[ϕ,m]],
Triad, 
GeometricTransformation[s1,T1f[ϕ,γ]],
GeometricTransformation[s2,T2f[ϕ,m]]
},
PlotRange->{{-3,3},{-6,1},{-2,2}},
ImageSize->700
]
],
{ϕ,ϕi,ϕf}
]

Import[FileNameJoin[{NotebookDirectory[],"tooth flank assembly.tif"}]]

Implicit

What happens if the constraint cannot be solved explicitly?

We will look at the case of generating a face gear surface again, but this time we introduce another parameter that makes the input surface skewed, so that the input body is a helical spur gear.

f=Function[{rb,ϵ,η},Function[{θ,z},{rb (-θ Cos[ϵ-z η+θ]+Sin[ϵ-z η+θ]),rb (Cos[ϵ-z η+θ]+θ Sin[ϵ-z η+θ]),z}]];

Block[
{rb=1,ϵ=π/12,η=15Degree},
Show[
ParametricPlot3D[f[rb,ϵ,η][θ,z],{θ,θi,θf},{z,0,1}],
Graphics3D[Triad],
Graphics3D[{Opacity[0.5],Cylinder[{{0,0,0},{0,0,1}},rb]}],
PlotRange->All,
Axes->None,
ImageSize->500
]
]

Note that the transformations remain the same.

g=Function[{rb,γ,ϵ,m,η},
Function[{θ,z,ϕ},
T12[ϕ,γ,m][f[rb,ϵ,η][θ,z]]
]
];

If we try to solve the constraint, we will fail to find a solution:

Block[
{constraint},
constraint=
Solve[
Simplify@Det[
D[g[rb,γ,ϵ,m,η][θ,z,ϕ],{{θ,ϕ,z}}]
]==0,
z
]
]

Instead of direct substitution, we use a numerical technique to find a substitution value.  Note that this yields only one solution, not multiple branches, and is dependent on the quality of the intial guess and the ease of finding the root.

w=Function[{rb,γ,ϵ,m,η},
Evaluate@Module[{expr}, (* need to evaluate this expensive expression, and localize it *)
expr=Function[{θ,ϕ},
Evaluate[
Abs@Simplify@Det[D[g[rb,γ,ϵ,m,η][θ,z,ϕ],{{θ,ϕ,z}}]]
]
];
Function[{θ,ϕ},Last@Quiet@FindMinimum[expr[θ,ϕ],{z,0},Method->Automatic]]
]
];

Block[{rb=1,m=12/28,ϵ=π/12,γ=π/2,η=2Degree,zmin,zmax,pts,sol,sols,sel},
ListPointPlot3D[
Table[
{θ,ϕ,z/.Last@w[rb,γ,ϵ,m,η][θ,ϕ]},
{θ,θi,θf,(θf-θi)/50},{ϕ,ϕi,ϕf,(ϕf-ϕi)/50}
],
AxesLabel->{"θ","ϕ","z"},
PlotRange->All,
ImageSize->600,
PlotLabel->Style[Framed["Constraint Solutions for z"]]
]
]

pts=Block[
{rb=1,m=12/28,ϵ=π/12,γ=π/2,η=2Degree},
Table[
ReplaceAll[
g[rb,γ,ϵ,m,η][θ,z,ϕ],
w[rb,γ,ϵ,m,η][θ,ϕ]
],
{θ,θi,θf,(θf-θi)/20},
{ϕ,ϕi,ϕf,(ϕf-ϕi)/20}
]
];

ListPointPlot3D[pts,PlotRange->{Automatic,{-5,-2},Automatic},BoxRatios->Automatic,ImageSize->500]

Conversion

It’s a NURBS World Out There

Virtually all computer aided design tools are based on NURBS representations of surfaces.  Therefore, for our generated surfaces to be useful in the manufacturing world, we must translate our algebraic surfaces into NURBS with an adequate degree of accuracy.

pts=Block[
{points},
points=RandomReal[{0,1},{10,2}];
Part[points,Last@FindShortestTour[points]]
];

Show[
ListPlot[pts,PlotStyle->Directive[PointSize[Large],Red]],
Graphics[{BSplineCurve[pts]}]
]

Additional functionality is needed to interpolate points using NURBS functions, so we use a package written for this purpose:

Get[FileNameJoin[{NotebookDirectory[],"..","NURBS","NURBS Trunk Development.m"}]]

For simple functions, interpolation is well behaved:

Block[
{points,cpts},
points=Table[{u,N@Sin[u]},{u,0,2π,π/5}];
cpts=SplineInterpolate[points,Parameterization->"Centripetal"];
Show[
ListPlot[points,Axes->None,PlotStyle->Directive[PointSize[Medium]]],
Graphics[BSplineCurve[cpts]],
ListPlot[cpts,PlotStyle->Directive[PointSize[Large],Red]],
ImageSize->700,
PlotRange->All
]
]

However, for more difficult functions, like those with cusps, interpolation requires a lot more work:

Block[
{points,cpts,a=-1.05,b=0.65},
points=Table[{u,N@Abs[(u^(1/3))]},{u,a,b,0.07}];
cpts=SplineInterpolate[points,Parameterization->"Centripetal"];
Show[
Plot[Abs[(u^(1/3))],{u,a,b},PlotStyle->Gray],
ListPlot[points,Axes->None,PlotStyle->Directive[PointSize[Medium]]],
Graphics[BSplineCurve[cpts]],
ListPlot[cpts,PlotStyle->Directive[PointSize[Large],Red]],
PlotRange->All,ImageSize->700
]
]

Sampling

We sample the exact surface at regular intervals, with the goal of keeping the fidelity of the shape:

pts=Block[
{a=2,b=1,c=4,d=1,e=0.2,surf},Table[ConveyorFunction[z,ϕ],{z,-1,1,2/10},{ϕ,0,2π,2π/10}]
];

ListPointPlot3D[pts,BoxRatios->Automatic,ImageSize->800]

Interpolation

We interpolate a NURBS surface:

cpts=SurfaceInterpolate[pts,Parameterization->"Centripetal"];

Block[
{a=2,b=1,c=4,d=1,e=0.2},
Show[
ParametricPlot3D[ConveyorFunction[z,ϕ],{z,-1,1},{ϕ,0,2π},PlotStyle->Directive[Blue,Specularity[White,20],Opacity[0.1]]],ListPointPlot3D[cpts,BoxRatios->Automatic,PlotStyle->Black],
ImageSize->800,
PlotLabel->Style[Framed["Control Points and Exact Surface"]]
]
]

Troubleshooting Interpolation

Interpolation to recover a NURBS surface from an array of points is no picnic.  Where interpolation fails, it is sometimes helpful to do the following:

Change the parameterization

Interpolate a smaller segment of the surface, and reassemble later

Try approximation instead of interpolation

Error Control

Plotting the NURBS surface and the sampled points, we can get an idea of the goodness of fit:

Show[
Graphics3D[BSplineSurface[cpts]],
ListPointPlot3D[pts,PlotStyle->Directive[Red]],
ImageSize->800,
PlotLabel->Style[Framed["NURBS Surface and Sample Points"]]
]

For a better estimate of goodness of fit, it is easy to oversample the exact surface function and the NURBS function, so the fidelity of the NURBS representation can be judged.

Integration

Preparation for Export

Get[FileNameJoin[{NotebookDirectory[],"..","NURBS","IGES Trunk Development.m"}]]

s=BSplineFunction[cpts];

ParametricPlot3D[s[u,v],{u,0,1},{v,0,1}]

Exporting to IGES Format

Export[
FileNameJoin[{NotebookDirectory[],"surface.iges"}],
s,
"IGES"
]

References

Piegl, Les A. and Tiller, Wayne.  The NURBS Book.  Springer, 1997.

Utility

SetOptions[Solve,InverseFunctions->True];

SetOptions[ParametricPlot3D,Mesh->None,PlotStyle->Directive[Blue,Specularity[White,20]]];

Diad=Block[
{origin,u1,u2},
origin={0,0};
u1={1,0,0};
u2={0,1,0};
{
Arrowheads[0.02],
Arrow[{origin,1.0u1}],
Arrow[{origin,1.0u2}],
Text[Style["X",FontFamily->"Calibri"],1.1 u1],
Text[Style["Y",FontFamily->"Calibri"],1.1 u2]
}
];

Triad=Block[
{origin,u1,u2,u3},
origin={0,0,0};
u1={1,0,0};
u2={0,1,0};
u3={0,0,1};
{
Arrowheads[0.02],
Arrow[{origin,u1}],
Arrow[{origin,u2}],
Arrow[{origin,u3}],
Text[Style["X",FontFamily->"Calibri"],1.1 u1],
Text[Style["Y",FontFamily->"Calibri"],1.1 u2],
Text[Style["Z",FontFamily->"Calibri"],1.1 u3]
}
];