Startpage >> Main >> Emg

Emg

Se quiere estudiar una máquina estática electromagnética alimentada con corriente continua.

* Documentación explicativa pdf

* Código MATLAB MaquinaEstatica.m

campocampo mu

Gracias a Juan Álvaro Fuentes, del dpto. de Ingeniería Eléctrica UPCT juanalvaro@upct.es

Código FreeFem++ MaquinaEstaticaPresentacion.edp

// FREEFEM++
// Notas:
// * Las variables no pueden contener el caracter _
// * La vida de una variable es el bloque actual: {} excepto la variable fespace que es global. 
// Las variables locales son destruidas al terminar el bloque
// * Mucho cuidado con divisiones del tipo 1/2 que las realiza como números enteros (igual que C)  
//  y por lo tanto el resultado es cero y no 0.5 que es lo deseado!!!!!!

// Maquina estática con entrehierro largo dimensión perpendicular al plano es de 3 metros.

//////////////////////////////////////////////////////////////////////////////////////////////////
// Parámetros
//////////////////////////////////////////////////////////////////////////////////////////////////
// de los gráficos
bool esperar = true;

// del mallado
int nbobina   = 20;    // Número de divisiones de los contornos de las bobinas
int naire     = 40;    // Número de divisiones de los contornos del hueco interior
int ncontorno = 100;   // Número de divisiones de los contornos exteriores

// de la nolinealidad
int maxIter     = 20;
real tolerancia = 1e-4;
real residuo    = 0;
real avance     = 0.25;

// de la máquina
real longitudMaquina = 3;
real ladoBobina = 0.3;

// Función que para la ejecución del script hasta que no pulses la tecla de retorno
func string pausa()
{
	cout << "Pulsa una tecla para continuar" << endl; 
	string s;
	getline(cin,s);
	return s;
}

//////////////////////////////////////////////////////////////////////////////////////////////////
// Dominio y malla
//////////////////////////////////////////////////////////////////////////////////////////////////

///////////////// BOBINAS
// Contorno de la Bobina Inferior Derecha (BID)
// l -> límite, izq y der, b -> bobina, i -> inferior, d-> derecha
real x13 = 0.35;
real y13 = -0.15;
real lizqbid = x13;
real lderbid = x13+ladoBobina;
real lsupbid = y13;
real linfbid = y13-ladoBobina;
border CDerBID (t = linfbid, lsupbid) {x = lderbid; y = t;};
border CSupBID (t = lderbid, lizqbid) {x = t;       y = lsupbid;};
border CIzqBID (t = lsupbid, linfbid) {x = lizqbid; y = t;};
border CInfBID (t = lizqbid, lderbid) {x = t;       y = linfbid;};

// Contorno de la Bobina Superior Derecha (BSD)
real x14 = 0.35;
real y14 = 0.15;
real lizqbsd = x14;
real lderbsd = x14+ladoBobina;
real lsupbsd = y14+ladoBobina;
real linfbsd = y14;
border CDerBSD (t = linfbsd, lsupbsd) {x = lderbsd; y = t;};
border CSupBSD (t = lderbsd, lizqbsd) {x = t;       y = lsupbsd;};
border CIzqBSD (t = lsupbsd, linfbsd) {x = lizqbsd; y = t;};
border CInfBSD (t = lizqbsd, lderbsd) {x = t;       y = linfbsd;};

// Contorno de la Bobina Superior Izquierda (BSI)
real x15 = -0.35;
real y15 = 0.15;
real lizqbsi = x15-ladoBobina;
real lderbsi = x15;
real lsupbsi = y15+ladoBobina;
real linfbsi = y15;
border CDerBSI (t = linfbsi, lsupbsi) {x = lderbsi; y = t;};
border CSupBSI (t = lderbsi, lizqbsi) {x = t;       y = lsupbsi;};
border CIzqBSI (t = lsupbsi, linfbsi) {x = lizqbsi; y = t;};
border CInfBSI (t = lizqbsi, lderbsi) {x = t;       y = linfbsi;};

// Contorno de la Bobina Inferior Izquierda (BII)
real x16 = -0.35;
real y16 = -0.15;
real lizqbii = x16-ladoBobina;
real lderbii = x16;
real lsupbii = y16;
real linfbii = y16-ladoBobina;
border CDerBII (t = linfbii, lsupbii) {x = lderbii; y = t;};
border CSupBII (t = lderbii, lizqbii) {x = t;       y = lsupbii;};
border CIzqBII (t = lsupbii, linfbii) {x = lizqbii; y = t;};
border CInfBII (t = lizqbii, lderbii) {x = t;       y = linfbii;};

///////////////// VENTANA DEL CIRCUITO MAGNÉTICO

// Defino el hueco del interior por tres rectángulos => el izquierdo (Izq), el centro (Cen) y el
// derecho (Der) 
real x5 = 0.7;
real y5 = -0.5;
real x6 = 0.7;
real y6 = 0.5;
real x7 = -0.7;
real y7 = 0.5;
real x8 = -0.7;
real y8 = -0.5;
real x9 = 0.3;
real y9 = -0.1;
real x10 = 0.3;
real y10 = 0.1;
real x11 = -0.3;
real y11 = 0.1;
real x12 = -0.3;
real y12 = -0.1;
// Rectangulo de la derecha
real xIzquierdaDer = x9; 
real xDerechaDer = x5;
real ySuperiorDer = y6;
real yInferiorDer = y5;
// Rectangulo del centro
real xIzquierdaCen = x11; 
real xDerechaCen = x10;
real ySuperiorCen = y10;
real yInferiorCen = y12;
// Rectangulo de la izquierda
real xIzquierdaIzq = x7; 
real xDerechaIzq = x11;
real ySuperiorIzq = y7;
real yInferiorIzq = y8;
// Generamos la ventana del circuito magnético en sentido antihorario
// R -> Rectángulo, D -> Derecha, C -> Centro, I -> Izquierdo
border CDerRD  (t = yInferiorDer, ySuperiorDer)  {x = xDerechaDer;   y = t;};
border CSupRD  (t = xDerechaDer,  xIzquierdaDer) {x = t;             y = ySuperiorDer;};
border CIzqRD  (t = ySuperiorDer, ySuperiorCen)  {x = xIzquierdaDer; y = t;};
border CSupRC  (t = xDerechaCen,  xIzquierdaCen) {x = t;             y = ySuperiorCen;};
border CDerRI  (t = ySuperiorCen, ySuperiorIzq)  {x = xDerechaIzq;   y = t;};
border CSupRI  (t = xDerechaIzq,  xIzquierdaIzq) {x = t;             y = ySuperiorIzq;};
border CIzqRI  (t = ySuperiorIzq, yInferiorIzq)  {x = xIzquierdaIzq; y = t;};
border CInfRI  (t = xIzquierdaIzq, xDerechaIzq)  {x = t;             y = yInferiorIzq;};
border CDerRI2 (t = yInferiorIzq, yInferiorCen)  {x = xDerechaIzq;   y = t;};
border CInfRC  (t = xIzquierdaCen, xDerechaCen)  {x = t;             y = yInferiorCen;};
border CIzqRD2 (t = yInferiorCen, yInferiorDer)  {x = xIzquierdaDer; y = t;};
border CInfRD  (t = xIzquierdaDer, xDerechaDer)  {x = t;             y = yInferiorDer;};

///////////////// CONTORNO DEL DOMINIO DEL PROBLEMA

// DP = Dominio del Problema
real x1 = 1.1;
real yy1 = -1.1;  // Oculta un identificador global
real x2 = 1.1;
real y2 = 1.1;
real x3 = -1.1;
real y3 = 1.1;
real x4 = -1.1;
real y4 = -1.1;
real xInfDP = x1;
real yInfDP = yy1;
real xSupDP = x2;
real ySupDP = y2;
real xSupIP = x3;
real ySupIP = y3;
real xInfIP = x4;
real yInfIP = y4;

border CDerP(t = yInfDP, ySupDP) {x = xInfDP;   y = t;};
border CSupP(t = xSupDP, xSupIP) {x = t;        y = ySupDP;};
border CIzqP(t = ySupIP, yInfIP) {x = xSupIP;   y = t;};
border CInfP(t = xInfIP, xInfDP) {x = t;        y = yInfDP;};

//////////////////////////////////////////////////////////////////////////////////////////////////
// Malla del problema
//////////////////////////////////////////////////////////////////////////////////////////////////
mesh Th = buildmesh(CDerP(ncontorno) + CSupP(ncontorno) + CIzqP(ncontorno) + CInfP(ncontorno)
                    + CDerRD(2*naire) + CSupRD(naire) + CIzqRD(naire)
                    + CSupRC(naire)
                    + CDerRI(naire) + CSupRI(naire) + CIzqRI(2*naire) + CInfRI(naire) + CDerRI2(naire)
                    + CInfRC(naire)
                    + CIzqRD2(naire) + CInfRD(naire)
                    + CDerBID(nbobina) + CSupBID(nbobina) + CIzqBID(nbobina) + CInfBID(nbobina)
                    + CDerBSD(nbobina) + CSupBSD(nbobina) + CIzqBSD(nbobina) + CInfBSD(nbobina)
                    + CDerBSI(nbobina) + CSupBSI(nbobina) + CIzqBSI(nbobina) + CInfBSI(nbobina)
                    + CDerBII(nbobina) + CSupBII(nbobina) + CIzqBII(nbobina) + CInfBII(nbobina));

plot (Th, wait = esperar, cmm = "malla", ps="me_malla.png");
//////////////////////////////////////////////////////////////////////////////////////////////////
// Espacios de elementos finitos
//////////////////////////////////////////////////////////////////////////////////////////////////
fespace Vh(Th,P2); // Espacio para calcular el potencial vectorial magnético A
fespace Bh(Th,P1); // Espacio para calcular B
fespace Ch(Th,P0); // Espacio para definir ctes

// funciones de Vh para cálculo del potencial vectorial A
Vh u, v, uold;

// funciones de Bh para el cálculo de B
Bh Bx, By, mu, muold, murelativa, Bmodulo;

//////////////////////////////////////////////////////////////////////////////////////////////////
// Obtención de las regiones para definir las densidades de corriente y la permeabilidad magnética
//////////////////////////////////////////////////////////////////////////////////////////////////
Ch reg = region;
plot(reg, fill = 1, wait = esperar, value = 1, cmm = "regiones", ps="me_regiones.png");
int bobinaSuperiorDerecha = reg(lizqbsd+1e-3,  linfbsd+1e-3);
int bobinaSuperiorIzquierda = reg(-lderbsd+1e-3, linfbsd+1e-3);
int bobinaInferiorDerecha = reg(lizqbsd+1e-3, -lsupbsd+1e-3);
int bobinaInferiorIzquierda = reg(-lderbsd+1e-3,-lsupbsd+1e-3);
int aire = reg(0, 0);
int hierro = reg(-0.8,-0.8);

cout << "Region Bobina Superior Derecha: "   << bobinaSuperiorDerecha << endl;
cout << "Region Bobina Superior Izquierda: " << bobinaSuperiorIzquierda << endl;
cout << "Region Bobina Inferior Derecha: "   << bobinaInferiorDerecha << endl;
cout << "Region Bobina Inferior Izquierda: " << bobinaInferiorIzquierda << endl;
cout << "Region Aire: "                      << aire << endl;
cout << "Region Hierro: "                    << hierro << endl;

//////////////////////////////////////////////////////////////////////////////////////////////////
// Densidades de corriente
//////////////////////////////////////////////////////////////////////////////////////////////////
real densidadCorriente = 1e6;

