problem.cpp

一个2D电磁场FEM计算的VC++源程序

#include "stdafx.h"
#include "problem.h"
#include "fullmatrix.h"

#define ElementsPerSkinDepth 10

/////////////////////////////////////////////////////////////////////////////
// CNode construction

CNode::CNode()
{
	x=0.;
	y=0.;
	IsSelected=FALSE;
	BoundaryMarker=-1;
}

CComplex CNode::CC()
{
	return (x+I*y);
}

double CNode::GetDistance(double xo, double yo)
{
	return sqrt((x-xo)*(x-xo) + (y-yo)*(y-yo));
}

void CNode::ToggleSelect()
{
	if (IsSelected==TRUE) IsSelected=FALSE;
	else IsSelected=TRUE;
}

/////////////////////////////////////////////////////////////////////////////
// CMeshNode construction

CMeshNode::CMeshNode()
{
	x=0.; y=0.;
	A.re=0; A.im=0;
	msk=0;
}

CComplex CMeshNode::CC()
{
	return (x+I*y);
}

double CMeshNode::GetDistance(double xo, double yo)
{
	return sqrt((x-xo)*(x-xo) + (y-yo)*(y-yo));
}

/////////////////////////////////////////////////////////////////////////////
// CSegment construction

CSegment::CSegment()
{
	n0=0;
	n1=0;
	IsSelected=FALSE;
	BoundaryMarker=-1;
}

void CSegment::ToggleSelect()
{
	if (IsSelected==TRUE) IsSelected=FALSE;
	else IsSelected=TRUE;
}
/////////////////////////////////////////////////////////////////////////////
// CArcSegment construction

CArcSegment::CArcSegment()
{
	n0=0;
	n1=0;
	IsSelected=FALSE;
	MaxSideLength=-1;
	ArcLength=90.;
	BoundaryMarker=-1;
}

void CArcSegment::ToggleSelect()
{
	if (IsSelected==TRUE) IsSelected=FALSE;
	else IsSelected=TRUE;
}

/////////////////////////////////////////////////////////////////////////////
// CNode construction

CBlockLabel::CBlockLabel()
{
	x=0.;
	y=0.;
	MaxArea=0.;
	IsSelected=FALSE;
	InGroup=0;
	InCircuit=0;
	BlockType=-1;
	IsExternal=FALSE;

	Case=0;
	dVolts=0.;
	J=0.;
	FillFactor=1;
	o=0;
	mu=0;
}

void CBlockLabel::ToggleSelect()
{
	if (IsSelected==TRUE) IsSelected=FALSE;
	else IsSelected=TRUE;
}

double CBlockLabel::GetDistance(double xo, double yo)
{
	return sqrt((x-xo)*(x-xo) + (y-yo)*(y-yo));
}

CMaterialProp::CMaterialProp()
{		
		BlockName="New Material";
		mu_x=1.;
		mu_y=1.;			// permeabilities, relative
		BHpoints=0;		
		Bdata=NULL;
		Hdata=NULL;
		slope=NULL;
		H_c=0.;				// magnetization, A/m
		Jr=0.;
		Ji=0.;				// applied current density, MA/m^2
		Cduct=0.;		    // conductivity of the material, MS/m
		Lam_d=0.;			// lamination thickness, mm
		Theta_hn=0.;		// hysteresis angle, degrees
		Theta_hx=0.;		// hysteresis angle, degrees
		Theta_hy=0.;		// hysteresis angle, degrees
		NStrands=0;			// number of strands per wire
		WireD=0;
		LamFill=1.;			// lamination fill factor;
		LamType=0;			// type of lamination;
}

CMaterialProp::~CMaterialProp()
{
	if (Bdata!=NULL)	free(Bdata);
	if (Hdata!=NULL)	free(Hdata);
	if (slope!=NULL)	free(slope);
}

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;
	CComplexFullMatrix 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));
		}

	}

	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;
		}

	}

	free(bn);
	free(hn);
	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;
		double res=0;
		double Relax=1;
		CComplex mu,vo,vi,c,H;
		CComplex Md,Mo;
		BOOL Converged=FALSE;
		
		// 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(CComplex(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));