halbux / sparselizard

C++ FEM library | user-friendly | multi-physics | hp-adaptive | HPC
http://www.sparselizard.org
Other
334 stars 62 forks source link

gettotalforce results do not make sense. #62

Closed DerAlbiG closed 2 years ago

DerAlbiG commented 2 years ago

Hi, i have an axis-symmetric magnetostatic setup of an inductor and iron-core which is slightly offset from the inductor center. The inductor is defined by a current density and the iron core is non-linear. For this, i adapted multiple samples that you provide. (Thank you!). I am cross-checking the simulation-results with commercial software and so far i am immensely satisfied with the results for B and µr (µr due to saturation). I am now trying to:

gettotalforce(coil, h, mu); //176N (expect -97N)
gettotalforce(core, h, mu); //86N (expect 97N)

and here, my luck kind of runs out. In my sparselizard-solution the calculated forces do not make any sense. The coil and core-force should balance to 0 and the numerical values are all over the place.

I suspect that i am doing something fundamentally wrong. Unfortunately i do not know when gettotalforce is applicable and what the alternative would be. So far i have been unable to adapt the force-calculation of the rotating motor example.

Any hint is appreciated, Thank you.

halbux commented 2 years ago

Check this example for mag forces in axisym:

https://github.com/halbux/sparselizard/blob/master/examples/magnetic-force-sandbox-axisymmetry-2d/main.cpp

Gettotalforce is not correct for saturated cases. Boundary integral of Maxwell stress tensor is the way to go.

DerAlbiG commented 2 years ago

Sadly, i have to report, that i have only partial success with this. I did: 1) round the corners of the core, but as it is on-axis as a full cylinder, i can only round the edges that are not on the center-line. Rounded or not seems to be inconsequential, tbh. 2) I introduced the core-skin: mymesh.selectskin(core_skin, core); (and coil-skin too) 3) I adapted the force calculation:

expression T = nu * ( b*transpose(b) - 0.5*b*b * eye(3) );
double fsurcalc = 2*getpi()*compy(on(wholedomain, T)*normal(core)).integrate(core_skin, 5);

which gives -NAN.

Success: The integration over the inductor/coil works, and it actually gives a reasonable result (-93N, where -97 is expected). I did not round the corners of its geometry. (is this advisable?)

Partial success: When i put the core not on the center-line, but move it off-axis (make it a tube) the integration works. If I make the inner diameter of the core close to 0, the value of the core-force is reasonable close to the force on the inductor.

I suspect the fact that core_skin includes the portion of the outline that is on-axis creates the issue (is the normal there 0?). How do I exclude it? I am using the gmsh-api and have the relevant core-outline-part as wire.

I tried to give the outline its own physical region: gmsh::model::addPhysicalGroup(1, {core_outline_excl_axis}, core_skin); But the created entity is not usable in the force-calculation as it complains: "Did you try to compute normal(physreg) on a line not touching the argument region?"

Can you explain why including the axis creates this issue?

DerAlbiG commented 2 years ago

So far i have corrected the issue by making densemat::sum ignoring NANs.

double densemat::sum(void)
{
    double* myvaluesptr = myvalues.get();
    double val = 0;

    for (long long int i = 0; i < numrows*numcols; i++)
        val += std::isnan(myvaluesptr[i]) ? 0.0 : myvaluesptr[i];
    return val;
}

Feedback appreciated.