La función pdepe de MATLAB

Autor: O. A. Hernández-Flores ([email protected])

Afiliación: Escuela de Sistemas Biológicos e Innovación Tecnológica, Universidad Autónoma Benito Juárez de Oaxaca. Av. Universidad s/n C.P. 68120 Oaxaca de Juárez, Oaxaca, México.

Resumen:

Se muestra como utilizar el solver pdepe de MATLAB. El sistema de estudio es una placa infinita superconductora.

• Usage: Secondary publications and information retrieval purposes.
• Structure: You may use the description environment to structure your abstract; use the optional argument of the \item command to give the category of each item.

1. La función 'pdepe'

La función 'pdepe' resuelve sistemas de EDP parabólicas y elípticas en una variable espacial $x$ y el tiempo $t$, de la forma

$\begin{aligned} c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t} &= x^{-m}\frac{\partial}{\partial x}\left( x^m f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right) \\ &+ s\left(x,t,u,\frac{\partial u}{\partial x}\right). \end{aligned}$

Las PDE son válidas para $t_0 \le t \le t_f$ y $a \le x \le b$. El intervalo $[a, b]$ debe ser finito. La constante $m$ puede ser $0, 1 \mathrm{o} 2$, correspondientes a una simetría de placa, cilíndrica o esférica, respectivamente. Si $m> 0$ entonces $a \ge 0$ debe satisfacerse. En la ecuación anterior, $f$ es un término de flujo y $s$ es un término fuente. El término de flujo debe depender de $\partial u/\partial x$. El acoplamiento de las derivadas parciales con respecto al tiempo se restringe a la multiplicación por una matriz diagonal $c(x,t,u,\partial u/\partial x)$. Los elementos diagonales de esta matriz son idénticamente cero o positivos. Un elemento que es idénticamente cero corresponde a una ecuación elíptica y, de lo contrario, a una ecuación parabólica. Debe haber al menos una ecuación parabólica. Un elemento de c que corresponde a una ecuación parabólica puede anularse en valores aislados de $x$ si son puntos de la malla. Se permiten discontinuidades en $c$ y/o $s$ en las interfaces del material siempre que se coloque un punto de malla en cada interfaz.

En el tiempo inicial $t = t_0$, para todo $x$, los componentes de la solución satisfacen las condiciones iniciales de la forma:

$u(x,t_0)=u_0(x).$

En las fronteras $x = a$ y $x = b$, para todo $t$, los componentes de la solución satisfacen una condición de frontera de la forma:

$p(x,t,u) + q(x,t)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0,$

aquí $q(x,t)$ es una matriz diagonal con elementos que pueden ser idénticamente cero. Tenga en cuenta que las condiciones de contorno se expresan en términos de $f$ en lugar de la derivada parcial de $u$ con respecto a $x$, $\frac{\partial u}{\partial x}$, además, de los dos coeficientes, solo $p(x,t,u)$ puede depender de $u$.

2. Sintaxis

La PDE solver tiene la siguiente sintaxis:

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)

donde los argumentos de entrada son:

• m: Especifica la simetría del problema, placa, cilindro, esfera m=0,1,2.
• pdefun: Función que define el sistema de EDP. Se escribe en terminos de c,f y s, y tiene la forma

[c,f,s] = pdefun(x,t,u,dudx)

donde x y t son escalares y u y dudx son vectores. Por otro lado, c,f y s son vectores columna. En c se almacenan los elementos de la matriz $c(x,t,u,\partial u/\partial x)$.

• icfun: Es una funcion que evalua las condiciones iniciales. Es de la forma

u = icfun(x)

cuando es llamada con el argumento $x$, la función icfun evalúa y regresa los valores iniciales evaluados en $x$ en un vector columna u.

• bcfun: Funcion que evalua los terminos p y q de las condiciones a la frontera. Es de la forma