func j = (+densidadCorriente*(reg == bobinaSuperiorDerecha)
          -densidadCorriente*(reg == bobinaSuperiorIzquierda)
          +densidadCorriente*(reg == bobinaInferiorDerecha)
          -densidadCorriente*(reg == bobinaInferiorIzquierda));
Ch jh = j;
plot(jh, fill = 1, wait = esperar, value = 1, cmm = "densidades de corriente", ps="me_densidades_corriente.png");

//////////////////////////////////////////////////////////////////////////////////////////////////
// Permeabilidades
//////////////////////////////////////////////////////////////////////////////////////////////////
// Permeabilidad del vacío
real mu0 = 4*pi*1e-7;
// Permeabilidad inicial del hierro
real murelhierroini = 1000;
// Permeabilidad inicial de todo el dominio
// Defino de forma conjunta la región de Cobre y Aire ya que tienen la misma permeabilidad
func regionCobreAire = (+1*(reg == bobinaSuperiorDerecha)
                        +1*(reg == bobinaSuperiorIzquierda)
                        +1*(reg == bobinaInferiorDerecha)
                        +1*(reg == bobinaInferiorIzquierda)
                        +1*(reg == aire));          
// Defino la región de hierro
func regionHierro = (reg == hierro);
func muini = mu0*regionCobreAire+murelhierroini*mu0*regionHierro;

///////////////// DISTINTOS AJUSTES DE PERMEABILIDAD DEL MATERIAL

// // 1) Problema lineal
// func mu = muini;
// // 2) Ajuste original de César
// Tiene el inconveniente de que para B = 0 da permeabilidades nulas de ahí que para evitar problemas sumemos 1e-10
// func mu = mu0*regionCobreAire + (abs(sqrt(Bx^2+By^2)/(0.4367*(exp(5.42*(sqrt(Bx^2+By^2)))-0.999)))+1e-10)*regionHierro;
// // 3) Función de César y el algoritmo sqp.
// real a = 2.81765548102844;
// real b = 3.70241760042811;
// real c = -16.55391737676212;
// func mu = mu0*regionCobreAire + (abs(Bmodulo/(a*(exp(b*(Bmodulo))-c)))+1e-10)*regionHierro; 
// 4) Función de maplesoft ajustando el valor de BmyMax ajusta mejor para valores altos (podemos
// aproximar mejor los valores altos incrementando B_mu_max en el ajuste de maquina_estatica.m)
real mui = 1210.0;
real BmyMax = 1.0;
real ca = 24630.0008070838;
real cb = 2.80738309289622;
real n = 8.09712848132294;
func BN = sqrt(Bx^2+By^2)/BmyMax;
func muRelativa = 1+(mui-1+ca*BN)/(1+cb*BN+BN^n);
func muNolineal = mu0*regionCobreAire + mu0*muRelativa*regionHierro;

