Select Git revision
gmsh.h_cwrap
-
Christophe Geuzaine authoredChristophe Geuzaine authored
F_Hysteresis.cpp 176.95 KiB
// GetDP - Copyright (C) 1997-2018 P. Dular and C. Geuzaine, University of Liege
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <getdp@onelab.info>.
//
// Contributor(s):
// Johan Gyselinck
// Kevin Jacques
//
#include <math.h>
#include <string.h>
#include "ProData.h"
#include "F.h"
#include "Message.h"
#include <iostream>
#define SQU(a) ((a)*(a))
#define CUB(a) ((a)*(a)*(a))
#define MU0 1.25663706144e-6
#define FLAG_WARNING_INFO_INV 1
#define FLAG_WARNING_INFO_ROOTFINDING 2
#define FLAG_WARNING_INFO_MIN 3
#define FLAG_WARNING_DISP_INV 11
#define FLAG_WARNING_DISP_ROOTFINDING 12
#define FLAG_WARNING_STOP_INV 101
#define FLAG_WARNING_STOP_ROOTFINDING 102
#define FLAG_WARNING_ITER 100
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))
extern struct CurrentData Current ;
/* ------------------------------------------------------------------------ */
/*
Vectorized Jiles-Atherton hysteresis model
J. Gyselinck, P. Dular, N. Sadowski, J. Leite and J.P.A. Bastos,
"Incorporation of a Jiles-Atherton vector hysteresis model in
2D FE magnetic field computations. Application of the Newton-Raphson method",
Vol. 23, No. 3, pp. 685-693, 2004.
*/
double F_Man(double He, double Ms, double a)
{
// Anhysteretic magnetisation
if(fabs(He) < 0.01*a)
//return Ms*He/(3.*a) ; // Aprox. up to 1st order
return Ms*(He/(3.*a)-1/45*CUB(He/a)) ; // Approx. up to 3rd order
else return Ms*(cosh(He/a)/sinh(He/a)-a/He) ;
}
double F_dMandHe(double He, double Ms, double a)
{
// Derivative of the magnetisation Man with respect to the effective field He
if(fabs(He) < 0.01*a)
//return Ms/(3.*a) ; // Aprox. up to 1st order
return Ms/(3.*a)-Ms/(15*a)*SQU(He/a) ; // Approx. up to 3rd order
else return Ms/a*(1-SQU(cosh(He/a)/sinh(He/a))+SQU(a/He)) ;
}
void FV_Man(double He[3], double Ms, double a, double Man[3])
{
double nHe = sqrt(He[0]*He[0]+He[1]*He[1]+He[2]*He[2]) ;
if( !nHe ) {
Man[0] = Man[1] = Man[2]= 0. ;
}
else {
double auxMan = F_Man(nHe, Ms, a) ;