Virtual Work implementation not clear
I want to use the method of virtual works to calculate a force. I know there is a function in GetDP called F_VirtualWork
in file functions/F_Misc.cpp
so I tried to understand it in order to be sure of its usability. Nevertheless, I'm not sure if I do understand how it is implemented. I made some assumptions but they don't match with the code. I couldn't find something that matched with the code. Therefore, I write here my assumptions to see if someone finds an error in them or something in the implementation. To sum up, I don't understand why the DetJac_dx
is computed that way and I don't understand the signs in the final calculation. Down below, what I understood:
The principle of virtual work for the electric or magnetic field can be expressed given the energy:
W = \frac{\epsilon}{2}\int\limits_V \bf{E}\cdot \bf{E} dV
+ \frac{\mu}{2} \int\limits_V \bf{H}\cdot \bf{H} dV
Then virtual force in a direction p
can be expressed as
F_p = \frac{\epsilon}{2}\frac{d}{d \bf{p}} \int_V \bf{E}\cdot \bf{E} d V
+ \frac{\mu}{2}\frac{d }{d \bf{p}} \int_V \bf{H}\cdot \bf{H} d V
If we take the integral term of the electric field (E) and apply the chain rule:
F_{E,p} = \frac{\epsilon}{2}\frac{d }{d\bf{p}} \int\limits_V \bf{E}^2 d V =
\frac{\epsilon}{2}
\left[
\int\limits_V 2\frac{d }{d \bf{p}}\bf{E} d V
+
\int\limits_V \bf{E}^2\frac{d }{d \bf{p}} d V
\right]
Consider then, without coefficients:
\frac{d }{d\bf{p}} \int\limits_V \bf{E}^2 d V =
\int\limits_V 2\frac{d }{d \bf{p}}\bf{E} d V
+
\int\limits_V \bf{E}^2\frac{d }{d \bf{p}} (d V)
= A + B
I understand that the function F_VirtualWork
only takes into account B, so now we consider:
F_{E,p} = \int\limits_V \bf{E}^2\frac{d }{d \bf{p}}( d V )
The contribution of a reference element K_e
to the force in direction \bf{p}
F_{K_e} = \int\limits_{K_e }\bf{E}^2\frac{d }{d \bf{p}} (d K_e)
All articles and info that I could find about virtual work in the electric field take a paper from J. Coulomb and G. Meunier as reference (cited and the end of this text). In that article, the following derivative is defined (and not explained):
\frac{d }{d \bf{s}} (dv) = |G|^{-1} \frac{d |G|}{d \bf{s}} \cdot dv
I suppose this expression is used in the code to express \frac{d }{d \bf{p}} (d K_e)
using DetJac_dx[0]/DetJac
, DetJac_dx[1]/DetJac
, DetJac_dx[2]/DetJac
. But I don't really get how DetJac_dx[i]
come from. I can deduce that this quantity can come from the solution of a system like:
grad(DetJac) = grad(node\_basis\_functions) Adj(Jac)
Where grad()
is the gradient and Adj()
is the classical adjoint. But I can't understand where this comes from.
Then we need to compute square of the field, {E}=(E_{1},E_{2},E_{3})=E_{1}i+E_{2}j+E_{3}k
. If we expressed it in its coordinates:
{E}^2 = (E_{1}i+E_{2}j+E_{3}k)(E_{1}i+E_{2}j+E_{3}k) =E_{1}^2ii+E_{2}^2jj+E_{3}^kk+2E_{1}E_{2}ij+2E_{1}E_{3}ik+2E_{2}E_{3}jk
Considering DetJac_dx/DetJac = D
, we have {D}=(D_{1},D_{2},D_{3})=D_{1}i+D_{2}j+D_{3}k
. To get the force contribution F =(F_{1},F_{2},F_{3}) = {E}^2 D
:
F = (D_{1}i+D_{2}j+D_{3}k)(E_{1}^2ii+E_{2}^2jj+E_{3}^kk+2E_{1}E_{2}ij+2E_{1}E_{3}ik+2E_{2}E_{3}jk)
For F_a = D_{1}i{E}^2
:
F_a = E_{1}^2 D_{1}iii+ E_{2}^2 D_{1}ijj+E_{3}^2 D_{1}ikk+ 2E_{1}E_{2} D_{1}iij + 2E_{1} E_{3} D_{1} iik + 2E_{2} E_{3} D_{1} ijk
Computing the directions:
F_a = E_{1}^2 D_{1}i+E_{2}^2 D_{1}i+E_{3}^2 D_{1}i+2E_{1}E_{2}D_{1}j+2E_{1}E_{3}D_{1}k
For F_b =D_{2}i{E}^2
:
F_b = E_{1}^2D_{2}jii+E_{2}^2D_{2}jjj+E_{3}^2 D_{2}jkk+2E_{1}E_{2}D_{2}jij+2E_{1}E_{3}D_{2}jik+2E_{2}E_{3}D_{2}jjk
Computing the directions:
F_b = E_{1}^2 D_{1}j+E_{2}^2 D_{1}j+E_{3}^2 D_{1}j+2E_{1}E_{2}D_{2}i+2E_{2}E_{3}D_{2}k
For F_c =D_{3}i{E}^2
:
F_c = E_{1}^2D_{3}kii+E_{2}^2D_{3}kjj+E_{3}^2 D_{3}kkk+2E_{1}E_{2}D_{3}kij+2E_{1}E_{3}D_{3}kik+2E_{2}E_{3}D_{3}kjk
Computing the directions:
F_c = E_{1}^2 D_{1}j+E_{2}^2 D_{1}j+E_{3}^2 D_{1}j+2E_{1}E_{3}D_{3}i+2E_{2}E_{3}D_{3}j
Then we arrange F = F_a+F_b+F_c = F_{1}i+F_{2}j+F_{3}k
and we have:
F_1i = E_{1}^2 D_{1}i+E_{2}^2 D_{1}i+E_{3}^2 D_{1}i+2E_{1}E_{2}D_{2}i+2E_{1}E_{3}D_{3}i
F_2j = E_{1}^2 D_{1}j+E_{2}^2 D_{1}j+E_{3}^2 D_{1}j+2E_{1}E_{2}D_{1}j+2E_{2}E_{3}D_{3}j
F_3k = E_{1}^2 D_{1}k+E_{2}^2 D_{1}k+E_{3}^2 D_{1}k+2E_{1}E_{3}D_{1}k+2E_{2}E_{3}D_{2}k
The code in F_VirtualWork
is similar, except for the signs. The only crazy idea that comes to my mind is that we take the sign of i,j,k
as negative so we have (i,j,k) = (-i_{new},-j_{new},-k_{new})
:
s[0] = (DetJac_dx[0] * (-squF[0] + squF[1] + squF[2]) -
2 * DetJac_dx[1] * squF[3] - 2 * DetJac_dx[2] * squF[5]) /
DetJac;
s[1] = (DetJac_dx[1] * (squF[0] - squF[1] + squF[2]) -
2 * DetJac_dx[0] * squF[3] - 2 * DetJac_dx[2] * squF[4]) /
DetJac;`
s[2] = (DetJac_dx[2] * (squF[0] + squF[1] - squF[2]) -
2 * DetJac_dx[0] * squF[5] - 2 * DetJac_dx[1] * squF[4]) /
DetJac;
J. Coulomb and G. Meunier, "Finite element implementation of virtual work principle for magnetic or electric force and torque computation," in IEEE Transactions on Magnetics, vol. 20, no. 5, pp. 1894-1896, September 1984, doi: 10.1109/TMAG.1984.1063232.