Aimant et bobine : flux et force

2.a. Définition

Le principe d'un modèle bidimensionnel en géométrie plane et les équations associées sont exposés dans Flux magnétique généré par un aimant permanent.

Tous les champs sont supposés invariants par translation dans la direction (Oz) donc ils ne dépendent que des variables x et y.

Figure pleine page

La bobine est modélisée par deux parties en cuivre : l'une est le siège d'un courant électrique de densité $\overrightarrow{J}=J\overrightarrow{u_z}$, l'autre est le siège du courant opposé. La bobine comporte éventuellement un noyau en fer doux (ou en ferrite doux), supposé linéaire et de perméabilité magnétique relative ${\mu }_r$. La bobine et son noyau sont centrés sur l'origine du repère. Ce système est censé modéliser une bobine mais la variation dans l'espace du champ ne peut être identique à celle du champ généré par une bobine de forme cylindrique axisymétrique (le champ décroît en 1/r² à grande distance au lieu de 1/r³).

L'aimant est modélisé par un rectangle uniformément aimanté : son aimantation est $\overrightarrow{M}=M\overrightarrow{u_z}$. La position du centre C de l'aimant est variable. L'aimant est supposé idéal, ce qui fait que son aimantation est constante (indépendante de sa distance à la bobine). L'aimant est équivalent à deux courants de surfaces localisés sur ses faces perpendiculaires à (Oy). Le vecteur densité de courant surfacique est :

$$\overrightarrow{S_m}=\overrightarrow{M}\wedge \overrightarrow{n}=\pm M\overrightarrow{u_z}\text{(1)}$$

Le potentiel vecteur est de la forme :

$$\overrightarrow{A}=A_z(x,y)\overrightarrow{u_z}\text{(2)}$$

Il vérifie l'équation de Poisson :

$$\frac{{\partial }^2{A}_z}{\partial x^2}+\frac{{\partial }^2{A}_z}{\partial y^2}=-{\mu }_0{J}_z\text{(3)}$$

Dans le cuivre, on a $J_z=\pm J$. Partout ailleurs on a $J_z=0$, à l'exception des surfaces de l'aimant parcourues par le courant de surface $\overrightarrow{S_m}$, que l'on traite comme un courant volumique réparti sur une épaisseur égale à la distance de la maille des différences finies (voir plus loin).

Le champ magnétique est obtenu par

$$\begin{array}{ccc}\hfill & B_x=\frac{\partial A_z}{\partial y}\hfill & \hfill \text{(4)}\\ \hfill & B_y=-\frac{\partial A_z}{\partial x}\hfill & \hfill \text{(5)}\end{array}$$

Les lignes de champ magnétique sont définies par la relation différentielle :

$$B_{x}dy-B_{y}dx=0\text{(6)}$$

qui devient :

$$d{A}_z=0\text{(7)}$$

Les lignes de champs sont donc les lignes d'égales valeurs de A_z.

À l'interface entre le fer et l'air ou le cuivre, on a continuité de la composante normale de $\overrightarrow{B}$ et continuité de la composante tangentielle de $\overrightarrow{B}/{\mu }_r$ :

$$\begin{array}{ccc}\hfill & {\left(B_x(x=\pm \frac{L}{2})\right)}_{air}={\left(B_x(x=\pm \frac{L}{2})\right)}_{fer}\hfill & \hfill \text{(8)}\\ \hfill & (B_y(x=\pm \frac{L}{2}){)}_{air}=(\frac{1}{{\mu }_r}{B}_y(x=\pm \frac{L}{2}){)}_{fer}\hfill & \hfill \text{(9)}\\ \hfill & {\left(B_y(y=\pm \frac{e}{2})\right)}_{air}={\left(B_y(y=\pm \frac{e}{2})\right)}_{cuivre}\hfill & \hfill \text{(10)}\\ \hfill & (B_x(y=\pm \frac{e}{2}){)}_{air}=(\frac{1}{{\mu }_r}{B}_x(y=\pm \frac{e}{2}){)}_{cuivre}\hfill & \hfill \text{(11)}\\ \hfill & \hfill \text{(12)}\end{array}$$

