📄 prob4big.cpp
字号:
#include<stdafx.h>
#include<stdio.h>
#include<math.h>
#include "fkn.h"
#include "fknDlg.h"
#include "complex.h"
#include "mesh.h"
#include "spars.h"
#include "FemmeDocCore.h"
#define Log log
// #define NEWTON
BOOL CFemmeDocCore::HarmonicAxisymmetric(CBigComplexLinProb &L)
{
int i,j,k,s,flag,sdi_iter,sdin,ww,Iter=0;
int pctr;
CComplex Mx[3][3],My[3][3],Mn[3][3],Me[3][3],be[3]; // element matrices;
double l[3],p[3],q[3]; // element shape parameters;
int n[3]; // numbers of nodes for a particular element;
double a,r,t,x,y,B,w,res,lastres,ds,R,rn[3],g[3],a_hat,R_hat,vol,Cduct;
CComplex K,mu,dv,B1,B2,v[3],u[3],mu1,mu2,lag,halflag,deg45,Jv;
CComplex **Mu,*V_sdi,*V_old;
double c=PI*4.e-05;
double units[]={2.54,0.1,1.,100.,0.00254,1.e-04};
CElement *El;
BOOL LinearFlag=TRUE;
BOOL SDIflag=FALSE;
res=0;
#ifndef NEWTON
CComplex murel,muinc;
#endif
extRo*=units[LengthUnits];
extRi*=units[LengthUnits];
extZo*=units[LengthUnits];
deg45=1+I;
w=Frequency*2.*PI;
CComplex *CircInt1,*CircInt2,*CircInt3;
// Can't handle LamType==1 or LamType==2 in AC problems.
// Detect if this is being attempted.
for(i=0;i<NumEls;i++)
{
if( (blockproplist[meshele[i].blk].LamType==1) ||
(blockproplist[meshele[i].blk].LamType==2) ){
AfxMessageBox("On-edge lamination not supported in AC analyses");
return FALSE;
}
}
// Go through and evaluate permeability for regions subject to prox effects
for(i=0;i<NumBlockLabels;i++) GetFillFactor(i);
V_old=(CComplex *) calloc(NumNodes+NumCircProps,sizeof(CComplex));
// check to see if any circuits have been defined and process them;
if (NumCircProps>0)
{
CircInt1=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
CircInt2=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
CircInt3=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
for(i=0;i<NumEls;i++){
if(meshele[i].lbl>=0)
if(labellist[meshele[i].lbl].InCircuit!=-1){
El=&meshele[i];
// get element area;
for(k=0;k<3;k++) n[k]=El->p[k];
p[0]=meshnode[n[1]].y - meshnode[n[2]].y;
p[1]=meshnode[n[2]].y - meshnode[n[0]].y;
p[2]=meshnode[n[0]].y - meshnode[n[1]].y;
q[0]=meshnode[n[2]].x - meshnode[n[1]].x;
q[1]=meshnode[n[0]].x - meshnode[n[2]].x;
q[2]=meshnode[n[1]].x - meshnode[n[0]].x;
a=(p[0]*q[1]-p[1]*q[0])/2.;
r=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;
// if coils are wound, they act like they have
// a zero "bulk" conductivity...
Cduct=blockproplist[El->blk].Cduct;
if (abs(labellist[El->lbl].Turns)!=1) Cduct=0;
// evaluate integrals;
// total cross-section of circuit;
CircInt1[labellist[El->lbl].InCircuit]+=a;
// integral of conductivity / R over the circuit;
CircInt2[labellist[El->lbl].InCircuit]+=a*Cduct/(0.01*r);
// integral of applied J over current;
CircInt3[labellist[El->lbl].InCircuit]+=
(blockproplist[El->blk].Jr+I*blockproplist[El->blk].Ji)*a*100.;
}
}
for(i=0;i<NumCircProps;i++)
{
if (circproplist[i].CircType==0) // specified current
{
if(CircInt2[i]==0){ //circuit composed of zero cond. materials
circproplist[i].Case=1;
if (CircInt1[i]==0.) circproplist[i].J=0.;
else circproplist[i].J=0.01*(
(circproplist[i].Amps_re+I*circproplist[i].Amps_im) -
CircInt3[i])/CircInt1[i];
}
else{
circproplist[i].Case=2; // need to include an extra
// entry in matrix to solve for
// voltage grad in the circuit
}
}
else{
// case where voltage gradient is specified a priori...
circproplist[i].Case=0;
circproplist[i].dV=circproplist[i].dVolts_re +
I*circproplist[i].dVolts_im;
}
}
}
// check to see if there are any SDI boundaries...
// lineproplist[ meshele[i].e[j] ].BdryFormat==0
for(i=0;i<NumLineProps;i++)
if(lineproplist[i].BdryFormat==3) SDIflag=TRUE;
if(SDIflag==TRUE){
// there is an SDI boundary defined; check to see if it is in use
SDIflag=FALSE;
for(i=0;i<NumEls;i++)
for(j=0;j<3;j++)
if (lineproplist[meshele[i].e[j]].BdryFormat==3){
SDIflag=TRUE;
printf("Problem has SDI boundaries\n");
i=NumEls;
j=3;
}
}
if (SDIflag==TRUE)
{
V_sdi=(CComplex *) calloc(NumNodes+NumCircProps,sizeof(CComplex));
sdin=2;
}
else sdin=1;
// compute effective permeability for each block type;
Mu=(CComplex **)calloc(NumBlockProps,sizeof(CComplex *));
for(i=0;i<NumBlockProps;i++) Mu[i]=(CComplex *)calloc(2,sizeof(CComplex));
for(k=0;k<NumBlockProps;k++){
if (blockproplist[k].LamType==0){
Mu[k][0]=blockproplist[k].mu_x*exp(-I*blockproplist[k].Theta_hx*PI/180.);
Mu[k][1]=blockproplist[k].mu_y*exp(-I*blockproplist[k].Theta_hy*PI/180.);
if(blockproplist[k].Lam_d!=0){
if (blockproplist[k].Cduct != 0){
halflag=exp(-I*blockproplist[k].Theta_hx*PI/360.);
ds=sqrt(2./(0.4*PI*w*blockproplist[k].Cduct*blockproplist[k].mu_x));
K=halflag*deg45*blockproplist[k].Lam_d*0.001/(2.*ds);
Mu[k][0]=((Mu[k][0]*tanh(K))/K)*blockproplist[k].LamFill +
(1.-blockproplist[k].LamFill);
halflag=exp(-I*blockproplist[k].Theta_hy*PI/360.);
ds=sqrt(2./(0.4*PI*w*blockproplist[k].Cduct*blockproplist[k].mu_y));
K=halflag*deg45*blockproplist[k].Lam_d*0.001/(2.*ds);
Mu[k][1]=((Mu[k][1]*tanh(K))/K)*blockproplist[k].LamFill +
(1.-blockproplist[k].LamFill);
}
else{
Mu[k][0]=Mu[k][0]*blockproplist[k].LamFill +
(1.- blockproplist[k].LamFill);
Mu[k][1]=Mu[k][1]*blockproplist[k].LamFill +
(1. - blockproplist[k].LamFill);
}
}
}
else{
Mu[k][0]=1;
Mu[k][1]=1;
}
}
do{
for(sdi_iter=0; sdi_iter<sdin; sdi_iter++)
{
TheView->SetDlgItemText(IDC_FRAME1,"Matrix Construction");
TheView->m_prog1.SetPos(0);
pctr=0;
if (Iter>0) L.Wipe();
// build element matrices using the matrices derived in Allaire's book.
for(i=0;i<NumEls;i++){
// update ``building matrix'' progress bar...
j=(i*20)/NumEls+1;
if(j>pctr){
j=pctr*5; if (j>100) j=100;
TheView->m_prog1.SetPos(j);
pctr++;
}
// zero out Me, be;
for(j=0;j<3;j++){
for(k=0;k<3;k++)
{
Me[j][k]=0;
Mx[j][k]=0;
My[j][k]=0;
Mn[j][k]=0;
}
be[j]=0;
}
// Determine shape parameters.
// l == element side lengths;
// p corresponds to the `b' parameter in Allaire
// q corresponds to the `c' parameter in Allaire
El=&meshele[i];
for(k=0;k<3;k++){
n[k]=El->p[k];
rn[k]=meshnode[n[k]].x;
}
p[0]=meshnode[n[1]].y - meshnode[n[2]].y;
p[1]=meshnode[n[2]].y - meshnode[n[0]].y;
p[2]=meshnode[n[0]].y - meshnode[n[1]].y;
q[0]=meshnode[n[2]].x - meshnode[n[1]].x;
q[1]=meshnode[n[0]].x - meshnode[n[2]].x;
q[2]=meshnode[n[1]].x - meshnode[n[0]].x;
g[0]=(meshnode[n[2]].x + meshnode[n[1]].x)/2.;
g[1]=(meshnode[n[0]].x + meshnode[n[2]].x)/2.;
g[2]=(meshnode[n[1]].x + meshnode[n[0]].x)/2.;
for(j=0,k=1;j<3;k++,j++){
if (k==3) k=0;
l[j]=sqrt( pow(meshnode[n[k]].x-meshnode[n[j]].x,2.) +
pow(meshnode[n[k]].y-meshnode[n[j]].y,2.) );
}
a=(p[0]*q[1]-p[1]*q[0])/2.;
R=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;
for(j=0,a_hat=0;j<3;j++) a_hat+=(rn[j]*rn[j]*p[j]/(4.*R));
vol=2.*R*a_hat;
for(j=0,flag=0;j<3;j++) if(rn[j]<1.e-06) flag++;
switch(flag)
{
case 2:
R_hat=R;
break;
case 1:
if(rn[0]<1.e-06){
if (fabs(rn[1]-rn[2])<1.e-06) R_hat=rn[2]/2.;
else R_hat=(rn[1] - rn[2])/(2.*log(rn[1]) - 2.*log(rn[2]));
}
if(rn[1]<1.e-06){
if (fabs(rn[2]-rn[0])<1.e-06) R_hat=rn[0]/2.;
else R_hat=(rn[2] - rn[0])/(2.*log(rn[2]) - 2.*log(rn[0]));
}
if(rn[2]<1.e-06){
if (fabs(rn[0]-rn[1])<1.e-06) R_hat=rn[1]/2.;
else R_hat=(rn[0] - rn[1])/(2.*log(rn[0]) - 2.*log(rn[1]));
}
break;
default:
if (fabs(q[0])<1.e-06)
R_hat=(q[1]*q[1])/(2.*(-q[1] + rn[0]*log(rn[0]/rn[2])));
else if (fabs(q[1])<1.e-06)
R_hat=(q[2]*q[2])/(2.*(-q[2] + rn[1]*log(rn[1]/rn[0])));
else if (fabs(q[2])<1.e-06)
R_hat=(q[0]*q[0])/(2.*(-q[0] + rn[2]*log(rn[2]/rn[1])));
else
R_hat=-(q[0]*q[1]*q[2])/
(2.*(q[0]*rn[0]*log(rn[0]) +
q[1]*rn[1]*log(rn[1]) +
q[2]*rn[2]*log(rn[2])));
break;
}
// Mr Contribution
// Derived from flux formulation with c0 + c1 r^2 + c2 z
// interpolation in the element.
K=(-1./(2.*a_hat*R));
for(j=0;j<3;j++)
for(k=j;k<3;k++)
Mx[j][k] += K*p[j]*rn[j]*p[k]*rn[k];
// need this loop to avoid singularities. This just puts something
// on the main diagonal of nodes that are on the r=0 line.
// The program later sets these nodes to zero, but it's good to
// for scaling reasons to grab entries from the neighboring diagonals
// rather than just setting these entries to 1 or something....
for(j=0;j<3;j++)
if (rn[j]<1.e-06) Mx[j][j]+=Mx[0][0]+Mx[1][1]+Mx[2][2];
// Mz Contribution;
// Derived from flux formulation with c0 + c1 r^2 + c2 z
// interpolation in the element.
K=(-1./(2.*a_hat*R_hat));
for(j=0;j<3;j++)
for(k=j;k<3;k++)
My[j][k] += K*(q[j]*rn[j])*(q[k]*rn[k])*
(g[j]/R)*(g[k]/R);
// Fill out rest of entries of Mx and My;
Mx[1][0]=Mx[0][1]; Mx[2][0]=Mx[0][2]; Mx[2][1]=Mx[1][2];
My[1][0]=My[0][1]; My[2][0]=My[0][2]; My[2][1]=My[1][2];
// contribution from eddy currents;
// induced current interpolated as constant (avg. of nodal values)
// over the entire element;
K = -I*R*a*w*blockproplist[meshele[i].blk].Cduct*c/6.;
// radially laminated blocks appear to have no conductivity;
// eddy currents are accounted for in these elements by their
// frequency-dependent permeability.
if((blockproplist[El->blk].LamType==0) &&
(blockproplist[El->blk].Lam_d>0)) K=0;
// if this element is part of a wound coil,
// it should have a zero "bulk" conductivity...
if(abs(labellist[El->lbl].Turns)!=1) K=0;
for(j=0;j<3;j++)
for(k=0;k<3;k++)
Me[j][k]+=K*4./3.;
// contributions to Me, be from derivative boundary conditions;
for(j=0;j<3;j++){
k=j+1; if(k==3) k=0;
r=(meshnode[n[j]].x+meshnode[n[k]].x)/2.;
if (El->e[j] >= 0){
if (lineproplist[El->e[j]].BdryFormat==2)
{
// conversion factor is 10^(-4) (I think...)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -