高斯激光热源下的n层材料瞬态热传导方程

xy112 · 发表于 2010-5-15 22:13:50

高斯激光热源下的n层材料瞬态热传导方程
高斯, 激光, 热源, 方程, 热传导高斯, 激光, 热源, 方程, 热传导
出于需要，写了个 n层材料的瞬态热传导方程，由表面激光加热。每层的热源由其光学性质N= n + ik决定. 通过解Maxwell方程获得激光在材料中的强度分布，由多项式 a0+a1*y+a2*y^2+a3*y^3+a4*y^4拟和获得。
P=Power/r^2/{PI}*exp(-({X}^2+({Z}-V*{TIME})^2)/r^2)*( A0+A1*(-{Y}/1e-9)+A2*(-{Y}/1e-9)^2+A3*(-{Y}/1e-9)^3+ A4*(-{Y}/1e-9)^4).

高斯分布体热源在程序中由HG.mac加载。不同层的a0,a1,a2,a3,a4决定了其内部的热源分布。

本程序可已根据问题选择所需的材料层数，所有参数都可以改变。

以下是命令流，如果有不合适的地方，欢迎大家提出、讨论！

FINISH  ! Make sure we are at BEGIN level
/Clear
/FILNAME,Media,0

/PREP7
/PNUM,VOLU,1

! In simulation, the X-Z is the media plane, while Y is the depth direction

Layer=5                      !Total layer number
nm= 1e-9                   !length unit
Power=0.005             ! Laser incident power W
r=5e-7                      ! Laser spot radius    m
V=0                         ! Laser moving veolcity, along z direction,  m/s

pw=Power*1000                !For display Title
rh=r/1e-6                         !For display Title

/TITLE,%Layer% layer media, Laser Power=%pw%mW, Heating radius=%rh%um, Velocity=%V%m/s

*DIM,Thickness,ARRAY,Layer ! Array input for each layer thickness
*DIM,k_L,ARRAY,Layer       ! Array input for each layer lateral thermal conductivity
*DIM,k_N,ARRAY,Layer       ! Array input for each layer normal thermal conductivity
*DIM,Cp,ARRAY,Layer          ! Array input for each layer heat capacity
*DIM,Den,ARRAY,Layer       ! Array input for each layer density
*DIM,div,ARRAY,Layer       ! Array input for each layer thickness division
*DIM,dscale,ARRAY,Layer    ! Array input for each layer thickness division scale
*DIM,A0,ARRAY,Layer          ! Array input for each layer a0 coefficient for heat generation curve
*DIM,A1,ARRAY,Layer          ! Array input for each layer a1 coefficient for heat generation curve
*DIM,A2,ARRAY,Layer          ! Array input for each layer a2 coefficient for heat generation curve
*DIM,A3,ARRAY,Layer          ! Array input for each layer a3 coefficient for heat generation curve
*DIM,A4,ARRAY,Layer          ! Array input for each layer a4 coefficient for heat generation curve

scale=10000                !for vovlap operation
Pulse_Duration=1e-8          !Laser pulse duration
End_Time=2e-8                !Final time for simulation
Step=20                      !Time interval
T0=50                      !Environment temperature

Length=4000*nm             !Length of the block in z direction
Width=2000*nm                !Width of the block in x direction

!Layer properties, from top to bottom

!First layer thermal property
Num=1
k_L(Num)=1.4
k_N(Num)=1.4
Den(Num)=4000
Cp(Num)=500
Thickness(Num)=50*nm
div(Num)=5
dscale(Num)=1                !last/first in -Y direction
A0(Num)=0
A1(Num)=0
A2(Num)=0
A3(Num)=0
A4(Num)=0

!Second layer thermal property
Num=2
k_L(Num)=16
k_N(Num)=16
Den(Num)=3000
Cp(Num)=900
Thickness(Num)=20*nm
div(Num)=2
dscale(Num)=1
A0(Num)=807.929E+6
A1(Num)=-43.896E+6
A2(Num)=934.529E+3
A3(Num)=-8.864E+3
A4(Num)=31.703

