Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Mtodos Numricos para la resolucin de Ecuaciones

Diferenciales en Derivadas Parciales en Electromagnetismo:
Magnetoesttica

Gonzalo Casaravilla
Universidad de la Repblica - Uruguay

Resumen: En el presente trabajo se analiza el planteo general de las

ecuaciones de Maxwell que intervienen en magnetoesttica y en particular
en aplicaciones 2D. El objeto final es tomar conocimiento del estado del
arte en cuanto a las herramientas de clculo de los esfuerzos estticos en
los bobinados de un transformador de potencia elctrica de mltiples
arroyamientos. Dentro de este marco se analiza el papel que juegan las
Diferencias Finitas, la forma de plantear el sistema de ecuaciones
implcitas que resuelven la ecuacin elptica de Poisson y las condiciones
de frontera de Neumann que aparecen en magnetoesttica. Asimismo se
analiza la estabilidad de la solucin y si existe a priori una forma de
seleccionar la grilla de clculo. Como resultado de este anlisis se intenta
contestar en primera instancia la pregunta qu es lo que se usa en la
actualidad?. Como ejemplo particular se analiza una herramienta de
clculo acadmica/comercial que desde el punto de vista prctico termina
contestando satisfactoriamente la motivacin y objeto de este trabajo.

1 Introduccin

Las ecuaciones de Maxwell [1][2] resumen las leyes del electromagnetismo clsico.
Bsicamente son ecuaciones en derivadas parciales (PDE) y en su forma mas general
involucra las tres dimensiones y el tiempo. Desde sus orgenes han surgido diversos
mtodos de anlisis y clculo basados en formulaciones matemticas de elevada
complejidad, ya sea en el dominio del tiempo o la frecuencia.

Sin embargo, con el advenimiento de las potentes computadoras y sus posibilidades de

clculo, la resolucin con mtodos numricos de las ecuaciones de Maxwell han
evolucionado explosivamente.

Como mtodos por excelencia se pueden individualizar el mtodo de los Elementos

Finitos (FE) y las Diferencias Finitas (FD). En los casos planos y en particular
estacionarios, la FE han dominado el mercado. En los casos tridimensionales, ya sea en
el dominio de las frecuencias o en el dominio del tiempo, las FD estn en franca
competencia con los FE.

En particular se puede decir que el mtodo de las FD en el dominio del tiempo

(FDTD)[3][4][5] [6] es ya un estndar. Basan su suceso en la simplicidad de
formulacin. Bsicamente son aproximaciones de segundo orden de las ecuaciones en
derivadas parciales en grillas desfasadas en el tiempo y en el espacio para las 2 variables
independiente del sistema electromarntico (el campo magntico B y el campo elctrico
E). El efecto de este deslocamiento temporal y espacial permite sistemas de ecuaciones

Gonzalo Casaravilla: [email protected] pag. 1 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

explcitos que a partir de condiciones de borde conocidas permite desarrollar los

clculos en forma secuencial.
El estado del arte actual indica que los esfuerzos se estn volcando hacia un mejor
aprovechamiento de la memoria durante el clculo y la adopcin de estrategias de
procesamiento paralelo adecuadas dada la multigrilla adoptada [4] [13].

Si nos atenemos al objeto inicial de este trabajo y se analizan, por ejemplo, los mtodos
de clculo de los esfuerzos estticos dentro de transformadores, encontramos cuatro
mtodos de clculo. El mtodo de Rogowski , el mtodo de Rabins (Serie Simple de
Fourier cuyos coeficientes son funciones de Bessel y Struve), el Mtodo de Roth (Serie
Doble de Fourier), el mtodo de las Imgenes y finalmente el ya mencionado FE [7][8].

Queda por tanto planteada las duda de si las FD son adecuadas o no para resolver
problemas magnetoestticos (o sus anlogos electroestticos)

Respecto a las PED que aparecen en Electromagnetismo se puede decir que son
generalmente ecuaciones Elpticas que se reducen a los tres tipos bsicos de Laplace,
Poisson y Helmholtz [9][10][11] [12][13].

En cuanto a las condiciones de borde, existen casos de los tres tipos Neumann, Dirichlet
y, Mixtas [9][10][11] [12][13]. En temas magnticos generalmente hay condiciones de
borde en las que el campo magntico es perpendicular a una determinada frontera, por
lo que la condicin de borde termina siendo que alguna derivada parcial se anule
(Neumann). En cambio para casos en los que intervienen temas Elctricos (potenciales
elctricos), las condiciones terminan siendo imponer un valor en la frontera (Dirichlet).
Finalmente en problemas mixtos en donde queremos considerar campos magnticos y
elctricos aparecen naturalmente condiciones de borde mixtas.