[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

donde ul es la solución aproximada en la frontera izquierda xl=a, mientras que ur es la solución aproximada en la frontera derecha xr=b.

• xmesh: Vector [x0, x1, ..., xn] que especifica los puntos en los que se solicita una solución numérica para cada valor en tspan. x0 y xn corresponden a y b, respectivamente.
• tspan: Vector [t0, t1, ..., tf] que especifica los puntos en los que se solicita una solución para cada valor en xmesh. t0 y tf corresponden a $t_0$ y $t_f$, respectivamente.

El argumento de salida sol es una matriz tridimensional, tal que:

• sol (:,:, k): Componente k de la solución u.
• sol (i,:, k): Componente k de la solución en el momento tspan (i) y los puntos de malla xmesh (:).
• sol (i, j, k): Componente k de la solución en el tiempo tspan (i) y el punto de malla xmesh(j)

3. La placa infinita

Nosotros consideramos una placa infinita de espesor $d$ definida por $\vert{}y,z\vert{}<\infty$ y $\vert{}x\vert{}\le d/2$. La placa esta sometida a la interacción de un campo magnético externo homogéneo incidiendo en forma paralela al eje $z$. La intensidad del campo es $B_a=\mu_0H_a=at$ donde $a$ es el sweet rate. En esta denominada geometría paralela, la dinámica del campo magnético se describe con la ley de inducción de Faraday,

$\frac{\partial Bz}{\partial t} = \frac{1}{\mu_0}\frac{\partial}{\partial x} \left(\rho(x,t,B_z)\frac{\partial B_z}{\partial x} \right),$

esta ecuación en derivadas parciales esta sujeta a las condiciones inicial (CI) y de frontera (CF):

$\begin{aligned} \text{CI:} \quad & B_z(x,t=0) = 0, \\ \text{CF:} \quad & B_z(x=0,t) = B_z(x=d,t) = B_a. \end{aligned}$

Para hallar la ecuación anterior se ha utilizado la relación constitutiva $\mathbf{E} = \rho \mathbf{j}$ y la ley de Ampere $\mu_0\mathbf{j} = -\partial B_z/\partial x\hat{\mathbf{y}}$. Aquí hemos considerado que el material es isotrópico, por lo que la resistividad $\rho = \rho(B_z,T)$ es una función escalar y en general depende del campo magnético y de la temperatura, por otro lado, se asume que el primer campo critico $H_{c1}\ll H_a$ lo que permite despreciar los efectos de las barreras geométricas, las barreras de Levingston y las corrientes de Meissner, con ello, la relación constitutiva $B_z=\mu_0H_z$ es valida. Se asumirá que el campo externo varia linealmente con el tiempo $H_a=\alpha t+\beta$.

La resistividad se modela como:

$\rho = \rho_0\frac{j_{c0}}{j_c}S.$

aquí $j_{c0}$ es la intensidad máxima que puede alcanzar la densidad corriente crítica, mientras que $S$ es una función interpolante, que describe de forma suave y continua la transición del régimen de flux creep al régimen de flux flow. En particular, se utiliza:

$S = \frac{(j/j_c)^m}{1+(j/j_c)^m}.$

donde $m$ es un parámetro de ajuste conocido como exponente de creep. La magnitud de la densidad de corriente se calcula con:

$j_y=-\frac{\partial B_z}{\partial x}.$

Finalmente, la densidad de corriente corriente crítica es vital en la descripción de las propiedades electromagnéticas de un material superconductor, ya que esta cantidad determina la dinámica transciente o estacionaria. Para nuestros propósitos, la densidad de corriente crítica se escribe como:

$j_c = \frac{j_{c0}}{(1+B/B^*)^n}.$

donde $B^*$ y $n$ son parámetros de ajuste.

4. Ecuaciones adimensionadas

Las tres cantidades básicas de nuestro estudio se pueden escribir en términos de cantidades características, tal y como se muestra a continuación:

$x = x_0\bar{x}, \quad t = t_0\bar{t}, \quad B_z = B_0\bar{B}_z.$

La elección de las cantidades características es un tanto arbitraria, sin embargo, es recomendable que nuestra elección refleje la física del sistema de estudio. En nuestro caso se eligieron: 1) el espesor de la muestra $x_0=d$, 2) $B_0=1\mathrm{T}$ y 3) $t_0=1\mathrm{s}$, un tiempo característico del orden del tiempo de observación o el tiempo que dura el experimento.

$\frac{\partial \bar{B}z}{\partial \bar{t}} = \frac{t_0}{t_1}\frac{\partial}{\partial \bar{x}} \left(\bar{\rho}(\bar{x},\bar{t},\bar{B}_z)\frac{\partial \bar{B}_z}{\partial \bar{x}} \right)$

