1. A computer-implemented method of modeling thermal problems using a non-dimensional finite element method, comprising the steps of: (a) establishing a steady state heat conduction equation in three dimensions for an orthotropic material as k x ∂ 2 T ∂ x 2 + k y ∂ 2 T ∂ y 2 + k z ∂ 2 T ∂ z 2 + G = 0 , wherein the orthotropic material has three surfaces having surface areas S 1 , S 2 and S 3 , respectively, wherein T=T b holds for S 1 , k x ∂ T ∂ x l ^ + k y ∂ T ∂ y m ^ + k z ∂ T ∂ z n ^ + q * = 0 holds for S 2 , and k x ∂ T ∂ x l ^ + k y ∂ T ∂ y m ^ + k z ∂ T ∂ z n ^ + h ( T - T ∞ ) = 0 holds for S 3 , where k x , k y , and k z are thermal conductivities in Cartesian x, y and z-directions, respectively, G is a rate of internal heat generation per unit volume, T is temperature, T b and T ∞ are specified and ambient temperatures, respectively, q* is specified heat flux, h is a coefficient of convective heat transfer, and {circumflex over (l)}, {circumflex over (m)} and {circumflex over (n)} are surface normals, respectively corresponding to the Cartesian X, y and z-directions; (b) establishing a scalar functional quantity Π=U+Ω G +Ω q +Ω h , wherein U represents internal energy, Ω G represents internal heat generation, Ω h represents heat convection, and Ω q represents heat conduction, and: U = 1 2 ∫ ∫ ∫ V [ k x ( ∂ T ∂ x ) 2 + k y ( ∂ T ∂ y ) 2 + k z ( ∂ T ∂ z ) 2 ]  V , Ω G = - ∫ ∫ ∫ V GT  V , Ω q = - ∫ ∫ ∫ S 2 q * T  S , and Ω h = 1 2 ∫ ∫ S 3 h ( T - T ∞ ) 2  S , wherein V represents volume; (c) establishing a set of non-dimensional parameters θ = T T ∞ , x _ = x L x , y _ = y L y , and z _ = z L z , where θ is non-dimensional temperature, x , y and z are non-dimensional spatial variables, and L x , L y and L Z are maximum dimensions of a domain to be discretized along the x, y and z-directions, respectively; (d) expressing the scalar functional quantity for no internal heat generation in matrix form as: Π hT ∞ 2 = 1 2 ∫ ∫ ∫ V { g _ } T [ L _ ] { g _ }  V - ∫ ∫ S 2 { θ } T [ N _ ] T q _ *  S + 1 2 ∫ ∫ S 3 ( { θ } T [ N _ ] T - 1 ) 2  S , where [ N ] is a shape function matrix, { g } is a non-dimensional temperature gradient, {θ} is a non-dimensional temperature vector, [ D ] is a non-dimensional material property matrix, and [ L ] is a geometric dimensions ratio matrix, wherein { g _ } = [ ∂ θ ∂ x _ ∂ θ ∂ y _ ∂ θ ∂ z _ ] T and { g _ } = [ B _ ] { θ } , where [ B _ ] = [ ∂ [ N _ ] ∂ x _ ∂ [ N _ ] ∂ y _ ∂ [ N _ ] ∂ z _ ] T , and : [ D _ ] = [ 1 Bi x 0 0 0 1 Bi y 0 0 0 1 Bi z ] , [ L _ ] = [ 1 0 0 0 L x L y 0 0 0 L x L z ] , where Bi x , Bi y and Bi z are Biot numbers given by Bi x = hL x k x , Bi y = hL y k y and Bi z = hL z k z and non-dimensional specified heat flux q * is given by q _ * = q * hT ∞ ; (e) minimizing a dimensionless scalar functional quantity given as Π=Π/hT ∞ 2 L y L z with respect to {θ} to yield: [ ∫ ∫ ∫ V _ [ B _ ] T [ L _ ] [ D _ ] [ B _ ]  V _ + ∫ ∫ S _ 3 [ [ N _ ] T [ N _ ] ]  S _ ] { θ } = ∫ ∫ S _ 2 [ N _ ] T q _ *  S _ + ∫ ∫ S _ 3 [ N _ ] T  S _ ; (f) calculating an element stiffness matrix as [ k ]=[ k q ]+[ k h ], wherein [ k _ q ] = ∫ ∫ ∫ V _ [ B _ ] T [ L _ ] [ D _ ] [ B _ ]  V _ and [ k _ h ] = ∫ ∫ S 3 ⌊ [ N _ ] T [ N _ ] ⌋  S _ ; and (g) calculating an element load vector as { f }={ f q }+{ f h }, wherein: [ f _ q ] = ∫ ∫ S _ 2 [ N _ ] T q _ *  S _ and [ f _ h ] = ∫ ∫ S _ 3 [ N _ ]  S _ .