2 Magnetoesttica

Las ecuaciones de Maxwell que intervienen en temas magnetoestticos son las

ecuaciones de Ampere y Gauss

D
HIERRO BT: Baja tensin
H = J+ Ampere AIRE AT: Alta tensin
t

B = 0 Gauss Bobinado BT 1

donde Bobinado BT 2

H es la induccin magntica Bobinado AT 1

B es el campo magntico
D es el desplazamiento elctrico Bobinado AT 2
J es la densidad de corriente elctrica

Figura 1: Corte de un transformador de varios bobinados
El vector H es el ROTOR de H

El escalar B = 0 es la DIVERGENCIA del campo magntico B

Gonzalo Casaravilla: [email protected] pag. 2 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

D
Si no hay desplazamiento de cargas libres, = 0 , por lo que la ecuacin de Ampere
t

queda H = J

Por otra parte, de la ecuacin de Gauss se tiene 1 que existe A tal que

B = A

Donde A se define como el potencial vector de B

Por otra parte en el aire se tiene que B = 0 J por lo tanto

( A) = 0 J

2.1.1 Caso particular a estudiar

Desarrollando esta expresin para problemas con simetra plana (2D), como por
ejemplo transformadores del tipo Core Type y Shell Type cuyo corte podra ser el
que nuestra la Figura 1, se tiene que
v r r
i j k
A( x, y ) r A( x, y ) r
( A) = = i j
x y z y x
0 0 A( x, y )

v r r
i j k
2 A( x, y ) 2 A( x, y ) r r
( A) = = + k = J ( x , y ) k
x y z x y 2
2 0

A( x, y ) A( x, y )
0
y x

Por lo que se obtiene as la Ecuacin de Laplace-Poisson en coordenadas cartesianas:

2 A( x, y ) 2 A( x, y )
+ = 0 J ( x, y )
x 2 y 2

que es una Ecuacin en Derivadas Parciales Elptica.

La densidad de corriente elctrica J es un dato del problema y se determina en funcin

de las corrientes que recorren los bobinados del transformador. Si J es nulo, la ecuacin

1
Se puede demostrar que ( A) = 0 A

Gonzalo Casaravilla: [email protected] pag. 3 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

se denomina de Laplace y para el caso de que J sea no nulo, la ecuacin se denomina

como de Poisson.

Luego de calculado A, el campo magntico B se determinar de la frmula ya vista

A( x, y ) r A( x, y ) r
B = A = i j
y x

por lo que si tenemos en cuenta que el GRADIENTE de A es

A( x, y ) r A( x, y ) r A( x, y ) A( x, y ) r A( x, y ) r
A = i+ j+ k= i+ j
x y z x y
resulta para este caso particular que

0 1
B= ( A)
1 0

2.1.2 Condiciones de borde para nuestro caso particular

Para el hierro, se supone que =, por lo cual se puede demostrar [1][2] a partir de las
propias ecuaciones de Gauss y Ampere que el campo magntico B del lado del aire en
las interfaces hierro-aire, es perpendicular a las mismas, hecho que debe ser tenido en
cuenta como condicin de borde en la resolucin de la ecuacin diferencial.

Por ejemplo, en el problema 2D planteado tenemos que

r A( x, y ) r A( x, y ) r
B = A = i j
y x

por lo que, por ejemplo, la condicin de borde en una frontera vertical ser que la
componente vertical de campo B sea nula, resultando

A( x, y )
By = =0
x

lo cual resulta en una condicin de borde de Neumann.

3 Diferencias finitas aplicadas a la magnetoesttica

Como ya se ha presentado, estudiaremos el caso particular de la ecuacin de Poisson

con condiciones de borde de Neumann en una grilla rectangular. Otras condiciones de
borde como ser las de Dirichlet son de ms fcil resolucin y no son las que
normalmente aparecen en magnetoesttica.

Como se ver, tendremos que discretizar en la grilla (n x m) de tres formas diferentes

q Dentro de la grilla. Tendremos (n-2)*(m-2) ecuaciones

Gonzalo Casaravilla: [email protected] pag. 4 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

q En los bordes sin contar los vrtices. Tendremos 2*(n-2)+2*(m-2) ecuaciones

q En los vrtices. Tendremos 4 ecuaciones

Por lo tanto tendremos (m-2)*(m-2)+2*(n-2)+2*(m-2)+4=nm ecuaciones

3.1 Discretizacin dentro de la grilla

j=m

y j

j=1
i=1 i=n
i
x

Figura2. Discretizacin dentro de la grilla

Tomando como clula de clculo la mostrada en la figura, la aproximacin de las