!Third layer thermal property
Num=3
k_L(Num)=1.5
k_N(Num)=1.5
Den(Num)=4000
Cp(Num)=500
Thickness(Num)=100*nm
div(Num)=5
dscale(Num)=1
A0(Num)=0
A1(Num)=0
A2(Num)=0
A3(Num)=0
A4(Num)=0

!Fourth layer thermal property
Num=4
k_L(Num)=100
k_N(Num)=100
Den(Num)=3900
Cp(Num)=700
Thickness(Num)=50*nm
div(Num)=2
dscale(Num)=2
A0(Num)=1.627E+9
A1(Num)=-32.146E+6
A2(Num)=238.195E+3
A3(Num)=-784.173
A4(Num)=967.716E-3

!Fifth layer thermal property
Num=5
k_L(Num)=0.2
k_N(Num)=0.2
Den(Num)=3900
Cp(Num)=450
Thickness(Num)=100*nm
div(Num)=5
dscale(Num)=2
A0(Num)=0
A1(Num)=0
A2(Num)=0
A3(Num)=0
A4(Num)=0

! Define the thermal element
ET,1,SOLID70
KEYOPT,1,2,0
KEYOPT,1,4,0
KEYOPT,1,7,0
KEYOPT,1,8,0

!For post-precessing
nPoint=200
display_r=1e-6

TSum1=0
TSum2=0
*do,Num,1,Layer,1
  UIMP,Num,DENS, , ,Den(Num),                                     !Define density of material Num
  UIMP,Num,KXX,KYY,KZZ,k_L(Num),k_N(Num),k_L(Num),                !Define thermal conductivity of material Num
  UIMP,Num,C, , ,Cp(Num),                                           !Define heat capacity of material Num
  TSum1=TSum1+Thickness(Num)
  BLOCK,-Width/2,Width/2,-TSum1,-TSum2,-Length/2,Length/2,          !Draw the layer of Num
  TSum2=TSum2+Thickness(Num)
*Enddo

d_total=TSum1                !Total thickness of the multilayer

VLSCALE,ALL, , , scale, scale, scale,  ,1,1,

*IF,Layer,GT,1,THEN
VGLUE,ALL
*ENDIF

VLSCALE,ALL, , , 1/scale, 1/scale, 1/scale, ,1,1,

TSum1=0
TSum2=0

Size=100*nm                            !Area mesh size, in X-Z plane
aesize,all,Size,

*do,Num,1,Layer,1
TSum1=TSum1+THICKNESS(Num)
lsel,all

!Determine the mesh scale and divison for each layer at Y direction
*do,i,-1,1,2
   *do, j,-1,1,2
      lsel,s,loc,y,-TSum1,-TSum2
      lsel,r,loc,x,i*Width/2,i*Width/2
      lsel,r,loc,z,j*Length/2,j*Length/2
      *GET,Line_Num,LINE, 0, NUM,MIN,
      *GET,KP1,LINE,Line_Num,KP,1,
      *GET,KP2,LINE,Line_Num,KP,2,
      *GET,KP1_Y,KP,KP1, LOC, Y,
      *GET,KP2_Y,KP,KP2, LOC, Y,
      *IF,KP1_Y,GT,KP2_Y,THEN
         lesize,all,THICKNESS(Num)/div(Num), , ,dscale(Num),1, , ,1,
      *ELSE
         lesize,all,THICKNESS(Num)/div(Num), , ,1/dscale(Num),1, , ,1,
      *ENDIF
   *ENDDO
*ENDDO

asel,all

esize,Size,0,

VSEL,S,LOC,Y,-TSum1,-TSum2

TYPE,1
MAT,Num
mshkey,1
mshape,0,3d
vmesh,ALL          !Mesh the layer Num

allsel,all

TSum2=TSum2+Thickness(Num)

*enddo

/REP,FAST
/psf,hflu,,2,0
/PNUM,LINE,1
finish
/solu

