移动热源中心旋转变换（子程序数学模型说明）

qx',y',z'=63Qffππafbcexp-3x'af2+y'b2+z'c2,x' > 0qx',y',z'=63Qfrππarbcexp-3x'ar2+y'b2+z'c2,x'< 0

针对此空间椭圆方程，可以利用旋转椭圆参数方程对此方程进行改进。

下面是二维平面，坐标轴转换。

在机械臂打印过程中，因为装有工具头的Kuka机械臂在尖端上，针对同一个点位，有且仅有8种姿态，对应就有8种不同的热源作用方向。

在金属增材制造中，机械臂挤出装置在打印过程中会不断旋转调整角度，怎样在仿真中实现热源的旋转并与机械臂挤出装置作用方向同时一致呢？

为后续机械臂加入传感器测热源做基于物理上的实时电弧增材制造仿真做准备。

我们需要了解的一个最大的问题就是：热源旋转或者行走是x,y,z在变化，不涉及到其他参数，在任何情况下包括接入实时数据和相关的自动化几何线生成代码，都是针对x,y,z做变化，直接放进x,y,z,公式中的其他参数不变（根据自己实际情况）

有了上述思想，我们就可以对我们的热源模型进行改进。

1)基本旋转矩阵

图中显示了相对于x、y和z轴的旋转角度α、β和γ，显示了围绕z轴旋转一个角度γ的情况。通过收集旋转坐标系的标准基向量，我们可以构造绕z轴的旋转矩阵。

 

同理我们有： 

绕x轴旋转：

 

 绕y轴旋转：

 

2)欧拉角a:

欧拉角相对于运动的坐标系旋转一个坐标系是很重要的。因此，一个旋转发生在另一个旋转之前，因此我们从左到右乘以相应的旋转矩阵(后乘法)。例如，当我们绕x轴旋转一个坐标系，然后绕当前的y轴旋转，最后绕当前的z轴旋转，我们得到XYZ顺序的旋转矩阵:

 有了以上的坐标点旋转变化，我们对其参数化转换，如下列旋转2pi-a之后的坐标变化，

倘若我们要对热源模型进行旋转变化，我们需要对进行逆变化。

 以下给出类似思想的Dflux子程序例子：

      INCLUDE 'ABA_PARAM.INC'
      parameter(one=1.d0)
      DIMENSION COORDS(3),FLUX(2),TIME(2)
      CHARACTER*80 SNAME

               real:: PI,fr,ff,af,ar,b,c,v,x1,y1,L1,x2,y2,r,cita1,L2,x3,y3,cita2,L3,t1,t2,t3,Q,cita0,x0,y0,xx0,yy0
	x=COORDS(1)
	y=COORDS(2)
	z=COORDS(3)
C     
	PI=3.1415
C 热源参数
                fr=0.6
                ff=1.4
                af=0.04
                ar=0.08 
                b=0.04 
                c=0.010 
                v=0.0050  
C  L1、L2、L3三段路径长度、t1、t2、t3三段路径的结束时间
                x1=0.0
                y1=0.090
                L1=0.090
                x2=-0.039673
                y2=0.093603
                r=0.02
                cita1 = PI+atan( (y2-y1)/(x2+r) )
                L2 = cita1*r
                x3 = -0.305
                y3 = 0.15
                cita2=PI+atan( (y3-y2)/(x3-x2) )
                L3=sqrt( (x3-x2)**2+(y3-y2)**2 )
                t1=L1/v
                t2=(L1+L2)/v 
                t3=(L1+L2+L3)/v      

C     JLTYP＝1，表示为体热源 Q热源功率
C	JLTYP=1
               Q=0.8*30*250               

C 以时间作为判定，cita是热源局部坐标系与整体坐标系的夹角，x0,y0是热源中心坐标；xx0,yy0是将整体坐标系的坐标变换到局部坐标系
  
                if ( TIME(2).LE.t1 ) then
                   cita0 = 0                 
                   x0 = 0  
                   y0 = v*TIME(2)
                   xx0 =  (x-x0)*cos(cita0)+(y-y0)*sin(cita0) 
                   yy0 =-(x-x0)*sin(cita0)+(y-y0)*cos(cita0) 
                     if ( yy0.LE.0) then
                       FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*fr*Q*exp(-3*( xx0/b)**2-3*( yy0/ar)**2-3*( z/c)**2)
                   else
                       FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*ff*Q*exp(-3*( xx0/b)**2-3*( yy0/af)**2-3*( z/c )**2)
                    endif

                  elseif ( TIME(2).GT.t1.and.TIME(2).LE.t2) then
                     cita0 = v*( TIME(2) - t1)/r
                     x0 = r*( -1+cos(cita0) )  
                       y0 = y1+r*sin(cita0)
                       xx0 =  (x-x0)*cos(cita0)+(y-y0)*sin(cita0) 
                     yy0 =-(x-x0)*sin(cita0)+(y-y0)*cos(cita0)                    
                       if ( yy0.LE.0) then
                         FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*fr*Q*exp(-3*( xx0/b)**2-3*( yy0/ar)**2-3*( z/c)**2)
                       else
                         FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*ff*Q*exp(-3*( xx0/b)**2-3*( yy0/af)**2-3*( z/c)**2)
                       endif
    
                  elseif ( TIME(2).GT.t2.and.TIME(2).LE.t3) then
                       cita0 = cita2
                       x0 = x2+v*(TIME(2)-t2)*cos(cita0)
                       y0 = y2+v*(TIME(2)-t2)*sin(cita0)  
                       xx0 =  (x-x0)*cos(cita0)+(y-y0)*sin(cita0) 
                       yy0 =-(x-x0)*sin(cita0)+(y-y0)*cos(cita0)          
                       if ( yy0.LE.0) then
                         FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*fr*Q*exp(-3*( xx0/b)**2-3*( yy0/ar)**2-3*( z/c)**2)
                       else
                         FLUX(1) = 6*sqrt(3.0)/(PI*sqrt(PI)*ar*b*c)*ff*Q*exp(-3*( xx0/b)**2-3*( yy0/af)**2-3*( z/c)**2)
                       endif
               
                endif

      RETURN
      END