derivadas parciales de segundo orden resulta, por ejemplo para x

2 A Ai 1, j 2 Ai , j + Ai +1, j
= + O(x 2 )
x i , j
2
x 2

Sustituyendo esta expresin en la ecuacin de Poisson, resulta

Ai 1, j 2 Ai , j + Ai +1, j Ai , j 1 2 Ai , j + Ai , j +1
+ = 0 J ( x, y )
x 2
y 2

Esta ecuacin como ya se ha establecido se puede aplicar en los [(n-2)*(m-2)] puntos

interiores de la grilla resultando una serie de ecuaciones que vinculan las incgnitas A i,j.
Por tanto se tiene un sistema de ecuaciones implcito (resultado de la forma en que se
construyen las clulas de clculo de diferencias finitas).

Si definimos el vector S como una columna de dimensin mn con las incgnitas Aij
ordenados por las filas de la grilla tal como sigue

S1 A1, j
. .

S=S con S j = Ai , j

. .
S m An, j

Por tanto podremos ir formando la matriz Q (mn x mn) que determina el sistema de
ecuaciones a resolver tal que QS = G siendo G el trmino independiente de la
ecuacin de Poisson (-oJ(x,y)).

Gonzalo Casaravilla: [email protected] pag. 5 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Por ejemplo, para cada punto interior de la grilla (i,j), se aplica la discretizacin vista.
Tal cmo fue definido S, resulta que para cada punto (i,j) queda definida la fila k=i+n(j-
1) en la que se calcularn los trminos de la matriz Q.

Para simplificar la notacin definiremos

1 1
= = = -2( + )
x 2 y 2

Si observamos la discretizacin se puede inferir que, por ejemplo, tendremos la

incgnita A i-1,j multiplicada por , por lo que en la posicin h=i-1+n(j-1) de la fila k de
la matriz Q tendremos dicho valor . Luego la incgnita A i,j que corresponde a la
coordenada (k,k) de la matriz Q est multiplicada por . En el caso de Ai+1,j estar
multiplicada por y corresponde a la columna h=i+1+n(j-1). Paralelamente, el
multiplicador de A i,j-1 ser y estar en la columna h=i+n(j-1-1)=k-n y el multiplicador
de Ai,j+1 ser nuevamente y estar en la columna h=i+n(j+1-1)=k+n.

En resumen, para cada elemento interior de la grilla (i,j) tendremos que se definen
algunos elementos de la fila k=i+n(j-1) de la matriz Q, siendo los dems nulos.

k-n k k+n
[1..n]...[1....i....n ][1....i - 1 i i + 1.....n ][1....i....n ]...[1..n]
A 1,1


k 0............0 0.......0 0.........0 0.............0 Ai , j = 0 J (i, j )


An ,m

Al recorrer todos los elementos de la grilla, se ve como se va delineando una matriz Q

tridiagonal.

3.2 Discretizacin de las condiciones de borde

Analicemos el caso de la frontera inferior tal como muestra la Figura 3.

En este caso podemos hacer una aproximacin de primer orden del valor de A i,2 de la
siguiente forma

A 2 A y 2
Ai ,2 = Ai ,1 + .y + 2 . + O(y 3 )
y i ,1 y i ,1 2
123
=0

A
2
(Ai,2 Ai,1 ) ( )
= 2 + O y
y 2 i ,1 y 2

Luego, la discretizacin segn x es la misma ya vista, por lo que la ecuacin de Poisson

discretizada aproximada queda

Gonzalo Casaravilla: [email protected] pag. 6 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Ai 1,1 2 Ai ,1 + Ai +1,1 (A Ai ,1 )
+2 = 0 J ( x, y ) = 0 J (i,1)
i, 2

x 2 y 2

j=m

j=1
j=1
i=1 i=n

x i
Figura 3. Condiciones en la frontera j=1

Por tanto podemos determinar para los (n-2) puntos en consideracin perteneciente a
dicho borde, las filas de la matriz Q que quedan determinadas. Como j=1, las filas que
estaremos determinando son las filas k=2 .. (n-2) del vector S1.

k k+n
[1..........i - 1 i i + 1.....n ] [1....i....n ]........[1..n]
A1,1


k 0......0 0.........0 2 0.............0 Ai ,1 = 0 J (i,1)


An ,1

En el Apndice 1 se puede ver como quedan y donde contribuyen en la matriz Q las

otras tres condiciones de borde.

q La frontera j=m a las filas interiores de Sm

q La frontera i=1 a las primeras filas de S 2, S3..., Sm-1.
q La frontera i=n a las ltimas filas de S 2, S3..., Sm-1.