donde $t_1= \mu_0d^2/\rho_0$. El tiempo intrínseco $t_1$ delimita el momento en que la respuesta magnética comienza a relajarse. Las respectivas condiciones inicial y de frontera lucen como sigue:

$\begin{aligned} \text{CI:} \quad & \bar{B}_z(\bar{x},\bar{t}=0)=0, \\ \text{CF:} \quad & \bar{B}_z(0,\bar{t})=\bar{B}_z(1,\bar{t})=\bar{B_a}, \end{aligned}$

donde $\bar{B}_a=\bar{\alpha} \bar{t}+\bar{\beta}$. El resto de las ecuaciones lucen como sigue:

$\begin{aligned} \bar{\rho} &= \left(1+\dfrac{\bar{B}}{\bar{B}^*}\right)^n\frac{\left(\dfrac{\bar{j}}{\bar{j}_c}\right)^m}{1+\left(\dfrac{\bar{j}}{\bar{j}_c}\right)^m} \\ \bar{j}_c &= \frac{\bar{j}_{c0}}{(1+\bar{B}/\bar{B}^*)^n}, \\ \bar{j} &= \vert{}\bar{j}_y\vert{} = \left\vert{}-\frac{\partial \bar{B}_z}{\partial \bar{x}}\right\vert{}. \end{aligned}$

En la tabla se resumen los valores de los parámetros del modelado.

Recordar que $\mu_0=4\pi\cdot 10^{-7}\mathrm{Tm/A}=4\pi\cdot 10^{-7}\mathrm{H/m}$ donde $\mathrm{H} = \mathrm{Tm^2/A} = \mathrm{\Omega s}$.

5. Implementación

Figura 1: Señal triangular

El campo aplicado

Es nuestro objetivo trazar una curva de histeresis con un comportamiento lineal del campo aplicado, este campo es un requerimiento para las condiciones a la frontera. Para nuestros propósitos se utiliza la siguiente señal triangular:

$\bar{B}_a = 4\dfrac{A}{P}\left\vert{}\mathrm{mod}\left(\bar{t}-\dfrac{P}{4},P\right)-\dfrac{P}{2}\right\vert{}-A;$

donde $A$ es la amplitud y $P$ el periodo. En este caso la amplitud corresponde al máximo valor del campo aplicado $\bar{B}_a=\bar{B}_{max}$. La figura 1 muestra una señal triangular, para su construcción se eligio $\bar{\alpha} = 10^{-3}$ y un campo máximo $\bar{B}_{max}=1$, con esta información el tiempo necesario para alcanzar el campo máximo es $t_{max}=\bar{B}_{max}/\bar{\alpha}$, mientras que el periodo sera entonces $P = 4t_{max}$, por lo que para completar la histeresis se requiere un tiempo $t_H = 5t_{max}$. El siguiente código en Matlab produjo la grafica de la señal triangular.

Bmax  = 1;      % Valor maximo del campo 
alpha = 1e-3;   %Sweet rate
tmax  = Bmax/alpha; %Tiempo maximo
NoPuntos = 1000; 
t = linspace(0,5*tmax,NoPuntos);
A = Bmax;   % Amplitud
P = 4*tmax; % Periodo
Ba = (4.*A./P).*abs(mod((t-P./4),P)-P./2)-A;
plot(t,Ba,'LineWidth',2)
xlabel('Tiempo $\bar{t}$','Interpreter','Latex')
ylabel('Campo aplicado $\bar{B}_a$','Interpreter','Latex')

Rescribir la EDP en la forma estandar

Ahora bien, para utilizar el solver pdepe se requiere escribir nuestra ecuación diferencial parcial adimensionada en la forma de la ecuación estándar, esto es:

$\frac{t_1}{t_0}\frac{\partial \bar{B}z}{\partial \bar{t}} = \frac{\partial}{\partial \bar{x}} \left(\bar{\rho}(\bar{x},\bar{t},\bar{B}_z)\frac{\partial \bar{B}_z}{\partial \bar{x}} \right)$

y con ello identificamos términos:

$\begin{aligned} c(x,t,u,\dfrac{\partial u}{\partial x}) &= \dfrac{t_1}{t_0} \\ f(x,t,u,\dfrac{\partial u}{\partial x}) &= \bar{\rho}(\bar{x},\bar{t},\bar{B}_z)\frac{\partial \bar{B}_z}{\partial \bar{x}}, \\ s(x,t,u,\dfrac{\partial u}{\partial x}) &= 0. \end{aligned}$