Le potentiel vecteur subit donc une discontinuité de son gradient dans la direction perpendiculaire à la surface, l'autre composante du gradient étant continue. Partout ailleurs, le potentiel vecteur est continu.

2.b. Résolution par différences finies

La résolution numérique de l'équation de Poisson par la méthode des différences finies est expliquée dans Équation de Poisson à deux dimensions. La discrétisation de l'équation est expliquée dans Discrétisation de l'équation de Poisson vectorielle.

Le module python utilisé est présenté dans Équation de Poisson : programme Python.

Le domaine de calcul est un carré de côté L_d. La condition limite sur le bord du domaine est de type dirichlet (annulation du potentiel vecteur). Le domaine devra être assez grand par rapport à la taille du système pour que les effets de bords soient négligeables.

Le nombre de mailles dans chacune des deux directions (Ox) et (Oy) est 2ⁿ. Pour définir les objets, on utilise un maillage réduit de taille 2^p par 2^p avec p<n. Ce maillage réduit est le plus grossier utilisé dans l'algorithme multigrille. Il doit tout de même être relativement fin car le déplacement de l'aimant ne peut se faire que par saut égal à la longueur de la maille réduite.

Le système est défini par deux rectangles pour la bobine, un rectangle pour le noyau et un rectangle pour l'aimant. Les sommets de ces rectangles doivent être sur des nœuds de la maille réduite. Il est donc commode de définir les dimensions par des nombres entiers qui représentent les dimensions divisées par la largeur de maille réduite. Les longueurs des côtés des rectangles devront être des nombres entier pairs.

Les courants surfaciques sur les faces $y=y_c\pm \frac{b}{2}$ de l'aimant sont des sources localisées sur deux rangées de nœuds. Soit $\Delta x$ la taille de la maille (maillage fin). On définit sur une rangée de nœuds une densité volumique J_m telle que $J_{m}a\Delta x=Ma$ soit $J_m=M/\Delta x$.

Les calculs de flux et de forces se font pour une longueur $\ell$ dans la direction (Oz).

Le flux magnétique dans la bobine est :

$$\Phi =\frac{\ell }{S}({\iint }_{+}{A}_z(x,y)dxdy-{\iint }_{-}{A}_z(x,y)dxdy)\text{(13)}$$

où S est l'aire d'un rectangle de cuivre. Pour calculer le flux généré par l'aimant dans la bobine, il faut poser J=0.

La résultante des forces magnétique sur l'aimant est (Flux magnétique généré par un aimant permanent) :

$$\begin{array}{ccc}\hfill & F_x=\ell M{\int }_{x_c-a/2}^{x_c+a/2}(-B_y(x,y_c+b/2)+B_y(x,y_c-b/2))dx\hfill & \hfill \text{(14)}\\ \hfill & F_y=\ell M{\int }_{x_c-a/2}^{x_c+a/2}(B_x(x,y_c+b/2)-B_x(x,y_c-b/2))dx\hfill & \hfill \text{(15)}\end{array}$$

L'énergie du champ magnétique est calculée par l'intégrale suivante :

$$E_c=\frac{1}{2}\iint \overrightarrow{j}\cdot \overrightarrow{A}dv=\frac{1}{2}\ell \iint j_z{A}_{z}dxdy\text{(16)}$$

Cette énergie permet de calculer la force qui s'exerce sur un élément du système : il suffit de déplacer cet élément et de calculer le gradient par rapport à ce déplacement :

$$\overrightarrow{F}=\overrightarrow{grad}{E}_c\text{(17)}$$

Cette méthode permet en particulier de calculer les forces sur les pièces en ferromagnétique doux. On peut même calculer directement la force sur un ensemble formé d'un aimant et d'une pièce en fer doux, en déplaçant ces deux éléments ensemble. L'inconvénient de cette méthode est bien sûr l'erreur qui provient d'un déplacement insuffisamment petit.

Pour définir la valeur de la densité de courant dans la bobine, on choisit un diamètre de fil d_f, une intensité de courant dans ce fil I_f, puis on calcule la densité de courant :