3.3 Discretizacin de los vrtices

Tomemos por ejemplo el vrtice 1,1 tal como muestra la Figura 4.

En este caso se puede hacer una aproximacin de primer orden en x e y obtenindose

como resultado

(A A1,1 ) (A A1,1 )
2 2 ,1
+2 1, 2
= 0 J ( x, y ) = 0 J (1,1)
x 2
y 2

que contribuye a la primer fila de Q

Gonzalo Casaravilla: [email protected] pag. 7 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

j=m

j=1
j=1
i=1 i=n
i=1
x
Figura 4. Discretizacin del vrtice 1,1

k k+n
[1 2........n][1.......n]......... ..[1..n]
k = 1 2 0...0 2 0......................0 A 1,1 0 J (1,1)


=


An ,1

En el Apndice 2 se pueden ver las contribuciones a Q de los restantes tres vrtices

resultando:

q El vrtice i=1, j=m ser la primera fila de Sm

q El vrtice i=n, j=1 ser la ltima fila de S 1
q El vrtice i=n, j=m ser la ultima fila de la matriz Q

Como resumen se podra ver a que lugares de la matriz Q contribuye cada zona de la
grilla

VERTICE (1,1) A1,1 0 J1,1

1
BORDE ( j = 1)
k =i Ai ,1 0 J i ,1
An ,1 0 J n ,1
n VERTICE (n,1)

n +1 BORDE (i = 1) A1, 2 0 J1, 2
GRILLA INTERIOR A
k =i+n i,2 0 J i ,2

2n BORDE (i = n) An , 2 0 J n , 2
M M M M
k = 1 + n ( j 1) BORDE (i = 1) Ai, j 0 Ji, j
k = i + n( j 1)
GRILLA INTERIOR x A = J
i, j 0 i, j
k = n + n( j 1) BORDE (i = n)
An, j
0 Jn, j


M M M M
k = 1 + n( m 2) BORDE (i = 1) A1, m 1 0 J1, m 1
k = i + n( m 2) GRILLA INTERIOR A
i , m 1 0 J i , m 1
k = n + n(m 2) BORDE (i = n) An , m 1 0 J n , m 1
nm n + 1 VERTICE (1, m)
A1, m 0 J1, m
k = i + n(m 1) BORDE ( j = m)


Ai , m


nm

0 J i ,m
VERTICE ( n, m) 0 J n , m
n , m
A

Gonzalo Casaravilla: [email protected] pag. 8 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

3.4 Resolucin del sistema de ecuaciones tridiagonales

Como ya se ha establecido, la resolucin del problema pasa por obtener S a partir de la

ecuacin QS=G. Una primera aproximacin al problema sera invertir la matriz Q y as
obtener S. Si se considera la discretizacin adoptada donde dentro de la grilla se
obtuvieron aproximaciones de segundo orden, y se tiene en cuenta que la inversin de
matrices convencional por el mtodo de Gauss-Jordan necesita para matrices
tridiagonales aproximaciones de tercer orden [12], se concluye que esta no es la forma
adecuada de resolver este problema. Una forma podra ser aumentar el orden de la
discretizacin, pero eso har aumentar exponencialmente el tamao del sistema de
ecuaciones a resolver. Para resolver esta dificultad existen complicados mtodos para
resolver sistema de ecuaciones con matrices tridiagonales que no sern analizados en
este trabajo.

En este punto surge la pregunta: Hay alguna forma de seleccionar a priori la grilla y
tener asegurada la convergencia numrica del mtodo de diferencias finitas?. Para
intentar contestar esta pregunta se intentar analizar si el criterio de Von Neumann
responde la pregunta planteada.

3.5 Criterio de estabilidad de Von Neumann

Supongamos que se puede expresar

N
A( xi , y j ) = Ai( j ) = E w( j ) e si s = 1
w=1

El criterio de estabilidad de Von Neumann establece que secuencia de valores E debe

ser montona decreciente por lo que se debe cumplir que

E w( j +1) w
G = ( j)
1 =
Ew N

Aplicando esta definicin a la discretizacin de la ecuacin de Poisson resulta

E w( j ) e s (i 1) 2 E w( j ) e si + E w( j ) e s ( i+1) Ew( j 1) e si 2 Ew( j ) e si + E w( j +1) e si

+ =0
x 2 y 2
que dividiendo por E w( j ) e si resulta
E w( j 1) E w( j ) E w( j +1)
( j 1) ( j +1) 2 +
e s
2+e + s
E 2E + E
( j)
2 cos 2 E ( j)
E ( j)
E w( j )
+ w w w
= + w w
=0
x 2 E w( j ) y 2 x 2 y 2

