marc焊接仿真二次开发-以几何位置激活生死单元的方法？ 50

二次开发指南上和网上关于生死单元的二次开发uactive（）大多都是用单元号作的。这样在遍历单元格时很容易判断所遍历的单元生死状态，但对模型分网的网格单元号是有要求的。二次开发指南上网格单元号沿焊接路径方向有规律增加。可是对于不在marc上建模分网的焊接模型，网格号很难满足上述要求。

此时想用遍历单元与移动的热源中心相对位置来决定是否激活单元。我们看到二次开发指南中将生死单元状态另写入以文件中，仿照此方法，在初始化网格生死状态后，将所有会出现生死状态变化的单元的单元号和生死状态写入一.txt文件。进入焊接工况下，为判断遍历单元生死状态以决定生死操作，先将.txt文件读入二维向量，进行生死操作后删除原.txt文件，再将更新后的二维数组数据重新写入.txt文件。

这种方法效率极低，运行极慢。各位仿真大佬有没有效率更高的方法，十分感谢。另附现有代码。

另外，marc程序中是否有与fluent相似的数据结构，即用单元的结构体链表存储各个单元的各向参数。我看到elevar（）函数可以调出的项里没有单元的生死状态。

还有，弱弱问一句，死单元里有温度和热流等参数吗？

代码如下：

subroutine uactive(m,n,mode,irststr,irststn,inc,time,timinc) 

   include 'D:\MSC.Software\Marc\2016.0.0\marc2016\common\implicit'

   dimension m(2)

   real,dimension(3)::coordin                      !定义坐标    

   integer,dimension(35354,2)::modee

c    

c   用来控制单元状态的子程序 

c 

c   输入参数: 

c   M(1) 单元号 

c   M(2) 在自适应分析中为母单元号 

c   N 积分点号 

c   INC 增量步号 

c   TIME 增量步开始的时间 

c   TIMINC 时间增量 

c    

c   输出参数: 

c   MODE: 单元的死活状态 

c  =-1,定义死单元 

c  =1,定义活单元 

c  =2,保持单元状态不变 

c   IRSTSTR:等于1时表示重新设置单元应力为0  

c   IRSTSTN:等于1时表示重新设置单元应变为0  

c* * * * * 

c 

cc  添加三个公共变量, 

cc  foreInc保存增量步 

cc  iopen文件打开标志 

cc  iend计算完成标志 

   common/check/foreInc,iopen,iend 

   data iopen,iend/0,0/ 

   real x0,y0,z0,x1,y1,z1,alpha,v,r,omega

c   (x0,y0,z0)为圆心位置,(x1,y1,z1)为热源中心位置,v为激光扫描线速度，r为扫描路径半径，omega为激光扫描角速度，alpha为参数方程参量

c   在0-19.5s阶段，即第一层熔覆过程，设置圆心位置为修复层上表面圆心，

c   在19.5-43.5s阶段，即第二层熔覆过程，设置圆心位置为修复层上表面圆心    

   integer::ie,mode,icode,kcus,nn !ie 为单元号，icode=0直接返回坐标，连续体单元kcus=1

   real::dist,dire !dist为遍历单元到热源中心距离，dire为数量积

   ie=m(1);

   if(ie.gt.25712.and.ie.lt.502035) return  !set直接跳过

   icode=0;kcus=1;nn=n

   call elmvar(icode,m,nn,kcus,coordin)                       

cc    

   if(time.le.19.5) then

   x0=175.54; y0=187.363;z0=-187.388

   else if (time.gt.19.5.and.time.le.43.5)then  

   x0=178.72;y0=190.471; z0=-185.311

   end if

c   设置激光扫描的线速度常数10,定义角速度omega

   v=10;omega=v/r;pi=3.1415926

c   设置激光扫描热源中心位置参数方程重要参量alpha

   if(time.le.13)then 

   r=19;alpha=omega*time

   else if (time.gt.13.and.time.le.19.5)then

   r=9;alpha=omega*time-13*omega

   else if(time.gt.19.5.and.time.le.35.5)then 

   r=24;alpha=omega*time-19.5*time

   else if (time.gt.35.5.and.time.le.43.5)then 

   r=12;alpha=omega*time-35.5*time

   end if

