In FEM classes, it's usually taken for granted that the stiffness matrix is positive definite, but I just can't understand why. Could anyone give some explanation?
For instance, we can consider the Poisson problem: $$ -\nabla^2 u = f,$$ whose stiffness matrix is: $$K_{ij} = \int_\Omega\nabla\varphi_i\cdot\nabla\varphi_j\, d\Omega,$$ which is symmetric and positive definite. Symmetry is an obvious property, but the positive definiteness is not so explicite to me.
The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods.
To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary conditions here): Find $u\in H^1_0(\Omega)$ such that $$ a(u,v):= \int_\Omega \nabla u\cdot \nabla v \,dx = \int_\Omega fv\,dx \qquad\text{for all }v\in H^1_0(\Omega).$$ The important property here is that $$ a(v,v) = \|\nabla v\|_{L^2}^2 \geq c \|v\|_{H^1}^2 \qquad\text{for all }v\in H^1_0(\Omega). \tag{1}$$ (This follows from Poincaré's inequality.)
Now the classical finite element approach is to replace the infinite-dimensional space $H^1_0(\Omega)$ by a finite-dimensional subspace $V_h\subset H^1_0(\Omega)$ and find $u_h\in V_h$ such that $$ a(u_h,v_h):= \int_\Omega \nabla u_h\cdot \nabla v_h \,dx = \int_\Omega fv_h\,dx \qquad\text{for all }v_h\in V_h.\tag{2}$$ The important property here is that you are using the same $a$ and a subspace $V_h\subset H^1_0(\Omega)$ (a conforming discretization); that means that you still have $$ a(v_h,v_h) \geq c \|v_h\|_{H^1}^2 >0 \qquad\text{for all }v_h\in V_h\setminus\{0\}. \tag{3}$$
Now for the last step: To transform the variational form to a system of linear equations, you pick a basis $\{\varphi_1,\dots,\varphi_N\}$ of $V_h$, write $u_h =\sum_{i=1}^N u_i\varphi_i$ and insert $v_h=\varphi_j$, $1\leq j\leq N$ into $(2)$. The stiffness matrix $K$ then has the entries $K_{ij}=a(\varphi_i,\varphi_j)$ (which coincides with what you wrote).
Now take an arbitrary vector $\vec v=(v_1,\dots,v_N)^T\in \mathbb{R}^N\setminus\{0\}$ and set $v_h:=\sum_{i=1}^Nv_i \varphi_i\in V_h$. Then we have by $(3)$ and the bilinearity of $a$ (i.e., you can move scalars and sums into both arguments) $$ \vec v^T K \vec v = \sum_{i=1}^N\sum_{j=1}^N v_iK_{ij} v_j = \sum_{i=1}^N\sum_{j=1}^N a(v_i\varphi_i,v_j\varphi_j) = a(v_h,v_h) >0.$$ Since $\vec v\neq 0$ was arbitrary, this implies that $K$ is positive definite.
TL;DR: The stiffness matrix is positive definite because it comes from a conforming discretization of a (self-adjoint) elliptic partial differential equation.
If stiffness of element is not positive, then the system is not stable. So the model is most likely not correct. Look at the most basic equation of harmonic oscillator
$$m x''(t) + k x(t) = f(t)$$
The solution is unstable if $k$ is negative (look at the roots of the characteristic equation). It means the solution will blow up. The stiffness has to be a restoring force. At least for a physical spring. The stiffness matrix extends this to large number of elements (global stiffness matrix). That is all. But it is the same basic idea. FEM basis is in the stiffness matrix method for structural analysis where each element has a stiffness associated with it.