E w( j +1) E w( j )
Si nos atenemos a la definicin de G y asumimos que resulta
E w( j ) E w( j 1)

2 cos 2 G 1 2 + G y 2
+ =0 G 2 2(a + 1)G + 1 = 0 con a = (cos 1)
x 2 y 2 x 2

Gonzalo Casaravilla: [email protected] pag. 9 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

La solucin de la ecuacin en segundo grado en G es G = a + 1 a (a + 2)

por lo que se puede ver que no hay forma de acotar a (y por tanto los incrementos en x e
y) de tal modo que el mdulo de G sea siempre menor que 1 para todo w. Resulta por
tanto que no se puede seleccionar a priori la grilla que asegure convergencia al mtodo
de diferencias finitas para la resolucin la ecuacin de Poisson con discretizacin de
segundo orden de aproximacin.

Revisando los criterios de seleccin de grilla de los mtodos numricos como el FDTD,
resulta evidente, y posiblemente lo nico que se puede hacer, la adopcin de criterios
prcticos como es el de Courant, por el cual se selecciona una grilla que no avance ms
que una fraccin de longitud de onda de la frecuencia mas alta presente en el fenmeno
estudiado [14][15][16].

4 Qu es lo que se usa en la actualidad?

A los efectos de saber como se resuelven los problemas magnetoestticos en algunas

herramientas de clculo del mercado se analizaron algunas de ellas.

4.1 Partial Differential Equations Toolbox de Matlab

Esta herramienta cuenta con un interesante conjunto de funciones que resuelven

ecuaciones en derivadas parciales, ya sean hiperblicas, parablicas o elpticas. En
particular estas ltimas son las mas desarrolladas contndose incluso con una funcin
explcita para resolver la ecuacin de Poisson en 2D en coordenadas rectangulares por
diferencias finitas para el caso de fronteras de Dirichlet. Usa para este caso mtodos de
resolucin numricos combinados de transformadas de senos y resolucin de sistemas
tridiagonales.

Para otros casos, en particular con condiciones de borde de Neumann, resuelve

numricamente usando Elementos Finitos.

En el Anexo 3 se muestran algunas de las potencialidades de la herramienta.

4.2 Filds: http://www.crbond.com/fields.htm

Usa Diferencias Finitas (las define como avanzadas) con Interpolacin Cubica.

4.3 Alston: http://www.technology.alstom.com/en/rdpolicy/consultancy/slim

Usa Elementos Finitos.

4.4 Ansoft: http://www.ansoft.com/products/em/max2d/index.cfm

Usa Elementos Finitos.

Gonzalo Casaravilla: [email protected] pag. 10 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

5 Conclusiones

Se han planteado las ecuaciones fsicas de la magnetoesttica derivadas de dos de las

cuatro ecuaciones de Maxwell con sus particulares condiciones de borde. En problemas
planos se arriba a PDE elpticas (Poisson) con condiciones de borde de Neumann.

Se analiz como se resolvera el problema con Diferencias Finitas arribndose a la

conclusin de que si bien el sistema de ecuaciones tridiagonal es relativamente sencillo
de plantear, su resolucin numrica puede presentar problemas. En particular el criterio
de estabilidad de Von Neumann no ofrece respuestas concluyentes.

Del anlisis del estado del arte actual se concluye que Difrencias Finitas son adecuadas
y son utilizadas casi exclusivamente para resolver problemas electromagnticos
transitorios en los que interviene el tiempo mediante el mtodo Diferencias Finitas en el
Dominio del Tiempo. Respecto a como se resuelven problemas de magnetoesttica (o
electroesttica), el mtodo de Diferencias Finitas es el mas usado.

Finalemente se analizararon las herramientas Matlab para la resolucin de ecuaciones

diferenciales en derivadas parciales, obtenindose resultados mas que satisfactorios que
terminaron contestando la motivacin original de este trabajo.

Gonzalo Casaravilla: [email protected] pag. 11 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

6 Bibliografa

[1] Fundamentos de la Teora Electromagntica. J. R. Reitz, F. J. Milford

[2] Electromagntics. J.D. Kraus.

[3] Numerical Solutions of Initial Boundary Value Problems Involving Maxwells

Equations in Isotropic Media. K. S. Yee. IEEE Transactions on antennas and
propagation, May 1966.

[4] A Memory-Efficient Formulation of the Finite-Difference Time-Domain Method for

the Solution of Maxwell Equations. G. D. Kondylis, et all. IEEE Tranasactions on
Microwave Theory and Thecniques, July 2001.

[5] Finite-Difference Time-Domain Literature Database. http://www.fdtd.org/