Dado que $x^{-m}=1$ se tiene que $m=0$. Identificados $c,f,s,m$ estamos listos para crear un código en Matlab.

Código

En vista de que tenemos varios parámetros de ajuste es conveniente llamar al solver pdepe dentro de una función madre, esto con la finalidad de definir los parámetros en ella y se comporten como variables globales para las subfunciones que se utilizan. El código completo se muestra a continuación:

function [x,sol]=pdepePlaca
m  = 0;
t0 = 1;
t1 = 1e-3;
jc0= 1;
Be = 1;
mm = 1;
nn = 40;
Sweep_Rate = 1e-3;
Bmax = 1;
tmax = Bmax/Sweep_Rate;
tsim = 5*tmax;
meshpoints = 65;
x = linspace(0,1,meshpoints);
t = linspace(0,tsim,500);
%
sol = pdepe(m,@pde_pde,@pde_ic,@pde_bc,x,t);
%
%% Ecuacion diferencial parcial
function [c,f,s] = pde_pde(~,~,Bz,DBzDx)
c = t1/t0;
rho = resistividad(Bz,DBzDx);
f = rho*DBzDx;
s = 0;
end
%% Condiciones iniciales
function u0 = pde_ic(x)
u0 = zeros(size(x));
end
%% Condiciones a la frontera
function [pl,ql,pr,qr] = pde_bc(~,Bzl,~,Bzr,t)
P = 4*tmax;
Ba = (4*Bmax/P)*abs(mod((t-P/4),P)-P/2)-Bmax;
pl = Bzl-Ba;
ql = 0;
pr = Bzr-Ba;
qr = 0;
end
%% Resistividad 
function rho = resistividad(Bz,DBzDx)
jc  = jc0/(1+abs(Bz)/Be)^mm;
j   = abs(DBzDx);
S   = ((j/jc)^nn)/(1+(j/jc)^nn);
rho = ((1+abs(Bz)/Be)^mm)*S;    
end
%
end

Código trabajando

En esa sección se muestran gráficas obtenidas con el solver pdepe. En la gráfica se observa un conjunto de perfiles de Bz que se obtuvieron con las siguientes instrucciones:

% Perfiles de Bz
u = sol(:,:,1);
subplot(2,2,[1 3]),plot(x,u(10:20:end,:),'-');
title('Perfiles de B_z');
xlabel('$x/d$','FontSize',16,...
		'Interpreter','LaTeX');
ylabel('$B_z/B_0$','FontSize',16,...
		'Interpreter','LaTeX');

La instrucción pdeval evalúa numéricamente la solución de la EDP utilizando la salida de pdepe. El siguiente código despliega la solución y la derivada de la solución:

m=0;
[Bz,DBzDx] = pdeval(m,x,u(115,:),x);
subplot(2,2,2),plot(x,Bz,x,-DBzDx)
xlabel('$x/d$','FontSize',16, ...
		'Interpreter','LaTeX');
ylabel('$B_z/B_0\;(j_y/j_0)$','FontSize',16,...
		'Interpreter','LaTeX');
axis([0 1 -0.5 1.5])

La magnetización se puede hallar evaluando numéricamente lo siguiente:

$M_y = -H_a + \dfrac{1}{d}\int_0^d\dfrac{B_z}{\mu_0}dx$

o bien, en su versión adimensional:

$\bar{M}_y = -\bar{B}_a + \int_0^1\bar{B}_zd\bar{x}.$

El siguiente código despliega en un gráfico la magnetización:

dx = x(2)-x(1);
L=length(u(:,1));
Mz = zeros(L,1); 
Ba = u(:,1);
for i=1:L
[Bz,~] = pdeval(m,x,u(i,:),x);
Mz(i) = -Ba(i) + dx*trapz(Bz);
end
subplot(2,2,4),plot(Ba,Mz,'-');
xlabel('$B_a/B_0$','FontSize',16,...
		'Interpreter','LaTeX');
ylabel('$M_z/H_0$','FontSize',16,...
		'Interpreter','LaTeX');

Figura 2: Perfiles de $B_z$, densidad de corriente $j_y$ y una curva de histeresis