Working wonders with APDL Math – Ep 02 : Thermal Harmonic Analysis

 Thermal Harmonic Analysis using APDL Math

To help developing an intuition, let me review the similarities and differences between structural and thermal frequency response analysis.

Mathematically, both look very much the same, i.e. what we are seeking are solutions to dynamic equilibrium equations:

Domain	Equation solved	Terms explained
Structural		stiffness damping mass force vector
Thermal		conductivity capacity power vector

Both analyses will yield complex solution vectors, giving real (in phase) and imaginary (out of phase) part of the structural response, for a harmonic excitation at circular frequency.

Now, the fundamental difference is the following: while inertia forces are 180° out of phase with elastic forces, there will be frequencies where cancellation occurs and hence the structure will have zero apparent stiffness, that is, resonance will occur, only limited by damping forces.

In the thermal case, the capacitance terms will remain 90° out of phase with the conductivity terms, hence cancellation will never occur. Consequently, dynamic temperature response amplitudes will always have amplitude smaller than that of the static equilibrium case, and phase will always be negative.

Needless to say, if the loading consists of enforced temperature at the system boundary, the response temperature field should be uniform for quasi-static response (i.e. for temperature variations with periods much longer than the first thermal mode time constant, see first post on the subject).

Now let’s turn to the APDL input deck.

1. It reads in conductivity, capacitances matrices, and one forcing (RHS) vector,
2. It repeatedly solves for the entire problem at frequencies ranging from F_begin to F_end.
3. Then, the response at user-selected key nodes are stored into two matrices(TsolR and TsolI), and written to a text file
4. Optionally, the entire solution field can be stored into two .RTH files (one for real part, the other for imaginary part). This complex solution can be either displayed using POST 1 capabilities or used as inputs for a downstream structural analysis (using LDREAD).

Mind you, we are facing a bit of bookkeeping, because at a later point we might want to switch to a structural analysis. It is therefore my recommendation to change job names for each step of the analysis (generation of conductivity and capacitance matrices, load vector, thermal FRF, etc), following the same prophylaxis rules than for superelements.

Note that if your problem involves prescribed temperatures instead of prescribed nodal power, the script need not to be changed, the only requirement being that prescribed DOF values (D command) are applied during as soon as the matrix generation step (actual DOF values are not used at that step, but prescribed DOFs values are eliminated from the solution phase).

To emphasize, unless prescribed temperature node set are strictly identical, load cases won’t solve using the same solver engine (i.e. they need to be written to separate .FULL files), because the set of unknown degrees of freedom won’t be the same (just like the a-set and u-set of degree of freedom, for those of you who are familiar with Nastran DMAP. And yes, I happened to use that tool too, back in the days. Go on, hiss!). Have a look at Chapter 14.6 of the Mechanical APDL “Theory Reference” Manual if you need details.

! Setup Model
GenerationJobName=’Gene_Pass’
/filnam,%GenerationJobName%,ON

…
! Make THREE dummy transient solve to separately write conductivity, capacitances and load vector matrices to .FULL files named Run_Thermal_MatK, Run_Thermal_MatC, Run_Thermal_VecF 

fini

Also, it is good practice to store into a variable the jobname used at the generation pass, so that will be saved in the database. This way, you might interoperate scripts much more easily.

Next, we need to specify solver options. In this case, only the frequency bandwidth needs to be specified. For thermal problems, I suggest sticking to the regular spacing (in a logarithmic sense). As engineering unit, I personally use cycles per day (cpd) rather than cycles per second (Hz), but that is a matter of personal taste and of course of the physics you need to address. Secondly, we want to select a few nodes at which response temperature is of particular significance for the problem at hand.

! Change jobname to avoid contaminating .FULL or results file
/filnam,Run_FRF_Thermal,ON

! Load parameters
F_begin=.1			! lowest frequency [cpd]
F_end=1440			! highest frequency [cpd]
F_Nb_Steps=30		! Number of steps to be used for freq. sweep [1]

! Output Parameters
NbOutputNodes=3
ListOutputNodes(1)=MyPreferredNode,AnotherNode,OneLastNode

! STEP1: IMPORT REQUIRED DATA
! Conductivity and capacitance matrices
*SMAT,MatK,D,IMPORT,MAT,Run_Thermal_MatK
*SMAT,MatC,D,IMPORT,MAT,Run_Thermal_MatC

! RHS vector
*DMAT,VecF_real,D,IMPORT,MAT,Run_Thermal_VecF	! VecF is real
*DMAT,VecF,Z,COPY,VecF_real		! but it must be used as complex valued

! Input mapping matrice/ renumbering vector
*SMAT,Nod2Bcs,D,IMPORT,FULL,strcat(GenerationJobname,'.full'),NOD2BCS
*VEC,MapForward,I,IMPORT,FULL,strcat(GenerationJobname,'.full'),FORWARD

! Allocation for the solution vector (XBCs), same size as forcing vector
*DMAT,XBcs,Z,COPY,VecF 

! Allocate matrices to separately store real and imaginary part of solution for nodes at keypoints
*dim,TsolR,array,F_nb_steps,NbOutputNodes
*dim,TsolI,array,F_nb_steps,NbOutputNodes