[6] FDTD Method. A. Dunn

http://www.nmr.mgh.harvard.edu/~adunn/papers/dissertation/node31.html

[7] Clculo de las Fuerzas de Cortocircuito en Transformadores. A. Azziz, A. Portillo.

Mak S.A.

[8] Leakage Flux and Force Calculation on Power Transformer Winding under Short-
circuit: 2D and 3D Methods based on the Theory of Images and the Finite Element
Method Compared to Measurements. A. Kladas etc all. IEEE Transactions on
Magnetics, Sept. 1994.

[9] Applications of Parallel Computers. D. H. Bailey.

http://www.nersc.gov/~dhbailey/cs267/Lectures/Lect_11_2000.pdf

[10] Differential Equations. M. Guidry, M. Strayer.

http://csep1.phy.ornl.gov/guidry/phys594/lectures/diff_equations/pdeelliptic/pdeelliptic.
html

[11] Mtodos Numricos para la Resolucin de Ecuaciones Diferenciales. S.

Nesmachnow. http://www.fing.edu.uy/inco/cursos/numerico/mned.html

[12] Introduction to Computational Physics. R. Fitzpatrick

http://farside.ph.utexas.edu/teaching/329/lectures/lectures.html

[13] Partial Differential Equations. Anam.

http://www.ccs.uky.edu/~douglas/Classes/cs521-s02/pde/pde.pdf

[14] Modelling requirements for irregular grids

http://trigrid.sourceforge.net/tgdocs/node4.htm

[15] Charlie Chen, Tae-Woo Lee, Narayanan Murugesan, Susan Hagness. Generalized
FDTD-ADI: An Unconditionally stable Full-Wave Maxwells equations. Solver for
VLSI Interconnect Modelinghttp://vlsi.ece.wisc.edu/research/2000iccad01.pdf

[16] Brunel University, UK. The Finite Difference Time Domain Algorithm
http://www.nmr.mgh.harvard.edu/~adunn/papers/dissertation/node32.html

Gonzalo Casaravilla: [email protected] pag. 12 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Apndice 1: Contribucin del las fronteras j=m, i=1 e i=n a la matriz Q.

De la misma forma que se analiz j=1, para la frontera j=m se discretiza igual con la
salvedad que la discretizacin en y sufre un cambio en el orden en los subndices (se
calcula el valor de Ai,m-1 en funcin de A i,m) resultando

A 2 A y 2
Ai ,m1 = Ai ,m + .( y ) + 2 . + O (y 3 )
y i ,m y i , m 2
123
=0

A 2
(Ai,m1 Ai ,m )
2
y i ,m
2
y 2

Ai 1,m 2 Ai ,m + Ai +1,m
+2
(A i ,m 1 Ai ,m )
= 0 J ( x, y ) = 0 J (i, m)
x 2 y 2

por lo que quedan determinadas las filas de Q correspondientes a Sm con k de (mn-n+2)

a (mn-1).

k-n k
[1..n].........[1....i....n ]. [1..........i - 1 i i + 1.....n ]
A1, m


k = i + n(m - 1) 0...............0 2 0.........0 0....0 Ai ,m = 0 J (i, m)


An, m

Con razonamientos similares se puede ver como quedan las filas de Q correspondientes
a las condiciones de borde verticales, i=1 e i=n.

Por ejemplo para i=1, con j variando de 2 a (m-1) tenemos

(A A1, j ) A1, j 1 2 A1, j + A1, j +1

+ = 0 J ( x, y ) = 0 J (1, j )
2, j
2
x 2
y 2

Por lo tanto podemos determinar para los (n-2) puntos en consideracin de dicho borde
las filas de la matriz Q. Como i=1, las filas que estaremos determinando son las filas
k=1+n(j-1) lo cual corresponde a las primeras filas de S 2, S3..., Sm-1.

Gonzalo Casaravilla: [email protected] pag. 13 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

kn k k+n
[1..n]....[1......n] [1. ..........n ] [1......n]...[1..n]
k = 1 + n(j - 1) 0...0 0.....0 2 0....0 0.............0 A1, j 0 J (1, j )


=


An , j

Luego en la otra frontera i=n, queda

(A
n 1, j An, j ) An , j 1 2 An , j + An , j +1
2 + = 0 J ( x, y ) = 0 J (n, j )
x 2
y 2

cual corresponde a las ltimas filas de S 2, S3..., Sm-1.

kn k k+n
[1..n]....[1......n] [1. ..........n ] [1......n]...[1..n]
A1, j 0 J (1, j )


=


k = n + n(j - 1) 0............0 0....0 2 0....0 0.....0 An, j

