The Geometry of Intersecting Tubes Applied to Controlling a Robotic Welding Torch

John Stockie
Department of Mathematics and Statistics, Simon Fraser University
Jan. 28/98

Figure 1: A photograph of a welding joint between

two cylindrical pipes intersecting at an angle of 90

> restart;

Introduction: A welding problem
Many problems arising in engineering and physics have a mathematical description that is too complicated for consideration in lower-level mathematics courses. It is often possible to make simplifying assumptions that reduce the problem to manageable proportions, but it is often the case that these simplifications are made at the expense of physical significance. Consequently, it is very rewarding to find an application that can serve as a useful illustrative example, while still retaining its full physical interpretation. The situation presented here is one such application.

This is a problem that arose from a consulting contract with the welding industry. The question is this: How does one control a robotic welding torch so that the welder follows precisely the joint between two cylindrical pipes? A typical pipe weld is pictured in the photograph in Fig. 1, where the joint between the two pipes is clearly seen to describe a closed curve. In the most general situation, the company must be able to form welds between pipes with differing cross-sectional radius and arbitrary joint angle.

*** START 1 HERE ***

We begin by defining the parametric equations for the two cylinders:

> cyl1 := { x=R1*cos(theta1),

> y=R1*sin(theta1), z=z1 };

> cyl2 := { x=R2*cos(theta2)*cos(Phi)

> -z2*sin(Phi), y=R2*sin(theta2),

> z=R2*cos(theta2)*sin(Phi)+z2*cos(Phi) };

The curves of intersection are found by equating the $x$, $y$ and $z$ components of each surface:

> eqns := subs( cyl2, cyl1 );

The ensuing equations can then be solved for the unknowns $\theta_1$, $z_1$ and $z_2$ in terms of $\theta_2$ and the cylinder parameters $R_1$, $R_2$ and $\Phi$:

> soln := solve( eqns, {theta1,z1,z2} );

(though we could equally well have solved for \emph{any} of the three unknowns in terms of the fourth).

Notice the presence of the \mcode{RootOf()} expression appearing in the label \mcode{{\%}1}, whose argument is a quadratic polynomial having two roots:

> ztemp := R2*sin(theta2) - 2*R1*_Z

> + R2*sin(theta2)*_Z^2:

> zsoln := solve( ztemp, _Z );

Consequently, there are \emph{two distinct solution curves}, corresponding to the fact that cylinder 2 passes through cylinder 1 in two places (refer to Fig.~\ref{fig:cylsproj}). Looking at the functional form of the roots above, it should now be clear why it was necessary to make the earlier assumption of $R_1 \geq R_2$: without this restriction, the square root terms in \mcode{zsoln} are undefined for some values of $\theta_2$ in the interval $[0,2\pi]$.

Our next step is to derive the parametric equations of the welding joint by taking the expressions in \mcode{soln} and substituting them into Eqs.~(\ref{eq:one}). Not surprisingly, the resulting expressions are quite lengthy and benefit from some simplification. The following function definition will prove very useful in this regard:

> #readlib( rationalize );

> trigsimp := proc( a )

> expand( a, trig ):

> rationalize( % ):

> simplify( %, trig );

> end;

The \mcode{rationalize()} procedure is employed because of the presence of square roots of trig functions in the denominator of the expressions in \mcode{soln}. We now simplify the two welding curves in a \mcode{for} loop:

>

> subs( RootOf(ztemp) = zsoln[i], soln ):

> subs( %, cyl1 ):

> weld := allvalues(trigsimp( % )):

>

The two sets of parametric equations describing the intersection curves on either side of cylinder 1 are

> weld[1];

and

> weld[2];

In practice, a pipe weld consists of a single joint (as pictured in Figs.~\ref{fig:photo} and~\ref{fig:cyls}), and so only one of the two sets of equations (either \mcode{weld[1]} or \mcode{weld[2]}) is required. With this in mind, we next generate a plot of a single pipe weld, consisting of surface plots of cylinder 1 and the ``physical half'' of cylinder 2, and a curve plot of the intersection \mcode{weld[1]}, for representative parameter values $R_1=1$, $R_2=\frac{9}{10}$ and $\Phi=\frac{\pi}{4}$:

> with(plots): #plotsetup(x11):

Warning, the name changecoords has been redefined

> vars := { R1 = 1, R2 = 9/10, Phi = Pi/4 }:

> plot1 := plot3d(

> subs( subs( vars, cyl1 ), [x,y,z] ),

> theta1=0..2*Pi, z1=-2..5 ):

> plot2 := plot3d(

> subs( subs( vars, cyl2 ), [x,y,z] ),

> theta2=0..2*Pi, z2=0..5 ):

> plot3 := spacecurve(subs( subs(vars,weld[2] ),[x,y,z] ),theta2=0..2*Pi, thickness=3, color=red ):

> display( { plot1, plot2, plot3 },

> orientation=[110,60],

> scaling=CONSTRAINED );

>

*** END 1 HERE ***

It is easy to modify the list of parameters in 'vars' above to explore the shape of the intersection curves for other values of the radii and intersection angle.

This code produces plots for the MapleTech paper.

(A) plot of the cylinder 1, half of cylinder 2, and the single weld curve.

> opts := `noborder`:

> ori := [110, 60]:

> #plotsetup( ps, plotoptions=opts,plotoutput=`both.ps` ):

> display( { plot1, plot2, plot3 },

> orientation=ori, scaling=CONSTRAINED );

(B) plot of both full cylinder and the two intersection curves.

> plot5 := plot3d(

> subs( subs( vars, cyl2 ), [x,y,z] ),

> theta2=0..2*Pi, z2=-2..5,

> shading=z ):

> plot4 := spacecurve(

> subs( subs( vars, weld[2] ), [x,y,z] ),

> theta2=0..2*Pi, thickness=3,

> linestyle=2, color=green ):

> #plotsetup( ps, plotoptions=opts,plotoutput=`cyls-bw.ps` ):

> display( {plot1, plot3, plot5, plot4}, orientation=ori, scaling=CONSTRAINED,style=patch );

(C) a plot only of the weld curves (same parameters)

> #plotsetup( ps, plotoptions=opts,plotoutput=`joints.ps` ):

> display( {plot3,plot4}, orientation=ori,scaling=CONSTRAINED );

(D) plots of the weld curves only for the two extreme cases
- equal radii, and
- R1 >> R2.

> rrr := [1, 99/100, 1/3]:

> fff := [`weld1.ps`, `weld2.ps`, `weld3.ps`]:

> ori2 := [ 70, 65 ]:

> for i from 1 to 3 do

> #plotsetup( ps, plotoptions=opts,plotoutput=fff[i] );

> vars78 := {R1=1, R2=rrr[i], Phi=Pi/3};

> plot7[i] := spacecurve(

> subs( subs( vars78, weld[1] ), [x,y,z] ),

> theta2=0..2*Pi, numpoints=500, scaling=CONSTRAINED, color=black ):

> plot8[i] := spacecurve(

> subs( subs( vars78, weld[2] ), [x,y,z] ),

> theta2=0..2*Pi, numpoints=500, scaling=CONSTRAINED, color=black ):

>

> od:

> for i from 1 to 3 do
display( {plot7[i], plot8[i]} ):
od;

>

>

...

*** START 2 (anim) HERE ***

> with(plots): #plotsetup(x11):

> vars3 := {R1=1,R2=1/2}:

> animate3d( { subs( subs( vars3,

> subs( {theta1=theta2, z1=2*z2-5/2},

> cyl1 ) ), [x,y,z] ),

> subs( subs(vars3,cyl2), [x,y,z] ),

> subs( subs(vars3,weld[1]), [x,y,z] ),

> subs( subs(vars3,weld[2]), [x,y,z] )},

> theta2=0..2*Pi, z2=0..5/2,

> Phi=Pi/16..15*Pi/16,

> scaling=CONSTRAINED, frames=12,

> view=[-2..2,-2..2,-3..3] );

Change the parameter R2 above to do the animation with cylinders of differing radius.

Another way to animate is to vary the radius, $R_2$, which we do in the following code:

> vars4 := {R1=1,Phi=Pi/4}:

> animate3d( { subs( subs( vars4,

> subs( {theta1=theta2, z1=2*z2-5/2},

> cyl1 ) ), [x,y,z] ),

> subs( subs(vars4,cyl2), [x,y,z] ),

> subs( subs(vars4,weld[1]), [x,y,z] ),

> subs( subs(vars4,weld[2]), [x,y,z] )},

> theta2=0..2*Pi, z2=0..5/2, R2=1/4..1,

> scaling=CONSTRAINED, frames=12,

> view=[-2..2,-2..2,-3..3] );