! Allocate matrices to separately store real and imaginary part of solution for ALL nodes on the model 
! Prepare component containing nodes without prescribed Boundary conditions, this will set the length of the solution vector
cmsel,s,MyNodeComponent
*get,Nb_Nd_Res,NODE,0,COUNT

*dim,TsolMatR,array,F_nb_steps,Nb_Nd_Res
*dim,TsolMatI,array,F_nb_steps,Nb_Nd_Res

Now we are all set and we can begin solving the problem. For each frequency, we build a complex valued matrix (MatA) that needs to be solved. As usual with large numerical problems, we do not invert the matrix but rather we factorize it (using *LSENGINE/*LSFACTOR) and then use back substitution (using *LSBAC).

! STEP3: SOLVE
! Create excitation frequency table: even spacing in logarithmic sense
*dim,F_table,array,F_nb_steps
*vfill,F_table,RAMP,log(F_begin),(log(F_end)-log(F_begin))/(F_nb_steps-1)
*VFUN, F_table, EXP, F_table

*DO,ind_F,1,F_nb_steps			! LOOP OVER FREQUENCY VALUES
  curr_F=F_table(ind_F)			! frequency [cpd]
  w=(2*3.14159)*curr_F/86400		! pulsation [rad/s]

  ! FORM THE COMPLEX SYSTEM  -> A =K+jw*C
  *SMAT,MatA,Z,COPY,MatC	! MatA=MatC
  *AXPY,1.,0.,MatK,0.,w,MatA	! MatA=(1xMatK)+jw*MatA

  ! FACTORIZE MATRIX A
  *LSENGINE,BCS,MyBcsSolver,MatA
  *LSFACTOR,MyBcsSolver

  ! SOLVE THE LINEAR SYSTEM, PUT RESULT INTO VECTOR TBCS
  *LSBAC,MyBcsSolver,VecF,TBcs  
  
  ! SEPARATE REAL AND IMAGINARY PART OF SOLUTION VECTOR
  *vec,TBcsR,D,COPY,TBcs,REAL
  *vec,TBcsI,D,COPY,TBcs,IMAG
	
  ! TRANSFER TO USER ORDERING
  *MULT,Nod2Bcs,TRAN,TBcsR,,TintR
  *MULT,Nod2Bcs,TRAN,TBcsI,,TintI
  
  ! Store complete solution
  cmsel,s,MyNodeComponent
  curr_node=0
  *do,ind_node,1,nb_nd_res
	curr_node=ndnext(curr_node)
	TsolMatR(ind_F,ind_node)=TintR(MapForward(curr_node))
	TsolMatI(ind_F,ind_node)=TintI(MapForward(curr_node))
  *enddo
  
  ! Store solution at key locations
  *do,ind_node,1,NbOutputNodes
    Local_index=MapForward(ListOutputNodes(ind_node))
    TsolR(ind_F,ind_node)=TintR(Local_index)
    TsolI(ind_F,ind_node)=TintI(Local_index)
  *enddo
*ENDDO

Note that for the sake of readability, the code is extremely crude, and performance is limited by loop need to store the complete solution. This step becomes significant for large models (>10⁵DOFs) and could be speeded up by usage of the *VGET / *VPUT commands (and guess what, this was already covered here: http://www.padtinc.com/blog/the-focus/changing-results-values-in-ansys-mechanical-apdl-dnsol-and-vput)

Lastly, here is a code snippet to arrange the complex temperature vectors into a single matrix that is suitable to use with *VWRITE. What is special about it is that it automatically adjusts the format field to the number of output nodes, using the old command file trick. It has nothing particularly difficult to write, but still it took me a few attempts to achieve. I don’t know about you, but there must be something about *VWRITE formats, and “I keep forgetting where it is” (Depardieu sait où, for those of you that are French connoisseurs). Going back to *VWRITE command and formats, probably I should consider buying another copy of PADT’s “Introduction to the ANSYS Parametric Design Language” so that I have one spare.

Enough nonsense, here is the code:

! CREATE A GLOBAL MATRIX CONTAINING REAL AND IMAGINARY PARTS OF SOLUTION
*dim,TsolC,array,F_nb_steps,1+2*NbOutputNodes

*do,ind_f,1,F_nb_steps
	TsolC(ind_f,1)=F_table(ind_f)
*enddo

*do,ind_node,1,NbOutputNodes
	ind_r=2+2*(ind_node-1)
	ind_i=ind_r+1
	*do,ind_f,1,F_nb_steps
		TsolC(ind_f,ind_r)=TsolR(ind_f,ind_node)
		TsolC(ind_f,ind_i)=TsolI(ind_f,ind_node)
	*enddo
*enddo

*cfopen,command_file,inp
*vwrite,ExcType,NbOutputNodes
*mwrite,TsolC(1,1),Run_FRF_ExcType=%c_NbNd=%i,txt
*vwrite,'(',1+2*NbOutputNodes
%c F8.2,X,%i(F12.8))
*cfclose

/inpu,command_file,inp
/sys,del command_file,inp