//////////////////////////////////////////////////////////////////////////////////////////////////
// Resolvemos el problema antes de entrar al bucle
//////////////////////////////////////////////////////////////////////////////////////////////////
mu = muini;
macro Grad2(u) [dx(u),dy(u)] // EOM, gradiente
problem magnetostatica(u,v)=int2d(Th)((1/mu)*(Grad2(v)'*Grad2(u)))
  -int2d(Th)(j*v)
  +on(CDerP, CSupP, CIzqP, CInfP, u=0);
magnetostatica;

// Calculo de las componentes de B antes de entrar al bucle
Bx=dy(u);
By=-dx(u);

//////////////////////////////////////////////////////////////////////////////////////////////////
// Bucle para el cálculo de la pertimividad no-lineal
//////////////////////////////////////////////////////////////////////////////////////////////////
for (int i=1; i<=maxIter; i++){
  // Guardamos los valores de la iteracion anterior
  uold = u;
  muold = mu;

  //Calculo de las componentes de B para el cálculo de las nuevas permeabilidades
  Bx=dy(u);
  By=-dx(u);

  // Avanzamos poco a poco para posibilitar la convergencia
  mu = avance*(muNolineal) + (1-avance)*muold;
  // Volvemos a resolver el problema pero con las permeabilidades recalculadas
  magnetostatica;

  // Representación de los valores de mu relativos
  murelativa=mu/mu0;
  //if ((i/10.0-floor(i/10)) == 0)
  //  plot(murelativa,value=true, wait = esperar, cmm = "mu relativa i= "+i);

  cout << "Valor: " << (i/10.0-floor(i/10))*10 << endl;
  // Hemos convergido???
  residuo = int2d(Th)(abs(uold-u));
  cout << "Residuo: " << residuo << endl;
  if (residuo <= tolerancia) {
    cout << "Hemos convergido!!!" << endl;
    break;
  };
  if(i == maxIter){ // Nota hemos convergido
    cout << "Maxima iteracion alcanzada: " << i << endl;
  };
}

////////////////////////////////////////
// Representación gráfica de la solución
////////////////////////////////////////
// Potencial vector (componente z)
plot(u, value = true , wait = esperar, cmm ="Az (lineas de campo de B)", ps="me_Az.png");

// Densidad de flujo
plot([Bx,By], value = true, coef = -1.6, wait = esperar, cmm = "Vector de densidad de flujo B", ps="me_vector_B.png");
Bmodulo = sqrt(Bx^2+By^2);
plot(Bmodulo, value = true , wait = esperar, cmm ="Densidad de flujo B", ps="me_modulo_B.png");

// Permeabilidad relativa
plot(murelativa,value=true, wait = esperar, cmm = "mu relativa", ps = "me_mu relativa.png");

// Intensidad de campo magnético
Bh Hx, Hy, auxX, auxY;
Hx = Bx/mu;
Hy = By/mu;

auxX = Hx/1e6;
auxY = Hy/1e6;
plot([auxX,auxY], value = true, coef = -1.6, wait = esperar, cmm = "Vector H MegaA/m", ps = "me_vector_H.png");

// Dibujamos el campo vectorial de Magnetización
Bh Mx, My;
Mx = Bx*(murelativa-1)/(murelativa*mu0);
My = By*(murelativa-1)/(murelativa*mu0);
auxX = Mx/1e6;
auxY = My/1e6;
plot([auxX,auxY],value=true, coef = -1.6, wait = esperar, 
cmm = "Vector de Magnetizacion (MegaA/m) en y=0", ps = "me_vector_M.png");

// Calculos varios
cout << "Valor de Hx(0,0) (MA/m): " << Bx(0,0)/mu(0,0)/1e6 << endl;
cout << "Valor de Hy(0,0) (MA/m): " << By(0,0)/mu(0,0)/1e6 << endl;

///////////////// Sección por el eje y
// Dibujamos valores de By, Hy en la recta y = 0
real numpuntos = 100.0;
real[int] valoresX(numpuntos),valoresBy(numpuntos),valoresHy(numpuntos),valoresMy(numpuntos),valoresAz(numpuntos);
for (int i=0; i<numpuntos; i++)
{
  valoresX[i] = xSupIP+(xSupDP-xSupIP)/(numpuntos-1)*i;
  valoresBy[i] = By(valoresX[i],0); // valores de By en el punto (valoresX[i], 0.);
	valoresHy[i] = valoresBy[i]/(murelativa(valoresX[i],0)*mu0*1e6); // MegaAmperios/m
  valoresMy[i] = valoresHy[i]*(murelativa(valoresX[i],0)-1);       // MegaAmperios/m
	valoresAz[i] = u(valoresX[i],0);
}
// cout << xx << yy << endl;
plot([valoresX,valoresBy], wait = esperar, cmm = "By (T) seccion eje simetria (y=0)", ps = "me_seccion_By.png");
plot([valoresX,valoresHy], wait = esperar, cmm = "Hy (MegaA/m) en y=0", ps = "me_seccion_Hy.png");
plot([valoresX,valoresMy], wait = esperar, cmm = "My (MegaA/m) en y=0", ps = "me_seccion_My.png");

// Dibujamos Az y By en la recta y = 0
plot([valoresX,valoresAz], [valoresX, valoresBy], wait = esperar, 
cmm = "Az y By en y=0", ps = "me_seccion_Az_By.png");

///////////////// Circulación por el eje de antisimetría y el contorno derecho del dominio
// Cálculo de la circulación. La ley de Ampère dice que la circulación de H es igual a la intensidad abrazada
// Vamos aplicar la ley de Ampère a una curva que pasa por el eje de 
// antisimetría y por los contornos exteriores del problema
// Definimos el contorno de la curva
border C2DerP(t = yInfDP, ySupDP)           {x = xInfDP;   y = t;};
border C2SupP(t = xSupDP, 0)                {x = t;        y = ySupDP;};
border EjeAntisimetria (t = ySupDP, yInfDP) {x = 0;        y = t;};
border C2InfP(t = 0, xInfDP)                {x = t;        y = yInfDP;};
// Creamos una malla cuyos límites contengan el contorno de la curva sobre la que queremos integrar
mesh Th2 = buildmesh(C2DerP(ncontorno) + C2SupP(ncontorno) + EjeAntisimetria(8*ncontorno) + C2InfP(ncontorno));
fespace Bh2(Th2,P1);
// Dado vector normal (N.x,N.y) hacia el exterior el perpendicular en sentido horario será -N.y,N.x y en sentido antihorario será N.y,-N.x
cout << "Circulacion H sobre el eje de Antisimetria: " << int1d(Th2,EjeAntisimetria) (-Hx*N.y + Hy*N.x) << endl;
cout << "Circulacion H sobre los dos horizontales: " << int1d(Th2,C2SupP,C2InfP) (-Hx*N.y + Hy*N.x) << endl;
cout << "Circulacion H sobre el vertinal del nucleo: " << int1d(Th2,C2DerP) (-Hx*N.y + Hy*N.x) << endl;
real circulacionH = int1d(Th2,C2DerP,C2SupP,EjeAntisimetria,C2InfP) (-Hx*N.y + Hy*N.x);
cout << "Circulacion H sobre el contorno cerrado que pasa por el eje de Antisimetria en sentido horario: " << circulacionH << endl;

///////////////// Intensidad abrazada por la curva formada por el eje de antisimetría y el contorno derecho del dominio
real intensidadAbrazada = int2d(Th2) (j);
cout << "Intensidad abrazada por el eje de antisimetria y el contorno exterior derecho: " << intensidadAbrazada << endl;

// Ahora lo mismo pero para la densidad de flujo
cout << "Circulacion B sobre el eje de Antisimetria: " << int1d(Th2,EjeAntisimetria) (-Bx*N.y + By*N.x) << endl;
cout << "Circulacion B sobre los dos horizontales: " << int1d(Th2,C2SupP,C2InfP) (-Bx*N.y + By*N.x) << endl;
cout << "Circulacion B sobre el vertinal del nucleo: " << int1d(Th2,C2DerP) (-Bx*N.y + By*N.x) << endl;
real circulacionB = int1d(Th2,C2DerP,C2SupP,EjeAntisimetria,C2InfP) (-Bx*N.y + By*N.x);
cout << "Circulacion B sobre el contorno cerrado que pasa por el eje de Antisimetria en sentido horario: " << circulacionB << endl;

// Permeabilidad relativa media
cout << "Permeabilidad relativa media sobre el contorno del eje de Antisimetria: "
 << circulacionB/(mu0*intensidadAbrazada) << endl;

// Ahora lo mismo pero para la magnetización
cout << "Circulacion M sobre el eje de Antisimetria: " << int1d(Th2,EjeAntisimetria) (-Mx*N.y + My*N.x) << endl;
cout << "Circulacion M sobre los dos horizontales: " << int1d(Th2,C2SupP,C2InfP) (-Mx*N.y + My*N.x) << endl;
cout << "Circulacion M sobre el vertinal del nucleo: " << int1d(Th2,C2DerP) (-Mx*N.y + My*N.x) << endl;
real circulacionM = int1d(Th2,C2DerP,C2SupP,EjeAntisimetria,C2InfP) (-Mx*N.y + My*N.x);
cout << "Circulacion M sobre el contorno cerrado que pasa por el eje de
          Antisimetria en sentido horario: " << circulacionM << endl;

// Relación entre las corriente de magnetización encerrada por el contorno y la corriente libre
cout << "Relacion intensidad de magnetizaction frente a la libre encerrada por 
                   el contorno: " << circulacionM/circulacionH << endl;

///////////////// Cálculo del flujo que atraviesa el plano de simetría
cout << "Flujo magnetico entre 0 < x < 0.3: " << longitudMaquina*(u(0,0)-u(0.3,0)) << endl;
cout << "Flujo magnetico entre 0.3 < x < 0.7: " << longitudMaquina*(u(0.3,0)-u(0.7,0)) << endl;
cout << "Flujo magnetico entre 0.7 < x < 1.1: " << longitudMaquina*(u(0.7,0)-u(xInfDP,0)) << endl;
cout << "Suma de flujos: " << longitudMaquina*((u(0,0)-u(xInfDP,0))) << endl;

///////////////// Cálculo de fuerzas
// Campo de fuerzas magnéticas por unidad de volumen que actúa sobre los bobinados
Bh densidadFx = -j*By;
Bh densidadFy = j*Bx;
plot([densidadFx,densidadFy], coef = -7. , value = true, wait = esperar, 
                   cmm = "Vector de densidad de fuerzas en bobina", ps = "me_vector_fuerza_bobina.png");

// Fuerza total sobre la bobina superior derecha
real FxBSD = longitudMaquina*int2d(Th, bobinaSuperiorDerecha) (densidadFx);   // Fuerza total sobre la bobina, componente x
real FyBSD = longitudMaquina*int2d(Th, bobinaSuperiorDerecha) (densidadFy);    // Fuerza total sobre la bobina, componente y
cout << "Fx bobina superior derecha: " << FxBSD << endl;
cout << "Fy bobina superior derecha: " << FyBSD << endl;

// Punto de aplicación de la fuerza
// El par total será Suma(x*Fy)-Suma(y*Fx) =
//    x_aplicado*Ftotaly-y_aplicado*Ftotalx => x_aplicado 
//     = Suma(x*Fy)/Ftotaly; y_aplicado = Suma(y*Fx)/Ftotalx
real T1BSD = longitudMaquina*int2d(Th, bobinaSuperiorDerecha) (y*densidadFx);
real T2BSD = longitudMaquina*int2d(Th, bobinaSuperiorDerecha) (x*densidadFy);
cout << "T1 bobina superior derecha: " << T1BSD << endl;
cout << "T2 bobina superior derecha: " << T2BSD << endl;
cout << "(x,y) donde Ftotal es aplicada a la bobina superior derecha: (" <<   T1BSD/FxBSD << "," << T2BSD/FyBSD << ")" << endl;

///////////////// Energías magnéticas
Bh densidadEnergia = (Bx*Hx+By*Hy)/2;
plot(densidadEnergia, value = true , wait = esperar, cmm ="Densidad de energia", ps = "me_densidad_energia.png");

real energiaAire = longitudMaquina*int2d(Th, aire) (densidadEnergia);
real energiaHierro = longitudMaquina*int2d(Th, hierro) (densidadEnergia);
real energiaCobre = longitudMaquina*int2d(Th, bobinaSuperiorDerecha,
          bobinaSuperiorIzquierda,bobinaInferiorDerecha,bobinaInferiorIzquierda) (densidadEnergia);
real energiaTotal = energiaAire+energiaHierro+energiaCobre;

cout << "Energia magnetica del aire: " << energiaAire << endl;
cout << "Energia magnetica del hierro: " << energiaHierro << endl;
cout << "Energia magnetica del cobre: " << energiaCobre << endl;
cout << "Energia total: " << energiaTotal << endl;

Volver

Page last modified on May 21, 2014, at 07:38 PM