So the strong form of the heat diffusion and convection PDE is given as
$\rho c_m \mathbf{v}\cdot\nabla T – \nabla \cdot \nabla T = \dot{q}\\
T(\mathbf{x},t) = T_e(\mathbf{x},t)~~~ on ~~~\Gamma_e ~~~~(\text{Dirichlet-BC})\\
k \frac{\partial T}{\partial \mathbf{n}} = q_n~~~ on ~~~\Gamma_n ~~~~(\text{Neumann-BC})\\
\mathbf{v}….\text{velocity, given}$
Then I introduced a test function $w$ and derived the weak form of the PDE ($\Omega$ – volume, $\Gamma$ – surface):
$
\int_\Omega w \rho c_m \mathbf{v} \cdot\nabla T~ d\Omega + \int_\Omega (\nabla w)\cdot (k \nabla T) ~d\Omega = \int_\Omega w \dot{q} ~d\Omega + \int_{\Gamma_n} w q_n~d\Gamma
$
From that I want to derive Galerkin's finite element formulation by approximating the domain $\Omega$ and the function spaces $T(\mathbf{x},t)$ and $w(\mathbf{x})$ by known shape functions $N_i(\mathbf{x})$.
$
\int_\Omega (\nabla w)\cdot (k \nabla T)~ d\Omega \approx \sum_{m=1}^M \int_{\Omega^{(m)}} (\nabla w)\cdot (k \nabla T)~ d\Omega \\
w(\mathbf{x}) \approx \sum_{i=1}^{N^{(m)}} w_i N_i(\mathbf{x})\\
~~~~~~~~~~~~~~N_i(\mathbf{x_j}) = \delta_{ij}
$
When inserting this approximations in the weak form of the PDE, I can obtain a single equation for node i=k by assuming that $w_k = 1$ at the node and $w_k=0$ otherwise. I get the final result as:
$
\sum_{m=1}^{M_e} \sum_{m=1}^{N^{(m)}} [T_j \int_{\Omega^{(m)}} (N_i(\mathbf{x})\rho c_m \mathbf{v} \cdot \nabla N_j(\mathbf{x})+k(\nabla N_i(\mathbf{x})\cdot \nabla N_j(\mathbf{x}))) d\Omega]= \sum_{m=1}^{M_e} \int_{\Omega^{(m)}} N_i(\mathbf{x}) \dot{q} d\Omega + \sum_{m=1}^{M_n} \int_{\Gamma_n^{(m)}} N_i(\mathbf{x}) q_n d\Gamma
$
which I can write in Matrix-Vector form as the following:
$\sum_{j=1}^{N^{(m)}} K_{ij}^{(m)}T_j^{(m)} = f_i^{(m)}, ~~~~~~j = 1,…,N^{(m)} $
with the system matrix
$K_{ij}^{(m)} = \int_{\Omega^{(m)}} (N_i(\mathbf{x})\rho c_m \mathbf{v} \cdot \nabla N_j(\mathbf{x})+k(\nabla N_i(\mathbf{x})\cdot \nabla N_j(\mathbf{x}))) d\Omega$
But now my question: can this be right? Because I read that the system matrix normally should be symmetric and this is not fullfilled here?
Thanks for any help!
Best Answer
The calculation seems fine to me. I have not gone in details but from an overlook, it looks fine. Now to answer your main question, why isn't the matrix symmetric? The answer is the matrix can be non-symmetric! The Symmetricity of the matrix depends on the symmetric nature of the bilinear form. In fact, the bilinear form is symmetric if and only if the matrix is symmetric. Mostly standard books on FEM considers the Poisson problem, there we have the symmetric nature of the bilinear form and hence we get a symmetric matrix. The problem you have doesn't have a symmetric bilinear form