588 lines
14 KiB
C++
588 lines
14 KiB
C++
#include<stdafx.h>
|
|
#include<afxtempl.h>
|
|
#include<stdio.h>
|
|
#include<math.h>
|
|
#include "fkn.h"
|
|
#include "fknDlg.h"
|
|
#include "complex.h"
|
|
#include "mesh.h"
|
|
#include "spars.h"
|
|
#include "fullmatrix.h"
|
|
#include "FemmeDocCore.h"
|
|
|
|
#define ElementsPerSkinDepth 10
|
|
|
|
CMaterialProp::~CMaterialProp()
|
|
{
|
|
if (Bdata!=NULL) free(Bdata);
|
|
if (Hdata!=NULL) free(Hdata);
|
|
if (slope!=NULL) free(slope);
|
|
}
|
|
|
|
CMaterialProp::CMaterialProp()
|
|
{
|
|
slope=NULL;
|
|
Bdata=NULL;
|
|
Hdata=NULL;
|
|
}
|
|
|
|
void CMaterialProp::GetSlopes()
|
|
{
|
|
GetSlopes(0);
|
|
}
|
|
|
|
void CMaterialProp::GetSlopes(double omega)
|
|
{
|
|
if (BHpoints==0) return; // catch trivial case;
|
|
if (slope!=NULL) return; // already have computed the slopes;
|
|
|
|
int i,k;
|
|
BOOL CurveOK=FALSE;
|
|
BOOL ProcessedLams=FALSE;
|
|
CFullMatrix L;
|
|
double l1,l2;
|
|
CComplex *hn;
|
|
double *bn;
|
|
CComplex mu;
|
|
|
|
L.Create(BHpoints);
|
|
bn =(double *) calloc(BHpoints,sizeof(double));
|
|
hn =(CComplex *)calloc(BHpoints,sizeof(CComplex));
|
|
slope=(CComplex *)calloc(BHpoints,sizeof(CComplex));
|
|
|
|
|
|
// strip off some info that we can use during the first
|
|
// nonlinear iteration;
|
|
mu_x = Bdata[1]/(muo*abs(Hdata[1]));
|
|
mu_y = mu_x;
|
|
Theta_hx=Theta_hn;
|
|
Theta_hy=Theta_hn;
|
|
|
|
// first, we need to doctor the curve if the problem is
|
|
// being evaluated at a nonzero frequency.
|
|
if(omega>0)
|
|
{
|
|
// Make an effective B-H curve for harmonic problems.
|
|
// this one convolves B(H) where H is sinusoidal with
|
|
// a sine at the same frequency to get the effective
|
|
// amplitude of B
|
|
double munow,mumax=0;
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
hn[i]=Hdata[i];
|
|
for(k=1,bn[i]=0;k<=i;k++)
|
|
{
|
|
bn[i]+=Re((4.*(Hdata[k]*Bdata[k-1] -
|
|
Hdata[k-1]*Bdata[k])*(-cos((Hdata[k-1]*PI)/(2.*Hdata[i])) +
|
|
cos((Hdata[k]*PI)/(2.*Hdata[i]))) + (-Bdata[k-1] +
|
|
Bdata[k])*((Hdata[k-1] - Hdata[k])*PI +
|
|
Hdata[i]*(-sin((Hdata[k-1]*PI)/Hdata[i]) +
|
|
sin((Hdata[k]*PI)/Hdata[i]))))/ ((Hdata[k-1] - Hdata[k])*PI));
|
|
}
|
|
}
|
|
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
Bdata[i]=bn[i];
|
|
Hdata[i]=hn[i];
|
|
munow=Re(Bdata[i]/Hdata[i]);
|
|
if (munow>mumax) mumax=munow;
|
|
}
|
|
|
|
// apply complex permeability to approximate the
|
|
// effects of hysteresis. We will follow a kludge
|
|
// suggested by O'Kelly where hysteresis angle
|
|
// is proportional to permeability. This implies
|
|
// that loss goes with B^2
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
Hdata[i]*=exp(I*Bdata[i]*Theta_hn*DEG/(Hdata[i]*mumax));
|
|
}
|
|
|
|
MuMax = mumax/muo;
|
|
}
|
|
|
|
while(CurveOK!=TRUE)
|
|
{
|
|
// make sure that the space for computing slopes is cleared out
|
|
L.Wipe();
|
|
|
|
// impose natural BC on the `left'
|
|
l1=Bdata[1]-Bdata[0];
|
|
L.M[0][0]=4./l1;
|
|
L.M[0][1]=2./l1;
|
|
L.b[0]=6.*(Hdata[1]-Hdata[0])/(l1*l1);
|
|
|
|
// impose natural BC on the `right'
|
|
int n=BHpoints;
|
|
l1=Bdata[n-1]-Bdata[n-2];
|
|
L.M[n-1][n-1]=4./l1;
|
|
L.M[n-1][n-2]=2./l1;
|
|
L.b[n-1]=6.*(Hdata[n-1]-Hdata[n-2])/(l1*l1);
|
|
|
|
for(i=1;i<BHpoints-1;i++)
|
|
{
|
|
l1=Bdata[i]-Bdata[i-1];
|
|
l2=Bdata[i+1]-Bdata[i];
|
|
|
|
L.M[i][i-1]=2./l1;
|
|
L.M[i][i]=4.*(l1+l2)/(l1*l2);
|
|
L.M[i][i+1]=2./l2;
|
|
|
|
L.b[i]=6.*(Hdata[i]-Hdata[i-1])/(l1*l1) +
|
|
6.*(Hdata[i+1]-Hdata[i])/(l2*l2);
|
|
}
|
|
|
|
L.GaussSolve();
|
|
for(i=0;i<BHpoints;i++) slope[i]=L.b[i];
|
|
|
|
// now, test to see if there are any "bad" segments in there.
|
|
// it is probably sufficient to do this test just on the
|
|
// real part of the BH curve...
|
|
for(i=1,CurveOK=TRUE;i<BHpoints;i++)
|
|
{
|
|
double L,c0,c1,c2,d0,d1,u0,u1,X0,X1;
|
|
|
|
// it is probably sufficient to do this test just on the
|
|
// real part of the BH curve. We do the test on just the
|
|
// real part of the curve by taking the real parts of
|
|
// slope and Hdata
|
|
d0=slope[i-1].re;
|
|
d1=slope[i].re;
|
|
u0=Hdata[i-1].re;
|
|
u1=Hdata[i].re;
|
|
L=Bdata[i]-Bdata[i-1];
|
|
|
|
c0=d0;
|
|
c1= -(2.*(2.*d0*L + d1*L + 3.*u0 - 3.*u1))/(L*L);
|
|
c2= (3.*(d0*L + d1*L + 2.*u0 - 2.*u1))/(L*L*L);
|
|
X0=-1.;
|
|
X1=-1.;
|
|
|
|
u0=c1*c1-4.*c0*c2;
|
|
// check out degenerate cases
|
|
if(c2==0)
|
|
{
|
|
if(c1!=0) X0=-c0/c1;
|
|
}
|
|
else if(u0>0)
|
|
{
|
|
u0=sqrt(u0);
|
|
X0=-(c1 + u0)/(2.*c2);
|
|
X1=(-c1 + u0)/(2.*c2);
|
|
}
|
|
|
|
//now, see if we've struck gold!
|
|
if (((X0>=0.)&&(X0<=L))||((X1>=0.)&&(X1<=L)))
|
|
CurveOK=FALSE;
|
|
}
|
|
|
|
if(CurveOK!=TRUE) //remedial action
|
|
{
|
|
// Smooth out input points
|
|
// to get rid of rapid transitions;
|
|
// Uses a 3-point moving average
|
|
for(i=1;i<BHpoints-1;i++)
|
|
{
|
|
bn[i]=(Bdata[i-1]+Bdata[i]+Bdata[i+1])/3.;
|
|
hn[i]=(Hdata[i-1]+Hdata[i]+Hdata[i+1])/3.;
|
|
}
|
|
|
|
for(i=1;i<BHpoints-1;i++){
|
|
Hdata[i]=hn[i];
|
|
Bdata[i]=bn[i];
|
|
}
|
|
}
|
|
|
|
|
|
if((CurveOK==TRUE) && (ProcessedLams==FALSE))
|
|
{
|
|
// if the material is laminated and has a non-zero conductivity,
|
|
// we have to do some work to find the "right" apparent BH
|
|
// curve for the material for AC problems. Essentially, we have
|
|
// to solve a 1-D finite element problem at each B-H point to
|
|
// figure out how much flux the lamination is really carrying.
|
|
if((omega>0) && (Lam_d!=0) && (Cduct!=0))
|
|
{
|
|
// Calculate a new apparent b and h for each point on the
|
|
// curve to account for the laminations.
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
mu=LaminatedBH(omega,i);
|
|
bn[i]=abs(mu*Hdata[i]);
|
|
hn[i]=bn[i]/mu;
|
|
}
|
|
|
|
// Replace the original BH points with the apparent ones
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
Bdata[i]=bn[i];
|
|
Hdata[i]=hn[i];
|
|
}
|
|
|
|
// Make the program check the consistency and recompute the
|
|
// proper slopes of the new BH curve one more time;
|
|
CurveOK=FALSE;
|
|
}
|
|
|
|
// take care of LamType=0 situation by changing the apparent B-H curve.
|
|
if((LamType==0) && (LamFill!=1))
|
|
{
|
|
// this is changed from the previous version so that
|
|
// a complex-valued H can be properly accommodated
|
|
// in the algorithm.
|
|
for(i=1;i<BHpoints;i++)
|
|
{
|
|
mu=LamFill*Bdata[i]/Hdata[i]+(1.-LamFill)*muo;
|
|
Bdata[i]=abs(mu*Hdata[i]);
|
|
Hdata[i]=Bdata[i]/mu;
|
|
}
|
|
// Make the program check the consistency and recompute the
|
|
// proper slopes of the new BH curve one more time;
|
|
CurveOK=FALSE;
|
|
}
|
|
|
|
ProcessedLams=TRUE;
|
|
}
|
|
|
|
}
|
|
|
|
|
|
// FILE *fp;
|
|
// fp=fopen("c:\\temp\\bhjunk.txt","wt");
|
|
// for(i=1;i<BHpoints;i++)
|
|
// fprintf(fp,"%g %g %g\n",abs(Hdata[i]),Re(Bdata[i]*abs(Hdata[i])/Hdata[i]),Im(Bdata[i]*abs(Hdata[i])/Hdata[i]));
|
|
// fclose(fp);
|
|
|
|
free(bn);
|
|
free(hn);
|
|
return;
|
|
}
|
|
|
|
|
|
CComplex CMaterialProp::GetH(double B)
|
|
{
|
|
double b,z,z2,l;
|
|
CComplex h;
|
|
int i;
|
|
|
|
b=fabs(B);
|
|
|
|
if(BHpoints==0) return CComplex(b/(mu_x*muo));
|
|
|
|
if(b>Bdata[BHpoints-1])
|
|
return (Hdata[BHpoints-1] + slope[BHpoints-1]*(b-Bdata[BHpoints-1]));
|
|
|
|
for(i=0;i<BHpoints-1;i++)
|
|
if((b>=Bdata[i]) && (b<=Bdata[i+1])){
|
|
l=(Bdata[i+1]-Bdata[i]);
|
|
z=(b-Bdata[i])/l;
|
|
z2=z*z;
|
|
h=(1.-3.*z2+2.*z2*z)*Hdata[i] +
|
|
z*(1.-2.*z+z2)*l*slope[i] +
|
|
z2*(3.-2.*z)*Hdata[i+1] +
|
|
z2*(z-1.)*l*slope[i+1];
|
|
return h;
|
|
}
|
|
|
|
return CComplex(0);
|
|
}
|
|
|
|
CComplex CMaterialProp::GetdHdB(double B)
|
|
{
|
|
double b,z,l;
|
|
CComplex h;
|
|
int i;
|
|
|
|
b=fabs(B);
|
|
|
|
if(BHpoints==0) return CComplex(b/(mu_x*muo));
|
|
|
|
if(b>Bdata[BHpoints-1])
|
|
return slope[BHpoints-1];
|
|
|
|
for(i=0;i<BHpoints-1;i++)
|
|
if((b>=Bdata[i]) && (b<=Bdata[i+1])){
|
|
l=(Bdata[i+1]-Bdata[i]);
|
|
z=(b-Bdata[i])/l;
|
|
h=6.*z*(z-1.)*Hdata[i]/l +
|
|
(1.-4.*z+3.*z*z)*slope[i] +
|
|
6.*z*(1.-z)*Hdata[i+1]/l +
|
|
z*(3.*z-2.)*slope[i+1];
|
|
return h;
|
|
}
|
|
|
|
return CComplex(0);
|
|
}
|
|
|
|
CComplex CMaterialProp::Get_v(double B)
|
|
{
|
|
if (B==0) return slope[0];
|
|
|
|
return (GetH(B)/B);
|
|
}
|
|
|
|
CComplex CMaterialProp::Get_dvB2(double B)
|
|
{
|
|
if (B==0) return 0;
|
|
|
|
return 0.5*(GetdHdB(B)/(B*B) - GetH(B)/(B*B*B));
|
|
}
|
|
|
|
void CMaterialProp::GetBHProps(double B, double &v, double &dv)
|
|
{
|
|
// version to use in the magnetostatic case in
|
|
// which we know that v and dv ought to be real-valued.
|
|
CComplex vc,dvc;
|
|
|
|
GetBHProps(B,vc,dvc);
|
|
v =Re(vc);
|
|
dv=Re(dvc);
|
|
}
|
|
|
|
void CMaterialProp::GetBHProps(double B, CComplex &v, CComplex &dv)
|
|
{
|
|
double b,z,z2,l;
|
|
CComplex h,dh;
|
|
int i;
|
|
|
|
b=fabs(B);
|
|
|
|
if(BHpoints==0){
|
|
v=mu_x;
|
|
dv=0;
|
|
return;
|
|
}
|
|
|
|
if(b==0){
|
|
v=slope[0];
|
|
dv=0;
|
|
return;
|
|
}
|
|
|
|
if(b>Bdata[BHpoints-1]){
|
|
h=(Hdata[BHpoints-1] + slope[BHpoints-1]*(b-Bdata[BHpoints-1]));
|
|
dh=slope[BHpoints-1];
|
|
v=h/b;
|
|
dv=0.5*(dh/(b*b) - h/(b*b*b));
|
|
return;
|
|
}
|
|
|
|
for(i=0;i<BHpoints-1;i++)
|
|
if((b>=Bdata[i]) && (b<=Bdata[i+1])){
|
|
l=(Bdata[i+1]-Bdata[i]);
|
|
z=(b-Bdata[i])/l;
|
|
z2=z*z;
|
|
h=(1.-3.*z2+2.*z2*z)*Hdata[i] +
|
|
z*(1.-2.*z+z2)*l*slope[i] +
|
|
z2*(3.-2.*z)*Hdata[i+1] +
|
|
z2*(z-1.)*l*slope[i+1];
|
|
dh=6.*z*(z-1.)*Hdata[i]/l +
|
|
(1.-4.*z+3.*z*z)*slope[i] +
|
|
6.*z*(1.-z)*Hdata[i+1]/l +
|
|
z*(3.*z-2.)*slope[i+1];
|
|
v=h/b;
|
|
dv=0.5*(dh/(b*b) - h/(b*b*b));
|
|
return;
|
|
}
|
|
}
|
|
|
|
CComplex CMaterialProp::LaminatedBH(double w, int i)
|
|
{
|
|
int k,n,iter=0;
|
|
CComplex *m0,*m1,*b,*x;
|
|
double L,o,d,ds,B,lastres,res;
|
|
double Relax=1;
|
|
CComplex mu,vo,vi,c,H;
|
|
CComplex Md,Mo;
|
|
BOOL Converged=FALSE;
|
|
res=0;
|
|
|
|
// Base the required element spacing on the skin depth
|
|
// at the surface of the material
|
|
mu=Bdata[i]/Hdata[i];
|
|
o=Cduct*1.e6;
|
|
d=(Lam_d*0.001)/2.;
|
|
ds=sqrt(2/(w*o*abs(mu)));
|
|
n= ElementsPerSkinDepth * ((int) ceil(d/ds));
|
|
L=d/((double) n);
|
|
|
|
x =(CComplex *)calloc(n+1,sizeof(CComplex));
|
|
b =(CComplex *)calloc(n+1,sizeof(CComplex));
|
|
m0=(CComplex *)calloc(n+1,sizeof(CComplex));
|
|
m1=(CComplex *)calloc(n+1,sizeof(CComplex));
|
|
|
|
do{
|
|
// make sure that the old stuff is wiped out;
|
|
for(k=0;k<=n;k++)
|
|
{
|
|
m0[k]=0;
|
|
m1[k]=0;
|
|
b[k] =0;
|
|
}
|
|
|
|
// build matrix
|
|
for(k=0;k<n;k++)
|
|
{
|
|
if(iter!=0)
|
|
{
|
|
B=abs(x[k+1]-x[k])/L;
|
|
vi=GetdHdB(B);
|
|
vo=GetH(B)/B;
|
|
}
|
|
else
|
|
{
|
|
vo=1./mu;
|
|
vi=1./mu;
|
|
}
|
|
|
|
Md = ( (vi+vo)/(2.*L) + I*w*o*L/4.);
|
|
Mo = (-(vi+vo)/(2.*L) + I*w*o*L/4.);
|
|
m0[k] += Md;
|
|
m1[k] += Mo;
|
|
m0[k+1] += Md;
|
|
|
|
Md = ( (vi-vo)/(2.*L));
|
|
Mo = (-(vi-vo)/(2.*L));
|
|
b[k] += (Md*x[k] + Mo*x[k+1]);
|
|
b[k+1] += (Mo*x[k] + Md*x[k+1]);
|
|
}
|
|
|
|
// set boundary conditions
|
|
m1[0] = 0;
|
|
b[0] = 0;
|
|
b[n] += Hdata[i];
|
|
|
|
// solve tridiagonal equation
|
|
// solution ends up in b;
|
|
for(k = 0; k < n; k++)
|
|
{
|
|
c = m1[k]/m0[k];
|
|
m0[k+1] -= m1[k]*c;
|
|
b[k+1] -= b[k]*c;
|
|
}
|
|
b[n] = b[n]/m0[n];
|
|
for(k = n-1; k >= 0; k--)
|
|
b[k] = (b[k] - m1[k]*b[k + 1])/m0[k];
|
|
|
|
iter++;
|
|
|
|
lastres=res;
|
|
res=abs(b[n]-x[n])/d;
|
|
|
|
if (res<1.e-8) Converged=TRUE;
|
|
|
|
// Do the same relaxation scheme as is implemented
|
|
// in the solver to make sure that this effective
|
|
// lamination permeability calculation converges
|
|
if(iter>5)
|
|
{
|
|
if ((res>lastres) && (Relax>0.1)) Relax/=2.;
|
|
else Relax+= 0.1 * (1. - Relax);
|
|
}
|
|
|
|
for(k=0;k<=n;k++) x[k]=Relax*b[k]+(1.0-Relax)*x[k];
|
|
|
|
}while(Converged!=TRUE);
|
|
|
|
mu = x[n]/(Hdata[i]*d);
|
|
|
|
free(x );
|
|
free(b );
|
|
free(m0);
|
|
free(m1);
|
|
|
|
return mu;
|
|
}
|
|
|
|
// get incremental permeability of a nonlinear material for use in
|
|
// incremental permeability formulation about DC offset
|
|
void CMaterialProp::IncrementalPermeability(double B, double w, CComplex &mu1, CComplex &mu2)
|
|
{
|
|
// B == flux density in Tesla
|
|
// w == frequency in rad/s
|
|
|
|
double muinc,murel;
|
|
CComplex k;
|
|
double mu;
|
|
|
|
// get incremental permeability of the DC material
|
|
// (i.e. incremental permeability at the offset)
|
|
muinc=1./(muo*Re(GetdHdB(B)));
|
|
murel=1./(muo*Re(Get_v(B)));
|
|
|
|
// if material is not laminated, just apply hysteresis lag...
|
|
if ((Lam_d==0) || (LamFill==0))
|
|
{
|
|
mu1=muinc*exp(-I*Theta_hn*DEG*muinc/MuMax);
|
|
mu2=murel*exp(-I*Theta_hn*DEG*murel/MuMax);
|
|
return;
|
|
}
|
|
|
|
// crap. Need to make an equivalent permeability that rolls in the effects of laminated
|
|
// eddy currents, using the incremental permeability as the basis for creating the impedance.
|
|
// this can get annoying because we need to back out the iron portion of the permeability
|
|
// in the lamfill<1 case...
|
|
|
|
|
|
if (Cduct!=0)
|
|
{
|
|
CComplex deg45; deg45=1+I;
|
|
CComplex K,halflag;
|
|
double ds;
|
|
|
|
// incremental permeability direction
|
|
mu = (muinc - (1.-LamFill))/LamFill;
|
|
halflag=exp(-I*Theta_hn*DEG*mu/(2.*MuMax));
|
|
ds=sqrt(2./(0.4*PI*w*Cduct*mu));
|
|
K=halflag*deg45*Lam_d*0.001/(2.*ds);
|
|
mu1=(LamFill*mu*tanh(K)/K + (1.- LamFill));
|
|
|
|
// normal permeability direction
|
|
mu = (murel - (1.-LamFill))/LamFill;
|
|
halflag=exp(-I*Theta_hn*DEG*mu/(2.*MuMax));
|
|
ds=sqrt(2./(0.4*PI*w*Cduct*mu));
|
|
K=halflag*deg45*Lam_d*0.001/(2.*ds);
|
|
mu2=(LamFill*mu*tanh(K)/K + (1.- LamFill));
|
|
return;
|
|
}
|
|
else{
|
|
// incremental permeability direction
|
|
mu = (muinc - (1.-LamFill))/LamFill;
|
|
mu1=(mu*exp(-I*Theta_hn*DEG*mu/MuMax)*LamFill + (1.-LamFill));
|
|
|
|
// normal permeability direction
|
|
mu = (murel - (1.-LamFill))/LamFill;
|
|
mu2=(mu*exp(-I*Theta_hn*DEG*mu/MuMax)*LamFill + (1.-LamFill));
|
|
return;
|
|
}
|
|
}
|
|
|
|
void CMaterialProp::IncrementalPermeability(double B, double &mu1, double &mu2)
|
|
{
|
|
// B == flux density in Tesla
|
|
|
|
double muinc, murel;
|
|
|
|
// get incremental permeability of the DC material
|
|
// (i.e. incremental permeability at the offset)
|
|
muinc = 1. / (muo*Re(GetdHdB(B)));
|
|
murel = 1. / (muo*Re(Get_v(B)));
|
|
|
|
// if material is not laminated, just return
|
|
if ((Lam_d == 0) || (LamFill == 0)){
|
|
mu1 = muinc;
|
|
mu2 = murel;
|
|
return;
|
|
}
|
|
|
|
// incremental permeability direction
|
|
mu1 = (muinc*LamFill + (1. - LamFill));
|
|
|
|
// normal permeability direction
|
|
mu2 = (murel*LamFill + (1. - LamFill));
|
|
|
|
return;
|
|
}
|
|
|