c   设置激光热源中心位置，假定热源在法线（1，1，1）的平面上运动

    x1=x0+(cos(alpha)-sin(alpha))*r/1.4142135

    y1=y0+cos(alpha)*r/1.4142135

    z1=z0+(cos(alpha)+sin(alpha))*r/1.73205081             热源位置     

cc

   dist=sqrt((coordin(1)-x1)**2+(coordin(2)-y1)**2

   $+(coordin(3)-z1)**2)                       积分点热源距离 

    dire=3.18*(175.54-coordin(1))+3.108*(187.363-coordin(2))

   $+2.077*(-187.388-coordin(3))                   数量积确定位置

    

c   初始激活第一层热源周围0.5mm的单元 ,死单元存入for96 

    if (inc.eq.0) then                          

     if(dist.le.0.5)then 

       if(dire.ge.0)then 

         mode=1

        else mode=-1

        end if

     else mode=-1

     end if

c  死单元初始化

   open(95,file='mode.txt',position='append')

   write(95,*)ie,mode

   close(95)

   go to 11

   end if

   open(95,file='mode.txt',status='old')

   do im=1,35454

     read(95,*) modee(im,1),modee(im,2)

   if(ie.eq.modee(im,1)) then

     mode=modee(im,2)

   end if 

   end do

   close (95)

      

     if(time+timinc.le.19.5)then 

       if(mode.eq.-1)then

         if(dist.le.11)then 

           if(dire.ge.0)then

             mode=1

           else mode=-1

           end if

         else mode=-1 

         end if

       else 

       return

        end if

        end if

   if(time+timinc.le.43.5.and.time+timinc.gt.19.5)then  

      if(mode.eq.-1)then

       if(dist.le.12)  then 

        mode=1

       else mode=-1;

       end if

            else mode=-1;

            end if

   else 

   return

        end if

c    一轮生死单元判断过后，更新mode.txt，并记录死单元

    open(95,file='mode.txt',status='old')

   close(95,status='delete')

c 删除原mode.txt文件，再重建更新

    open(95,file='mode.txt',status='new')

    do im=1,35454

      write(95,*) modee(im,1),modee(im,2)

            end do

    close(95)

    

c   

cc 

c   将当前增量步的死单元号写入文件 fort.96 中 

c   以便查看或者进行其他的处理 

    

  11 if(mode.eq.-1) then 

      write(96,*) 'deactivating element ', ie, ' increment ', inc 

   endif 

cc 

cc  新建一个过程文件 PostHide.proc 用来辅助后处理 

     if(iopen.eq.0) then 

      iopen=1 

      foreInc=inc 

      open(50,file='PostHide.proc',status='unknown') 

      close(50,status='delete') 

      open(50,file='PostHide.proc',status='new') 

      write(50,*)'*select_elements' 

   endif 

   if(iend.eq.0) then !判断是否到最后增量步 

    if(inc.eq.foreInc) then !判断增量步是否变化 

    if(mode.eq.-1) write(50,*) ie !写入当前增量步的死单元号 

    else 

    write(50,*)'#' !选择结束标志 

    write(50,*)'*invisible_selected' !使选中单元不可见 

    write(50,200) (inc-1) !显示第 (inc-1) 增量步 

 200  format('*post_skip_to',i5,/,'*select_elements') 

    write(50,*)'*animation_save' !保存当前窗口图形 

cc 

    if(mode.eq.-1) write(50,*) ie 

    foreInc=inc 

    endif 

    if(inc.eq.435) then !到最后一步,置 iend=1  

    iend=1 

    write(50,*)'#' 

    write(50,*)'*invisible_selected' 

    write(50,*)'*animation_save' 

    endif 

   endif 

cc  

c* * * * * 

   return 

   end 

代码尚未跑通，见笑了