Apndice 2: Contribucin de los vrtices restantes a la matriz Q.

Para el vrtice i=1, j=m tendremos

(A A1,m ) (A 1, m 1 A1,m )
+2 = 0 J ( x, y ) = 0 J (1, m)
2 ,m
2
x 2
y 2

que ser la primera fila de Sm que quedar

k-n k
[1..n]............... [1 ..........n][1..............n]
k = 1 + n(m - 1) 0..................... 0.....0 2 0....0 A 1,m 0 J (1, m)


=


An ,m

Gonzalo Casaravilla: [email protected] pag. 14 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Luego para el vrtice i=n, j=1 tendremos

(A
n 1,1 An,1 ) (A An,1 )
2 +2 n, 2
= 0 J ( x, y ) = 0 J (n,1)
x 2
y 2

que ser la ltima fila de S1 que quedar

k k+n
[1 ...............n][1.............n]....[1..n]
A1,1


=


k = n 0...0 2 0.....0 2 0....0 An,1 0 J (n,1)

Finalmente tendremos el vrtice i=n, j=m cuya discretizacin queda

(A
n 1, m An, m ) (A n , m 1 An, m )
2 +2 = 0 J ( x, y ) = 0 J ( n , m )
x 2 y 2

que ser la ltima fila de Q que quedar sus elementos as dispuestos

k-n k
[1.....n]....[1.....n][1................n]
A 1,m


=


k = nm 0............0 2 0.....0 2 An ,m 0 J (n, m)

Gonzalo Casaravilla: [email protected] pag. 15 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Apndice 3: Partial Differential Equations Toolbox de Matlab

1.- PDETOOL

Bsicamente con esta herramienta se puede intentar resolver muchos problemas en 2D

estticos convencionales en donde aparecen ecuaciones en derivadas parciales. Llega al
extremo de poderse configurarse el tipo de problema a resolver (Generic Scalar, Generic
System, Structural Mechanics, Electrostatics, Magnetostatics, AC Power
Electromagnetics, Conductive Media DC, Heat Transfer, Diffusion) cambiando la forma
de pedir y mostrar los datos de entrada de acuerdo a la aplicacin seleccionada.

Permite determinar diversas fronteras (crculos, rectngulos, elipses, etc..) y para cada
una de las fronteras permite seleccionar si es de tipo Dirichlet o Neumann.

Por ejemplo para resolver un problema como el plantado en este trabajo, la mayor
dificultad pasa en los casos en que J es variable con las coordenadas x e y, dado que la
especificacin de la fuente J(x,y) debe ser una expresin analtica. Si por el contrario y
ms comn, J es constante en ciertas regiones, se pueden definir las mismas y
especificar cuanto vale J en cada una de ellas.

Por ejemplo supongamos que tenemos la ventana de transformador en donde hay dos
arrollamientos iguales tal como muestra la figura (las figuras son las pantallas de la
propia interface de PDETOOL)

Supongamos que el arrollamiento 3 es atravesado por corriente entrante con una

densidad de corriente tal que oJ=1 y que el arrollamiento 2 tiene oJ=-1 (en un
transformador la integral de J en un corte de ventana debe ser nula).
Por tanto una forma de especificar J analticamente podra ser:

J=(1/4)*(sign(x+0.5)-sign(x+0.3)).*(sign(y+0.3)-sign(y-0.5))+(-1/4)*(sign(x+0.1)-
sign(x-0.1)).*(sign(y+0.3)-sign(y-0.5))

Notar que esta forma permitira especificar tantos bobinados como sean necesarios.

Gonzalo Casaravilla: [email protected] pag. 16 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Sin embargo, para este caso es ms fcil especificar para la zona 1, J=0, para la zona 2
J=-1 y para la zona J=+1.

Luego la herramienta permite seleccionar la malla de los Elementos Finitos deseada,

pudindose regular la densidad de la misma.

Luego se obtiene el resultado buscado pudindose especificar que variables se desean

ver, en 2D ya sea A, H, B, contornos de cada una o diagrama de vectores absolutos o
normalizados de los campos B o H.

Gonzalo Casaravilla: [email protected] pag. 17 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Tambin dispone de herramientas de visualizacin en 3D ya predeterminadas de cualquiera de las

variables del sistema

2.- POISOLV

Esta es la funcin Matlab que tiene implementada la resolucin de la ecuacin de

Poisson en una grilla rectangular mediante diferencias finitas pero con condiciones de
borde de Dirichlet. La forma de especificarle el problema a la funcin est normalizado
de acuerdo a como trabaja toda la biblioteca PDE de Matlab siendo de gran ayuda el uso
de la propia PDETOOL a los efectos de generar las especificaciones geomtricas,
condiciones de borde y grillas del problema (PDETOOL exporta todas estas
definiciones al entorno Matlab).

