User login


You are here

Force Matrix for Transient Thermal Analysis


I have developed a 2D transient thermal analysis FE code which doesn't
seem to be correct. I am sure there is a conceptual problem in the
formulation itself which I have not been able to understand and could
not find concrete information in books.

I am putting up the question in a very simplified form so as to address the correctness of the formulation.

Let's say we have a square of unit length with 4 nodes at each corner
(Q4 element). Initial temperatures are specified as 100 for all nodes
(at t=0). The temperature of the nodes 1 and 2 are now specified as 200
and 300 and we investigate the temperature distribution on other 2
nodes with time say 0 to 100 seconds (with the specified temperatures
maintained). No other thermal load is applied.

-- Now, while solving the FE discreet equation, during the first time step would there be any force matrix on the R.H.S?

-- If no, then would the R.H.S be simply (1/dt)*C*tn (using 'theta=1').
'C' being the capacitance matrix, 'tn' is vector of initial dofs and
'dt' is time-step. (Reference: Ansys theory manual, transient thermal
analysis, eq 17-33)

I have tried to keep the problem as simple as possible without much
theory. Hope it is clear. Any suggestion or comment is deeply

Research Assistant
Structural and Multidisciplinary Optimization Lab.
University of Florida

Hi Shrimad,
I couldn't check what you are really asking. Could you please explain more accurately what the problem is?

"...with 4 nodes at each corner..." do you mean with one node at each corner? Seems you are talking about just one element?! strange! 
"...'C' being the capacitance matrix,..." I guess you mean the specific heat matrix, is it so?

As you said it is a "Transient Thermal Analysis", then what do you mean by theta=1? I mean the temperature field can be explained with the heat equation which involves geometry (x ,y , z), density  and specific heat. Now what is theta here.

I don't have access to Ansys help, certainly they have solved the heat equation numerically, finally sth. like m.T_dot+K.T=F ;T is temperature and T_dot the time rate of temp. change. Without the first term it would be a steady state problem. K has the same form as a stiffness matrix (B.C.B_transpose).

Just one more point please,  "... The temperature of the nodes 1 and 2 are now specified as 200 and 300..." Does now mean t=0?



Hi. Thanks for the reply. I am explaining the issues that you had mentioned.


-- Since I am just testing to check the correctnes of the formulation, I am using just 1 element i.e. a unit square with 4 nodes (Q4) element. It becomes easy to check the matrices etc.


-- C is the specific heat matrix


--  'theta' is the time integration parameter which lies between [0,1]. It's analogous to the Newmark's parameter for time integration. It's usually 1/2 and in that case it is called the Crank - Nicolson rule. I am using a value of 1. Please see eqtn 17-33 on page 894 at this link .This is precisely the equation I am solving.


--  The initial condition, i.e. temperature at t=0 is 100 at all nodes. At the first time step, the temperature on nodes 1 and 2 is specified as 200 and 300. We investigate the temperature change on other 2 nodes with time. (The specified temperature is maintained throughout).


 My question is: Will there be an additional force matrix in the RHS due to the specified temperature or it is simple the RHS of eq 17-33? How are the reaction loads calculated? 


Thanks for your interest!




Hi Shriran,
OK,I see! With theta =1  the equation 17.33 in your case changes to an Euler backward difference, for the approximation of Temperature derivative in time
(dT/dt at time n+1 also T_dot see above).
So set theta =1 and simplify you get:
C.(T_n+1 - T_n)/(delta t ) + K. T_n+1 = F_a.

Will there be an additional force matrix in the RHS?
No,  f_a is the vector of applied heat flows. But in Ansys help the discussion is more about the time approximation method as said in your case would be a backward Euler difference and unconditionally stable.
You are actually trying to solve the transient heat equation. If you use FEM, for the solution of the BVP you use let's say a linear approximating function N.  So you get N , B , C , K ... if you have the source term of the heat equation in your FEM Formulation then F_ a is actually the sum of two kinds of possible temperature load vectors or BC's.
I think you don't define a transient problem with the above conditions. Besides  T_0  which is temperature at t0 (initial condition 100 in your case) impose two  boundary conditions for example known temperature along a boundary say 200 and and no heat flow across another boundary. You may want to use two Q4-elements the global matrix would be a 6x6 matrix which is still not so larger than 4x4 case. Still you shouldn't expect a good result because of the course grid.

Hope this would help.


Hi Roozbeh,

Thanks for the reply. I still could not get the right idea about the force matrix. I want to give you a very simple example because I don't get the expected results which are obvious. Here I am not at all concerned with the mesh density as this is just a test problem. So I am just using 1 element with 4 nodes.Small 't' represents time and caps 'T' is temperature.


-- As I said before, the initial temperature  (t0) on all the nodes is 100 i.e. Tn = [100 100 100 100]. To further simplify let us assume that the conductivity matrix i.e. K = 0.


-- Now at t1 the temperature on node 1 is raised to 200. I want to find out the temperatures on other three nodes at t1. Obviously the temperatures on other 3 nodes would still be [100 100 100], as the conductivity matrix is 0. So, T1 should be [200 100 100 100]. This is what I get from ANSYS too and there is a reaction load or heat generated at node 1 and 0 on others (from ANSYS). This is also what we expect.


-- So if we try to solve this in the equation that you wrote, we will have C*[200 T2 T3 T4]' = C*[100 100 100 100]' + [R 0 0 0]. Where R is the reaction load at node 1 or the heat generated. Solving this leads to [T2 T3 T4] = [50 125 50] and [R=4692] which is not correct.


-- We don't need to concentrate on the values or the numbers. The above method seems perfectly fine to me but I get a ridiculously wrong result. The nodal solutions must be T(n+1) = [200 100 100 100] and not [200 50 125 50]. Do you see anything wrong in the above method of solution?


Thanks for your interest!


Hi Shriram,

With different C's you would get different solutions. May be the C matrix which you use is not the same as the one in Ansys. The formula above is about the time approximation and the the form of C , K and F_a are determined in your FEM formulation of the BVP. The C matrix has the form ∫∫s*N*N_transpose*dA, in which s would be the coefficient of the time derivative of Temperature in the heat equation, provided we use the same heat eq.. 
I am sure I have seen this in a FEM book . I will tell you if I find.




Yes you are right, I calculate the 'C' matrix by  ∫∫s*N*N_transpose*dA. Here 's' is the product of 'ρ' and 'c' where 'ρ' is the density and 'c' is the specific heat, both of them being constants.

And I define the C matrix exactly in ANSYS that I use by defining 'ρ' and 'c' . So the 'C' matrix cannot be different.



Research Assistant
Structural and Multidisciplinary Optimization Lab.
University of Florida

Hi shrimad,

The force vector is not zero at the node you have specified its temperature variation versus time. I think there is a conceptual mistake in the way you use the matrix equation from FEM model. As you know, we could not determine both temperature and heat flux on the boundary nodes. So if we determine temperature on a node of the 4-node element, then we should let the model to determine the heat flux on that node to produce the demanded temperature variation. The force vector on that node is nonzero and also has an unknown value which can be obtained from FEM model. But we should not take care of that unknown value in the force vector when solving system of equations (similar to the case that it is zero), because after imposing primary boundary conditions (like temperature) on the nodes and in the stage of solving the system of equations, we automatically neglect those equations which have predetermined values at the left hand side vector and unknown values at the right hand side vector. After solving system of equations and obtaining temperatures, then the unknown values in the force vector could be obtained by the equations firstly neglected.

Hope this helps.


Subscribe to Comments for "Force Matrix for Transient Thermal Analysis"

Recent comments

More comments


Subscribe to Syndicate