$$\begin{array}{ccc}\hfill & N=\frac{cd}{d_f^2}\hfill & \hfill \text{(18)}\\ \hfill & J=\frac{N{I}_f}{cd}=\frac{I_f}{d_f^2}\hfill & \hfill \text{(19)}\end{array}$$

Le terme de source qu'il faut placer dans l'équation de Poisson est $-{\mu }_{0}J$. Nous ferons le calcul numérique en système d'unités international car ce terme de source a un ordre de grandeur proche de 1 (pour un courant de 1 A dans un fil de diamètre 1 mm).

L'aimant est défini par son aimantation M. Le tableau suivant donne les valeurs typiques de champ rémanent, de champ coercitif de l'aimantation et de champ coercitif pour les quatre familles de ferromagnétiques :

Ferromagnétiques durs

Matériau	Br (T)	μ0Hcm (T)	μ0Hc (T)	(BH)max (kJ/m3)
AlNiCo	1.3	0.06	0.06	50
Ferrites	0.4	0.4	0.37	30
NdFeB	1.3	1.5	1.25	320
SmCo	1.1	1.3	0.97	220

Pour un aimant idéal, l'aimantation est donnée par ${\mu }_{0}M=B_r$.

Les calculs seront faits pour $\ell =1m$ : les grandeurs calculées sont donc par unité de longueur.

from matplotlib.pyplot import *
import numpy as np
from scipy.signal import convolve
import poisson.main
import gc
                

La fonction suivante effectue le calcul du potentiel vecteur, du champ magnétique, du flux dans la bobine, de la force sur l'aimant et de l'énergie du champ :

def calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc,yc,a,b,J,Br,noyau=False,extent=None,mur=1000,plaque=False,ep=2,lp=8):
    # n : taille du maillage fin
    # p : taille du maillage réduit
    # e,c,d : entier pairs, dimensions de la bobine
    # L : entier pair, longueur du noyau
    # xc,yc : entiers, position du centre de l'aimant
    # a,b : entiers pairs, dimensions de l'aimant 
    # J : densité de courant dans la bobine 
    # Br : champ B rémanent (mu0 * M)
    # extent = limites de la région intéressante
    print("xc = ",xc," yc = ",yc)
    levels = n-p+1
    echelle = Ld/2**p # échelle de la maille réduite
    Dx = Ld/2**n # taille de la maille du maillage fin
    ci,cj = 2**(p-1),2**(p-1) # centre du domaine
    systeme.init_array()
    systeme.laplacien()
    systeme.dirichlet_borders(0.0)
    if noyau:
        systeme.interface_polygon(ci-L//2,cj-e//2,[0,L,0,-L],[e,0,-e,0],0,1,1/mur)
    mu0 = 4*np.pi*1e-7
    systeme.source_rect(ci,cj+e//2+d//2,c//2,d//2,-J*mu0) # cuivre
    systeme.source_rect(ci,cj-e//2-d//2,c//2,d//2,J*mu0) # cuivre
    Jm = Br/(Dx*mu0)
    systeme.source_polygon(ci+xc-a//2,cj+yc+b//2,[a],[0],-Jm*mu0)
    systeme.source_polygon(ci+xc-a//2,cj+yc-b//2,[a],[0],Jm*mu0)
    systeme.source_polygon(ci+xc-a//2,cj+yc-b//2,[0],[b],0) # pour le dessin des bords
    systeme.source_polygon(ci+xc+a//2,cj+yc-b//2,[0],[b],0)
    if plaque:
        systeme.interface_polygon(ci+xc-a//2-ep,cj+yc-lp//2,[0,ep,0,-ep],[lp,0,-lp,0],0,1,1/mur)
    result=systeme.multigrid_full_norm(4,levels,200,2,2)
    #result = systeme.opencl_iterations_norm(200,100,omega=1.9)
    Az=systeme.get_array()
    Bx=systeme.get_derivY()
    By=-systeme.get_derivX()
    f = 2**(n-p)
    i1 = (ci-c//2)*f
    i2 = (ci+c//2)*f 
    j1 = (cj+e//2)*f 
    j2 = (cj+e//2+d)*f
    A1 = Az[j1:j2+1,i1:i2+1] # bobine +
    j1 = (cj-e//2-d)*f 
    j2 = (cj-e//2)*f
    A2 = Az[j1:j2+1,i1:i2+1] # bobine -
    flux = A1.mean()-A2.mean()
    Ec = 0.5*(A1-A2).sum()*J*Dx**2
    i1 = (ci+xc-a//2)*f
    i2 = (ci+xc+a//2)*f
    j1 = (cj+yc+b//2)*f
    j2 = (cj+yc-b//2)*f
    Fx = By[j2,i1:i2+1]-By[j1,i1:i2+1]
    Fx = Fx.sum()*Dx*Br/mu0
    Fy = Bx[j1,i1:i2+1]-Bx[j2,i1:i2+1]
    Fy = Fy.sum()*Dx*Br/mu0
    Ec += 0.5*(Az[j1,i1:i2+1]-Az[j2,i1:i2+1]).sum()*Jm*Dx**2
    mask = systeme.get_mask()
    if extent:
        i1 = 2**(n-1)+int(extent[0]/Dx)
        i2 = 2**(n-1)+int(extent[1]/Dx)
        j1 = 2**(n-1)+int(extent[2]/Dx)
        j2 = 2**(n-1)+int(extent[3]/Dx)
        return Az[j1:j2,i1:i2],Bx[j1:j2,i1:i2],By[j1:j2,i1:i2],flux,Fx,Fy,Ec,mask[j1:j2,i1:i2]
    else:
        return Az,Bx,By,flux,Fx,Fy,Ec,mask
        
                

La fonction suivante effectue le calcul de dérivation de l'énergie pour obtenir la composante de la force dans la direction du déplacement :

def gradientEc(xc,Ec,Dx):
    F = np.zeros(len(Ec)-2)
    x = np.zeros(len(Ec)-2)
    for k in range(1,len(Ec)-1):
        F[k-1] = (Ec[k+1]-Ec[k-1])
        x[k-1] = xc[k]
    F /= (2*Dx)
    return x,F
                

Voici la définition du système et un premier calcul :

n=10
p=8
Ld = 1.0
e,c,d = 16,16,4
L = 16
xc,yc = -20,0
a,b = 4,16
df = 1e-3
If = 5.0
J = If/df**2
Br = 0.4 # en teslas
extent = [-0.2,0.2,-0.2,0.2]

systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
Az,Bx,By,flux,Fx,Fy,Ec,mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc,yc,a,b,J,Br,noyau=True,extent=extent)
systeme.close()
figure(figsize=(8,8))
#extent=systeme.get_extent()
contour(Az,30,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('x')
ylabel('y')
grid()
                

fig1.pdf

print(flux)
--> 0.0068979273

print(Fx)
--> 131.30994211972106

print(Fy)
--> 0.12643557643274472

3.a. Déplacement longitudinal

L'aimant est déplacé sur l'axe de la bobine.

xc = np.arange(-30,-10)
yc=0
flux = np.zeros(len(xc))
Fx1 = np.zeros(len(xc))
Fy1 = np.zeros(len(xc))
Ec = np.zeros(len(xc))
for k in range(len(xc)):
    systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
    Az,Bx,By,flux[k],Fx1[k],Fy1[k],Ec[k],mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc[k],yc,a,b,J,Br,noyau=True,extent=extent)
    systeme.close()
del Az
del Bx
del By
del mask
gc.collect()
figure(figsize=(12,6))
plot(xc,flux)
ylim(0,flux.max())
grid()
xlabel('xc',fontsize=16)
ylabel('flux (Wb/m)',fontsize=16)
                

fig2.pdf

Voici les deux composantes de la force sur l'aimant, obtenues par calcul direct :

figure(figsize=(12,6))
plot(xc,Fx1,label='Fx')
plot(xc,Fy1,label='Fy')
grid()
legend(loc='upper left',fontsize=16)
xlabel('xc',fontsize=16)
ylabel('Force (N/m)',fontsize=16)
#ylim(0,600)
                

fig3.pdf

Voici l'énergie du champ $E_c$ en fonction de la position de l'aimant :

figure(figsize=(12,6))
plot(xc,Ec)
grid()
xlabel('xc',fontsize=16)
ylabel('Ec (J/m)',fontsize=16)
                

fig4.pdf

La composante F_x de la force sur l'aimant s'en déduit par dérivation par différence finie.

Dx = Ld/2**p
x1,F1 = gradientEc(xc,Ec,Dx)
figure(figsize=(12,6))
plot(x1,F1,'b--',label='dérivation de Ec')
plot(xc,Fx1,'b-',label='calcul direct')
grid()
legend(loc='upper left')
xlabel('xc',fontsize=16)
ylabel('Fx (N/m)',fontsize=16)
#ylim(0,600)
                

fig5.pdf

Afin de déterminer l'influence du noyau, on refait le calcul sans le noyau :

xc = np.arange(-30,-10)
yc=0
flux = np.zeros(len(xc))
Fx2 = np.zeros(len(xc))
Fy2 = np.zeros(len(xc))
Ec = np.zeros(len(xc))
for k in range(len(xc)):
    systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
    Az,Bx,By,flux[k],Fx2[k],Fy2[k],Ec[k],mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc[k],yc,a,b,J,Br,noyau=False,extent=extent)
    systeme.close()
x2,F2 = gradientEc(xc,Ec,Dx)
del Az
del Bx
del By
del mask
gc.collect()
figure(figsize=(12,6))
plot(xc,flux)
ylim(0,flux.max())
grid()
xlabel('xc',fontsize=16)
ylabel('flux (Wb/m)',fontsize=16)
                

fig6.pdf

figure(figsize=(12,6))
plot(xc,Fx2,label='Fx')
plot(xc,Fy2,label='Fy')
grid()
legend(loc='upper left',fontsize=16)
xlabel('xc',fontsize=16)
ylabel('Force (N/m)',fontsize=16)
                

fig7.pdf

On peut aussi calculer l'effet du noyau seul (sans la bobine) :

xc = np.arange(-30,-10)
yc=0
flux = np.zeros(len(xc))
Fx3 = np.zeros(len(xc))
Fy3 = np.zeros(len(xc))
Ec = np.zeros(len(xc))
for k in range(len(xc)):
    systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
    Az,Bx,By,flux[k],Fx3[k],Fy3[k],Ec[k],mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc[k],yc,a,b,0,Br,noyau=True,extent=extent)
    systeme.close()
x3,F3 = gradientEc(xc,Ec,Dx)
del Az
del Bx
del By
del mask
gc.collect()
figure(figsize=(12,6))
plot(xc,flux)
ylim(0,flux.max())
grid()
xlabel('xc',fontsize=16)
ylabel('flux (Wb/m)',fontsize=16)
                

fig8.pdf

figure(figsize=(12,6))
plot(xc,Fx3,label='Fx')
plot(xc,Fy3,label='Fy')
grid()
legend(loc='upper left',fontsize=16)
xlabel('xc',fontsize=16)
ylabel('Force (N/m)',fontsize=16)
                

fig9.pdf

Voici le tracé de la force F_x dans les trois cas. Pour chaque cas, la courbe en trait plein est la force sur l'aimant obtenue par calcul direct, la courbe en trait interrompu est la force obtenue par dérivation de l'énergie du champ.

figure(figsize=(12,6))
plot(xc,Fx1,'r-',label='bobine + noyau')
plot(x1,F1,'r--',label='bobine + noyau')
plot(xc,Fx2,'g-',label='bobine')
plot(x2,F2,'g--',label='bobine')
plot(xc,Fx3,'b-',label='noyau')
plot(x3,F3,'b--',label='noyau')
grid()
ylim(0,600)
legend(loc="upper left",fontsize=16)
legend(loc='upper left',fontsize=16)
xlabel('xc',fontsize=16)
ylabel('Fx (N/m)',fontsize=16)
                

fig10.pdf

Lorsque seul le noyau est présent, la force attractive est due à l'aimantation de celui-ci par le champ de l'aimant. Lorsque seule la bobine est présente, la force attractive est due aux courants dans celle-ci. Lorsque la bobine et le noyau sont présents, l'aimantation du noyau est plus forte qu'en l'absence de bobine; c'est pourquoi la force est plus grande que la somme des forces des deux cas précédents.

3.b. Déplacement transversal

L'aimant est déplacé perpendiculairement à l'axe de la bobine.

xc = -12
yc = np.arange(0,15)
flux = np.zeros(len(yc))
Fx1 = np.zeros(len(yc))
Fy1 = np.zeros(len(yc))
Ec = np.zeros(len(yc))
for k in range(len(yc)):
    systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
    Az,Bx,By,flux[k],Fx1[k],Fy1[k],Ec[k],mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc,yc[k],a,b,J,Br,noyau=True,extent=extent)
    systeme.close()
y1,F1 = gradientEc(yc,Ec,Dx)
figure(figsize=(8,8))
#extent=systeme.get_extent()
contour(Az,30,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('x')
ylabel('y')
del Az
del Bx
del By
del mask
gc.collect()

fig11.pdf

Voici les deux composantes de la force sur l'aimant. La composante dans la direction du déplacement est aussi calculée par dérivation de l'énergie du champ (courbe en trait interrompu) :

figure(figsize=(12,6))
plot(yc,Fx1,'r-',label='Fx')
plot(yc,Fy1,'b-',label='Fy')
plot(y1,F1,'b--',label='Fy')
grid()
legend(loc='upper left',fontsize=16)
xlabel('yc',fontsize=16) 
ylabel('Force (N/m)',fontsize=16)
                

fig12.pdf

La composante transversale de la force (F_y) est importante dans le cas d'un moteur dont le rotor comporte des aimants permanents. Pour maximiser cette force, il faut maximiser la composante B_x du champ (qui dans un moteur correspond à la composante radiale). En plaçant une plaque en fer doux derrière l'aimant, on devrait augmenter cette composante du champ :

xc = -12
yc2 = np.arange(0,15)
flux2 = np.zeros(len(yc))
Fx2 = np.zeros(len(yc))
Fy2 = np.zeros(len(yc))
Ec2 = np.zeros(len(yc))
for k in range(len(yc)):
    systeme = poisson.main.Poisson(p,p,n-p+1,Ld,Ld,x0=-Ld/2,y0=-Ld/2)
    Az,Bx,By,flux2[k],Fx2[k],Fy2[k],Ec2[k],mask = calculPotentiel(systeme,n,p,Ld,e,c,d,L,xc,yc[k],a,b,J,Br,noyau=True,extent=extent,plaque=True,ep=4,lp=32)
    systeme.close()
figure(figsize=(8,8))
#extent=systeme.get_extent()
contour(Az,30,extent=extent)
imshow(~mask,extent=extent,alpha=0.5, cmap=cm.gray)
xlabel('x')
ylabel('y')
del Az
del Bx
del By
del mask
gc.collect()


            

fig13.pdf

figure(figsize=(12,6))
plot(yc,Fy1,label='Sans plaque')
plot(yc2,Fy2,label='Avec plaque')
grid()
legend(loc='upper right',fontsize=16)
xlabel('yc',fontsize=16)
ylabel('Fy (N/m)',fontsize=16)
ylim(-600,0)
                

fig14.pdf

L'ajout d'une plaque en fer derrière l'aimant a bien un effet très important sur la force transversale. Il y a de plus une force sur la plaque en fer, qui n'est pas prise en compte dans le calcul. Dans un rotor de moteur, les forces sur les parties en tôle de fer contribuent au couple moteur.

Voici l'énergie pour différentes positions du bloc aimant-plaque :

figure(figsize=(12,6))
plot(yc2,Ec2)
grid()
xlabel('yc',fontsize=16)
ylabel('Ec (J/m)',fontsize=16)

                

fig15.pdf

Par dérivation, on obtient la force sur le bloc aimant-plaque :

Dx = Ld/2**p
ycc,Fyy = gradientEc(yc2,Ec2,Dx)
figure(figsize=(12,6))
plot(ycc,Fyy)
grid()
xlabel('yc',fontsize=16)
ylabel('Fy (N/m)',fontsize=16)
ylim(-600,0)
                

fig16.pdf