Por ejemplo, en el mismo problema recin visto, se puede ver la forma de especificar las
condiciones de frontera las cuales se pueden ver en la definicin de la matriz b que
sigue:

% Boundary Condition matrix

% Esto define las condiciones de borde en cada frontera
% u es el potencial vector que cumple -div(grad u)=mu J con
B=rot u
% n . grad(u) + qu = g (n es la normal a la superficie)
% h u = r
% Cada columna es para cada borde
% 1.- N = 1 si es una frontera
% 2.- M (1 Dirichlet, 0 Neumann)
% 3.- cantidad de caracteres de q
% 4.- cantidad de caracteres de g

Gonzalo Casaravilla: [email protected] pag. 18 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

% 5.- cantidad de caracteres de h

% 6.- cantidad de caracteres de r
% 7.- char equivalente del valor asignado a q (48 es 0)
% 8.- char equivalente del valor asignado a g
% 9.- char equivalente del valor asignado a h (49 es 1)
% 10.- char equivalente del valor asignado a r
% tener en cuanta que si se quiere h=22, son dos char etc..
% etc.

b =[1 1 1 1 % N (si es 1 es frontera)

1 1 1 1 % M (Si es 1 es Dirichlet)
1 1 1 1 % un caracter
1 1 1 1
1 1 1 1
1 1 1 1
48 48 48 48 % equivale al 0
48 48 48 48
49 49 49 49 % 1 de Newman pero igual r es 0.
48 48 48 48 ];

Luego se puede ver como se especifica geomtricamente la propia frontera

% Por cada columna se define una recta

% 1.- es 1 si es una frontera, o si es una subzona.
% 2.- x1
% 3.- x2
% 4.- y1
% 5.- y2
% 6.- O si es frontera
% 7.- 1 si es frontera
g =[2.0000 2.0000 2.0000 2.0000
-0.7000 0.3000 0.3000 -0.7000
0.3000 0.3000 -0.7000 -0.7000
0.6000 0.6000 -0.4000 -0.4000
0.6000 -0.4000 -0.4000 0.6000
0 0 0 0
1.0000 1.0000 1.0000 1.0000
];

Luego solo restara definir la grilla de clculo para lo cual se cuenta con una funcin
matlab (poimesh) que trabaja de la siguiente menera:

[p,e,t]=poimesh(g,nx,ny);
% nx y ny es la cantidad de subdivisiones en x e y
% de la grilla.
% p es un par de filas con los valores x,y de todos los
% nodos de
% la grilla ordenado desde la izquierda abajo recorriendo
% las filas
% desde izquierda a derecha etc..
% En cada columna de t estn los vrtices de cada tringulo
% de la mesh

Gonzalo Casaravilla: [email protected] pag. 19 de 20

Monografa del curso: Mtodos Numricos para la resolucin de Ecuaciones Diferenciales Noviembre 2002

Por ltimo slo resta especificar el trmino independiente de la ecuacin de Poisson que
como ya se ha dicho debe ser en forma analtica o por subzonas.

Finalmente solo quedara por invocar la propia funcin POISOLV

u=poisolv(b,p,e,t,f);

Si elegimos la grilla adecuadamente, dar el mismo resultado que la funcin

u=assempde(b,p,e,t,1,0,f) que no es otra que la funcin genrica de clculo por
Elementos Finitos que usa la propia PDETOOL para hacer los clculos.

Se destaca que existen funciones Matlab de visualizacin en 2D, contornos, en 3D etc,

que usan los mismos b, p, e, t etc. por lo que el entorno de trabajo, luego de
comprendido el protocolo, resulta sumamente fcil de utilizar.

Se contrastaron resultado de usar POISOLV (diferencias finitas) y PDETOOL

(elementos finitos) sobre un mismo caso, encontrndose resultados coincidentes luego
de achicadas las respectivas grillas adecuadamente. No se muestran los resultados por
carecer de significado fsico ya que se debi imponer condiciones de borde de Dirichlet
dada la limitante al respecto que tiene POISOLV.

Solo a los efectos de ilustrar respecto a la potencialidad de esta herramienta de clculo

se muestran resultados hechos con PDETOOL para un transformador de varios
arrollados. Por ejemplo la herramienta sirve para observar como por efecto de
asimetras constructivas existen campos tangenciales en las bobinas que traern como
consecuencias esfuerzos dinmicos que deben ser acotados durante el diseo de los
mismos.

Gonzalo Casaravilla: [email protected] pag. 20 de 20