how is a first variation of a functional calculated?

Assuming that $L$ is a linear operator, we have $$ frac{d}{dt} langle L (phi+teta), phi+teta rangle|_{t=0}= \ frac{d}{dt} langle L phi, phi rangle+t(langle L phi, eta rangle+ langle L eta, phi rangle) +t^2 langle L eta, eta rangle |_{t=0}= \ langle L phi, eta rangle+ langle L eta, phi rangle $$ At this point, you are already at your whit`s end - if you knew your operator (e.g. $L=- Delta$), you could probably perform an integration by parts. Adapting your notation, this means that the first term of $delta F$ is (with the factor of $frac{1}{2}$ in front) $$ frac{1}{2}(langle L phi, delta phi rangle+ langle L deltaphi, phi rangle) $$ for the linear terms, we have $$ frac{d}{dt} langle phi+teta, f rangle|_{t=0}=frac{d}{dt} langle phi,f rangle + t langle eta, f rangle|_{t=0}=langleeta,frangle $$ Switching back to your notation, we have with the factor of $frac{1}{2}$ in front $$ frac{1}{2}langledelta phi,frangle $$ The term $langle f, phi rangle$ is taken care of in the same manner. Adding all terms give you the desired variation $delta F$. However, I wanna note that I computed a Gateaux derivative of a functional (also called "variation" sometimes) and if you knew how $F$ depended on $phi, nabla phi$ explicitly, i.e. $F=F(phi, nabla phi)$, you would arrive at the weak form of the Euler-Lagrange equations. On another note, you still have to take care of the domain of your differential operator $L$, which is often NOT the whole $L^2$ space. I omitted this part, as it would require more informations on $L$.