antype,trans          !Transient simulation
tunif,T0                !Ambient temperature
TREF,T0
timint,on !time integration = on
kbc,1
OUTRES,ALL,ALL,       !Output all the result

/view, 1 ,1,1,1
/ang, 1
/rep,fast

*do,Num,1,Layer,1
*USE,HG.mac,'Layer%Num%',Power,r,V,A0(Num),A1(Num),A2(Num),A3(Num),A4(Num)    !Load heat generation for each layer
*enddo

!Start simulation
*do,tm,1e-15,End_Time,End_Time/Step

time,tm
autots,0
deltim, , , ,1
kbc,0

allsel,all
TSum1=0
TSum2=0
Num=0

*if,tm,LE,Pulse_Duration+1e-12,then

   !Apply the heat generation in each layer
   Num=1
   TSum1=TSum1+THICKNESS(Num)
   nsel,s,LOC,y,-TSum1,-TSum2
   BF,all,HGEN,%Layer1%
   TSum2=TSum2+Thickness(Num)

   Num=2
   TSum1=TSum1+THICKNESS(Num)
   nsel,s,LOC,y,-TSum1,-TSum2
   BF,all,HGEN,%Layer2%
   TSum2=TSum2+Thickness(Num)

   Num=3
   TSum1=TSum1+THICKNESS(Num)
   nsel,s,LOC,y,-TSum1,-TSum2
   BF,all,HGEN,%Layer3%
   TSum2=TSum2+Thickness(Num)

   Num=4
   TSum1=TSum1+THICKNESS(Num)
   nsel,s,LOC,y,-TSum1,-TSum2
   BF,all,HGEN,%Layer4%
   TSum2=TSum2+Thickness(Num)

   Num=5
   TSum1=TSum1+THICKNESS(Num)
   nsel,s,LOC,y,-TSum1,-TSum2
   BF,all,HGEN,%Layer5%
   TSum2=TSum2+Thickness(Num)

*else
   allsel,all
   BFDELE, ALL,HGEN
*endif

allsel,all
asel,s,LOC,y,-1e-9,0
! DA,ALL,TEMP,T0
allsel,all

/status,solu
solve
*enddo

finish
/POST1

!*ULIB, Post,mac
!*USE,TTV !Temperature decay at each layer thickness
!*USE,TRV,1  !Normal Temperature profile at each layer at given time
!*USE,TRL,1  !Plane Temperature profile at each layer at given time

!Laser power absorption in each layer
!Power/r^2/{PI}*exp(-({X}^2+({Z}-V*{TIME})^2)/r^2)*( A0+A1*(-{Y}/1e-9)+A2*(-{Y}/1e-9)^2+A3*(-{Y}/1e-9)^3+ A4*(-{Y}/1e-9)^4)

//UG.mac

!FileName=ARG1
!Power=ARG2
!r=ARG3
!V=ARG4
!A0=ARG5
!A1=ARG6
!A2=ARG7
!A3=ARG8
!A4=ARG9