>

*** END 2 HERE ***

>

*** START 3 HERE ***

The outward--pointing normal vector to cylinder 1 is

> diff( subs( cyl1, [x,y,z] ), R1 );

which in terms of the variable $\theta_2$, is given by

> n1 :=allvalues(trigsimp( subs( RootOf(ztemp)=zsoln[1], subs( soln, % ) ) )):n1;

The second normal vector is

> n2 := diff( subs( cyl2, [x,y,z] ), R2 );

>

from which we can use the \mcode{linalg[dotprod]} procedure to compute the angle between two cylinders as

> beta := ( Pi - arccos( linalg[dotprod](n1[1], n2) ) ) /2;

>

>

We can then plot a series of torch angle curves using the Maple \mcode{seq()} function:

> #plotsetup(x11):
i := 'i':

> plot( { seq( subs( { R1=1, R2=1/3,

> Phi=i*Pi/12 }, beta ), i=1..6 ) },

> theta2=0..2*Pi, 'beta'=0..Pi/2 );

>

*** END 3 HERE ***

A few plots for the paper from this previous example.

> #plotsetup( ps, plotoptions=`noborder`, plotoutput=`angleeq.ps` ):

> i:='i':

> plot( { seq( subs( {R1=1, R2=1,Phi=i*Pi/12},

> beta ), i=1..6 ) }, theta2=0..2*Pi, 'beta'=0..Pi/2,

> labelfont=[TIMES,ROMAN,24], axesfont=[TIMES,ROMAN,24] );

>

> #plotsetup( ps, plotoptions=`noborder`, plotoutput=`anglelt.ps` ):

> i:='i':

> plot( { seq( subs( {R1=1, R2=1/3,Phi=i*Pi/12},

> beta ), i=1..6 ) }, theta2=0..2*Pi, 'beta'=0..Pi/2,

> labelfont=[TIMES,ROMAN,24], axesfont=[TIMES,ROMAN,24] );

>

>

>

*** START 4 HERE ***

Arclength calculations.

First, differentiate the three components of the weld curve with respect to $\theta_2$

> diff( subs( weld[1], [x,y,z] ), theta2 );

and then set up the integrand for the arclength integral:

> arcInt := sqrt( %[1]^2 + %[2]^2 + %[3]^2 );

This is an integral that cannot be evaluated explicitly in terms of elementary functions. By performing some algebraic manipulations (beyond the capabilities of Maple) it can be re\-written in terms of an \emph{elliptic integrals}~\cite{abramowitz}, which are well-known functions that arise in many physical applications.

We are therefore forced to approximate the integral numerically, using the Maple construct \mcode{evalf(Int(...))}. To this end, we define an ``arclength function'' which depends on $R_2$ and $\Phi$ (with $R_1$ fixed at 1), and whose third argument is the number of digits of accuracy required in the floating point answer:

> arcFunc := (r,p,acc) ->

> Int( subs( {R1=1,R2=r,Phi=p}, arcInt ),theta2=0..2*Pi, acc );

To verify that the arclength formula returns a reasonable result, we perform a calculation with $R_2=\frac{1}{10}$ being significantly smaller than $R_1$:

> evalf( arcFunc( 1/10, Pi/2, 6 ) );

In this situation, the welding curve is approximately circular, and so we expect the arclength to be very close to the circumference of a circle of radius $\frac{1}{10}$. This is very easily tested as follows:

> evalf( 2*Pi/10, 6 );

Next, we produce a surface plot of the arclength as a function of $R_2$ and $\Phi$, which is pictured in Fig.~\ref{fig:arc}.

> #plotsetup(x11);

> plot3d( arcFunc( r2,phi,3 ), r2=1/10..1,

> phi=Pi/20..Pi/2, grid=[10,10],

> style=patchcontour, axes=frame );

The number of digits of accuracy was chosen to be 3, so that the cost of the computation is kept to a minimum while at the same time being sufficiently accurate for plotting purposes:

>

*** END 4 HERE ***

> #plotsetup(ps, plotoptions=opts, plotoutput=`arclength.ps` );

> plot3d( arcFunc( r2, phi, 3 ), r2=1/10..1,

> phi=Pi/20..Pi/2, grid=[10,10], axes=frame, style=patchcontour,

> labelfont=[TIMES,ROMAN,24], axesfont=[TIMES,ROMAN,24] );

>

>

>