In Wikipedia, the formula for the variance of the triangular distribution is given here. However, I don't know how to find it. I have tried a brute force method but the formula is quite complicated (polynomial of degree 5 in a, b, c) and I can't simplify it (I tried manually and with Xcas).
The following part is edited thanks to @Imaosome remark:
I came to this question with the following problem: Say $X$ and $Y$ are independent random variables with uniform distribution between $0$ and $1$. I want to study $Z = X-Y$. With convolution, I find that the distribution is triangular, centered in $0$ with extremities $-1$ and $1$ (the proof is also available in this pdf here). With Wikipedia notations, it gives $a=-1, b=1, c=0$. And in this case, we sum $2$ independent variables therefore the variance shoud be $Var(X)+Var(-Y) = Var(X)+Var(Y)=2 Var(X)$. The variance of $X$ is $1/12$ (see for instance formula here). Therefore the variance of $Z$ is $Var(Z) = 2 * 1/12 = 1/6$.
If I come back to Wikipedia formula, I find: $$ \frac{a^2+b^2+c^2-ab-ac-bc}{18}=\frac{1+1+0+1-0-0}{18}=1/6. $$
This shows, at least in that particular case that the formula is correct. But I would like to try and prove Wikipedia result. Once again, I know that it should be possible to prove it by integration but I did not succeed and I hope somebody has a simple way to get this formula.
Here are the details. By definition, we want to compute: \begin{align*} \sigma^2 &= \int_a^c \frac{2(x-a)}{(b-a)(c-a)} \left( x - \frac{a+b+c}{3} \right)^2 dx + \int_c^b \frac{2(b-x)}{(b-a)(b-c)} \left( x - \frac{a+b+c}{3} \right)^2 dx \\ &= \frac{2}{(b-a)(c-a)} \left[ \frac{1}{3} (x-a) \left( x - \frac{a+b+c}{3} \right)^3 - \frac{1}{12} \left( x - \frac{a+b+c}{3} \right)^4 \right]_a^c \\ & ~~~~~~ \frac{2}{(b-a)(b-c)} \left[ \frac{1}{3} (b-x) \left( x - \frac{a+b+c}{3} \right)^3 - \frac{1}{12} \left( x - \frac{a+b+c}{3} \right)^4 \right]_c^b \end{align*}
From there, one can see that terms in $\left( x - \frac{a+b+c}{3}\right)$ cancel out leaving us with: \begin{align*} \sigma^2 &= \frac{2}{12(b-a)(c-a)} \left( - \left( \frac{2c-a-b}{3} \right)^4 + \left( \frac{2a-b-c}{3} \right)^4 \right) \\ & ~~~~~~ \frac{2}{12(b-a)(b-c)} \left( - \left( \frac{2b-a-c}{3} \right)^4 + \left( \frac{2c-a-b}{3} \right)^4 \right) \\ &= \frac{((c-a)-(b-c))^5 +(b-c)((b-a)+(c-a))^4 - (c-a)((b-a)+(b-c))^4}{ 2 \times 3^5 (b-a)(b-c)(c-a) } \end{align*}
And from this point, I am stuck. After all, maybe the last line is not helping much. I don't know.
The computation is best performed on the location-transformed variable such that the mode is at $0$. In particular, suppose $$X \sim \operatorname{Triangular}(a,b,c), \\ f_X(x) = \begin{cases} 0, & x < a \\ \frac{2(x-a)}{(b-a)(c-a)}, & a \le x \le c \\ \frac{2(b-x)}{(b-a)(b-c)}, & c < x \le b, \\ 0, & x > b. \end{cases}$$ Then consider $$Y = X - c$$ which then has density $$f_Y(y) = \begin{cases} 0, & y < a - c \\ \frac{2(y-(a-c))}{(b-a)(c-a)}, & a-c \le y \le 0 \\ \frac{2((b-c)-y)}{(b-a)(b-c)}, & 0 < y \le b-c \\ 0, & x > b-c. \end{cases}$$ All we did was shift the density, which lets us now use the auxiliary parameters $\alpha = a-c$, $\beta = b-c$, to write $$f_Y(y) = \begin{cases} 0, & y < \alpha \\ \frac{2(y-\alpha)}{(\beta - \alpha)(-\alpha)}, & \alpha \le y \le 0 \\ \frac{2(\beta-y)}{(\beta - \alpha)\beta}, & 0 < y \le \beta \\ 0, & y > \beta. \end{cases}$$
Now the variance calculation performed on $Y$ proceeds as follows: $$\begin{align} \operatorname{E}[Y^k] &= \int_{y=\alpha}^0 y^k \frac{2(y-\alpha)}{(\beta-\alpha)(-\alpha)} \, dy + \int_{y=0}^\beta y^k \frac{2(\beta-y)}{(\beta-\alpha)\beta} \, dy \\ &= \frac{2}{(\beta - \alpha)(-\alpha)} \left[\frac{y^{k+2}}{k+2} - \alpha \frac{y^{k+1}}{k+1}\right]_{y=\alpha}^0 + \frac{2}{(\beta - \alpha)\beta} \left[\beta\frac{y^{k+1}}{k+1} - \frac{y^{k+2}}{k+2} \right]_{y=0}^\beta \\ &= \frac{2}{(\beta-\alpha)} (\beta^{k+1} - \alpha^{k+1}) \left(\frac{1}{k+1} - \frac{1}{k+2} \right) \\ &= \frac{2(\beta^{k+1} - \alpha^{k+1})}{(k+1)(k+2)(\beta-\alpha)}. \end{align}$$
For $k \in \{1, 2\}$ we easily get $$\operatorname{E}[Y] = \frac{\alpha + \beta}{3} \\ \operatorname{E}[Y^2] = \frac{\alpha^2 + \alpha\beta + \beta^2}{6}.$$
Consequently, $$\operatorname{Var}[X] = \operatorname{Var}[Y] = \frac{\alpha^2 - \alpha \beta + \beta^2}{18} = \frac{a^2 + b^2 + c^2 - (ab + bc + ca)}{18}.$$
Note that we have used the fact that a location transformation of a random variable does not change its variance. This calculation minimizes the number of integrations we need to perform and reduces the number of parameters we need to account for. We also get the additional result for the $k^{\rm th}$ raw moment of $Y$, which has a particularly convenient form.
Also note that we could have written the density of the transformed variable $Y$ simply by choosing $c = 0$, $a = \alpha$, $b = \beta$ from the density of $X$; however, I wanted to emphasize how the density of $Y$ is obtained by shifting $X$ by $c$.