*DEL,_FNCNAME
*DEL,_FNCMTID
*DEL,_FNC_C1
*DEL,_FNC_C2
*DEL,_FNC_C3
*DEL,_FNC_C4
*DEL,_FNC_C5
*DEL,_FNC_C6
*DEL,_FNC_C7
*DEL,_FNC_C8
*SET,_FNCNAME,ARG1
*DIM,_FNC_C1,,1
*DIM,_FNC_C2,,1
*DIM,_FNC_C3,,1
*DIM,_FNC_C4,,1
*DIM,_FNC_C5,,1
*DIM,_FNC_C6,,1
*DIM,_FNC_C7,,1
*DIM,_FNC_C8,,1
*SET,_FNC_C1(1),ARG2
*SET,_FNC_C2(1),ARG3
*SET,_FNC_C3(1),ARG4
*SET,_FNC_C4(1),ARG5
*SET,_FNC_C5(1),ARG6
*SET,_FNC_C6(1),ARG7
*SET,_FNC_C7(1),ARG8
*SET,_FNC_C8(1),ARG9
! /INPUT,HG.func
*DIM,%_FNCNAME%,TABLE,10,61,1
!
! Begin of equation: Power/r^2/{PI}*exp(-({X}^2+({Z}-V*{TIME})^2)/r^2)*(A0+A1*
! (-{Y}/1e-9)+A2*(-{Y}/1e-9)^2+A3*(-{Y}/1e-9)^3+ A4*(-{Y}/1e-9)^4)
*SET,%_FNCNAME%(0,0,1), 0.0, -999
*SET,%_FNCNAME%(2,0,1), 0.0
*SET,%_FNCNAME%(3,0,1), %_FNC_C1(1)%
*SET,%_FNCNAME%(4,0,1), %_FNC_C2(1)%
*SET,%_FNCNAME%(5,0,1), %_FNC_C3(1)%
*SET,%_FNCNAME%(6,0,1), %_FNC_C4(1)%
*SET,%_FNCNAME%(7,0,1), %_FNC_C5(1)%
*SET,%_FNCNAME%(8,0,1), %_FNC_C6(1)%
*SET,%_FNCNAME%(9,0,1), %_FNC_C7(1)%
*SET,%_FNCNAME%(10,0,1), %_FNC_C8(1)%
*SET,%_FNCNAME%(0,1,1), 1.0, -1, 0, 2, 0, 0, 18
*SET,%_FNCNAME%(0,2,1), 0.0, -2, 0, 1, 18, 17, -1
*SET,%_FNCNAME%(0,3,1), 0, -1, 0, 1, 17, 4, -2
*SET,%_FNCNAME%(0,4,1), 0.0, -2, 0, 3.14159265358979310, 0, 0, -1
*SET,%_FNCNAME%(0,5,1), 0.0, -3, 0, 1, -1, 4, -2
*SET,%_FNCNAME%(0,6,1), 0.0, -1, 0, 0, 0, 0, 0
*SET,%_FNCNAME%(0,7,1), 0.0, -2, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,8,1), 0.0, -4, 0, 1, -1, 2, -2
*SET,%_FNCNAME%(0,9,1), 0.0, -1, 0, 2, 0, 0, 2
*SET,%_FNCNAME%(0,10,1), 0.0, -2, 0, 1, 2, 17, -1
*SET,%_FNCNAME%(0,11,1), 0.0, -1, 0, 1, 19, 3, 1
*SET,%_FNCNAME%(0,12,1), 0.0, -5, 0, 1, 4, 2, -1
*SET,%_FNCNAME%(0,13,1), 0.0, -1, 0, 2, 0, 0, -5
*SET,%_FNCNAME%(0,14,1), 0.0, -6, 0, 1, -5, 17, -1
*SET,%_FNCNAME%(0,15,1), 0.0, -1, 0, 1, -2, 1, -6
*SET,%_FNCNAME%(0,16,1), 0.0, -2, 0, 1, -4, 3, -1
*SET,%_FNCNAME%(0,17,1), 0.0, -1, 0, 2, 0, 0, 18
*SET,%_FNCNAME%(0,18,1), 0.0, -4, 0, 1, 18, 17, -1
*SET,%_FNCNAME%(0,19,1), 0.0, -1, 0, 1, -2, 4, -4
*SET,%_FNCNAME%(0,20,1), 0.0, -1, 7, 1, -1, 0, 0
*SET,%_FNCNAME%(0,21,1), 0.0, -2, 0, 1, -3, 3, -1
*SET,%_FNCNAME%(0,22,1), 0.0, -1, 0, 0, 0, 0, 0
*SET,%_FNCNAME%(0,23,1), 0.0, -3, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,24,1), 0.0, -4, 0, 1, -1, 2, -3
*SET,%_FNCNAME%(0,25,1), 0.0, -1, 0, 1, -4, 3, 3
*SET,%_FNCNAME%(0,26,1), 0.0, -3, 0, 1e-9, 0, 0, -1
*SET,%_FNCNAME%(0,27,1), 0.0, -4, 0, 1, -1, 4, -3
*SET,%_FNCNAME%(0,28,1), 0.0, -1, 0, 1, 21, 3, -4
*SET,%_FNCNAME%(0,29,1), 0.0, -3, 0, 1, 20, 1, -1
*SET,%_FNCNAME%(0,30,1), 0.0, -1, 0, 0, 0, 0, 0
*SET,%_FNCNAME%(0,31,1), 0.0, -4, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,32,1), 0.0, -5, 0, 1, -1, 2, -4
*SET,%_FNCNAME%(0,33,1), 0.0, -1, 0, 1, -5, 3, 3
*SET,%_FNCNAME%(0,34,1), 0.0, -4, 0, 1e-9, 0, 0, -1
*SET,%_FNCNAME%(0,35,1), 0.0, -5, 0, 1, -1, 4, -4
*SET,%_FNCNAME%(0,36,1), 0.0, -1, 0, 2, 0, 0, -5
*SET,%_FNCNAME%(0,37,1), 0.0, -4, 0, 1, -5, 17, -1
*SET,%_FNCNAME%(0,38,1), 0.0, -1, 0, 1, 22, 3, -4
*SET,%_FNCNAME%(0,39,1), 0.0, -4, 0, 1, -3, 1, -1
*SET,%_FNCNAME%(0,40,1), 0.0, -1, 0, 0, 0, 0, 0
*SET,%_FNCNAME%(0,41,1), 0.0, -3, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,42,1), 0.0, -5, 0, 1, -1, 2, -3
*SET,%_FNCNAME%(0,43,1), 0.0, -1, 0, 1, -5, 3, 3
*SET,%_FNCNAME%(0,44,1), 0.0, -3, 0, 1e-9, 0, 0, -1
*SET,%_FNCNAME%(0,45,1), 0.0, -5, 0, 1, -1, 4, -3
*SET,%_FNCNAME%(0,46,1), 0.0, -1, 0, 3, 0, 0, -5
*SET,%_FNCNAME%(0,47,1), 0.0, -3, 0, 1, -5, 17, -1
*SET,%_FNCNAME%(0,48,1), 0.0, -1, 0, 1, 23, 3, -3
*SET,%_FNCNAME%(0,49,1), 0.0, -3, 0, 1, -4, 1, -1
*SET,%_FNCNAME%(0,50,1), 0.0, -1, 0, 0, 0, 0, 0
*SET,%_FNCNAME%(0,51,1), 0.0, -4, 0, 1, 0, 0, -1
*SET,%_FNCNAME%(0,52,1), 0.0, -5, 0, 1, -1, 2, -4
*SET,%_FNCNAME%(0,53,1), 0.0, -1, 0, 1, -5, 3, 3
*SET,%_FNCNAME%(0,54,1), 0.0, -4, 0, 1e-9, 0, 0, -1
*SET,%_FNCNAME%(0,55,1), 0.0, -5, 0, 1, -1, 4, -4
*SET,%_FNCNAME%(0,56,1), 0.0, -1, 0, 4, 0, 0, -5
*SET,%_FNCNAME%(0,57,1), 0.0, -4, 0, 1, -5, 17, -1
*SET,%_FNCNAME%(0,58,1), 0.0, -1, 0, 1, 24, 3, -4
*SET,%_FNCNAME%(0,59,1), 0.0, -4, 0, 1, -3, 1, -1
*SET,%_FNCNAME%(0,60,1), 0.0, -1, 0, 1, -2, 3, -4
*SET,%_FNCNAME%(0,61,1), 0.0, 99, 0, 1, -1, 0, 0
! End of equation: Power/r^2/{PI}*exp(-({X}^2+({Z}-V*{TIME})^2)/r^2)*(A0+A1*(-
! {Y}/1e-9)+A2*(-{Y}/1e-9)^2+A3*(-{Y}/1e-9)^3+ A4*(-{Y}/1e-9)^4)
!-->

		自动登录	找回密码
密码			注册