Resolution numerique EDP - FEM

[Nicolas666666]

Bonjour a tous,

Voila un moment que je n'avais pas poste, cette fois ce sera pour demander un petit peu d'aide !

Dans le cadre de mon stage je dois entre autre resoudre des EDP, et je commence par une simple : avec h dependant donc uniquement de x (1D) et de t.

Je dois utiliser les elements finis, je passe donc a la weak formulation ou je pose les weight functions egales aux shapes functions (galerkin method)

J'arrive donc a une equation matricielle reliant ma matrice de force, ma matrice h discretise sur l'espace a l'intant t+1 et ma matrice h a l'instant t.

Pour mes conditions initiales je peux imaginer un truc folklo, du type h(x,0)=constante.

La ou je bloque c'est sur mes conditions aux limites. J'ai au moins sur un cote quelque chose du genre et sur l'autre... disons la meme chose.
Mais je ne vois pas du tout comment integrer ceci a mon systeme matriciel et surtout comment le resoudre.

Je sais que ca risque d'etre difficile de m'aider avec si peu d'information et sans mes calculs sous les yeux, mais qui ne tente rien n'a rien et vous avez peut-etre des exemples types sous la main ?

Pour la petite histoire je suis etudiant en Genie Civil et je fais un stage sur matlab base sur les elements finis, autant dire pas trop mon domaine

Merci beaucoup !
-----

[KerLannais]

Bonjour,

Les conditions aux bord sont contenues dans ta formulation variationnelle et donc dans ton système matriciel normalement, et donc tu n'as pas à t'en préoccuper.

Mettons que l'on prenne comme domaine d'espace un segment , le temps étant naturellement pris dans . Soient des constantes. On cherche la solution du problème




Remarques:Il faut prendre pour avoir une équation du type équation de la chaleur, sinon on a l'équation de la chaleur avec le temps renversé et au lieu de se diffuser la solution se concentre en certains points et devient très singulière
En général est strictement positif mais la condition en est


Prenons une solution classique (i.e. une fonction continue sur une fois dérivable en temps et deux fois en espace et qui vérifie le système d'équation au sens classique). On suppose également que et sont continues sur . Prenons une fonction test de classe sur (attention j'ai pas dit à support compact) et continue sur .

On a

On peut sortir l'intégration en temps de la première intégrale et on fait une intégration par parties dans la seconde et on trouve


Si maintenant tu prends une famille de telles que si alors

(ce qui est en général le cas pour les éléments finis classiques P1, P2 ...) et si tu suppose

alors pour et pour tout


Soit

avec

une matrice qui est inversible dès que la famille est une base.

où

et enfin


Il faut aussi rajouter la condition initiale qui donne une condition sur les . On peut alors résoudre, pour assez grand on obtient une solution qui vérifie de façon approchée la formulation variationnelle et ont peut même (lorsque la forme bilinéaire est coercive) quantifier la distance de cette solution approchée à la solution de la formulation variationnelle qui vérifie au sens faible l'équation et les conditions au bords (cela découle directement de la formulation variationnelle) et si on montre qu'en fait la solution est plus régulière alors elle est solution au sens classique et vérifie les conditions au bord au sens classique.

Ce que les gens ne comprennent pas en général c'est que les conditions aux bords sont déjà incluse dans la formulation variationnelle et qu'il ne faut pas essayer de rajouter des équations en plus au système matriciel pour les vérifier. L'exercice le plus simple à faire pour s'en convaincre est de montrer en exercice que si est une fonction régulière qui vérifie la formulation variationnelle alors elle est solution du problème avec les conditions au bord (indication: commencer par considérer des fonctions test qui s'annulent en et en , montrer que est solution de l'équation puis avec des fonctions test qui valent soit en et en , soit le contraire, vérifier que vérifie les conditions au bord en et en ).

[Nicolas666666]

Bonjour !

Toujours aussi impressionne par ce forum...

Je suis en train d'essayer de comprendre toutes tes explications, mais j'ai quand meme mis en lien une partie de ce que j'avais fait, si jamais tu as 5 minutes pour y jeter un oeil...

Encore merci, et surtout a bientot

Nicolas

[Nicolas666666]

Bonjour,

Les conditions aux bord sont contenues dans ta formulation variationnelle et donc dans ton système matriciel normalement, et donc tu n'as pas à t'en préoccuper.

Mettons que l'on prenne comme domaine d'espace un segment , le temps étant naturellement pris dans . Soient des constantes. On cherche la solution du problème




Remarques:Il faut prendre pour avoir une équation du type équation de la chaleur, sinon on a l'équation de la chaleur avec le temps renversé et au lieu de se diffuser la solution se concentre en certains points et devient très singulière
En général est strictement positif mais la condition en est


Prenons une solution classique (i.e. une fonction continue sur une fois dérivable en temps et deux fois en espace et qui vérifie le système d'équation au sens classique). On suppose également que et sont continues sur . Prenons une fonction test de classe sur (attention j'ai pas dit à support compact) et continue sur .

On a

On peut sortir l'intégration en temps de la première intégrale et on fait une intégration par parties dans la seconde et on trouve

Jusque la tout va bien !

Si maintenant tu prends une famille de telles que si alors

(ce qui est en général le cas pour les éléments finis classiques P1, P2 ...) et si tu suppose

alors pour et pour tout

La ca va beaucoup moins bien... C'est le principe de la formule variationnelle de poser h(x,t) fonction de deux fonctions dependant respectivement du temps et de l'espace ?

Quand tu utilises la notation c'est pour symboliser la construction matricielle du decoupage en elements de la zone en question ?

Soit

avec

une matrice qui est inversible dès que la famille est une base.

où

et enfin


Il faut aussi rajouter la condition initiale qui donne une condition sur les . On peut alors résoudre, pour assez grand on obtient une solution qui vérifie de façon approchée la formulation variationnelle et ont peut même (lorsque la forme bilinéaire est coercive) quantifier la distance de cette solution approchée à la solution de la formulation variationnelle qui vérifie au sens faible l'équation et les conditions au bords (cela découle directement de la formulation variationnelle) et si on montre qu'en fait la solution est plus régulière alors elle est solution au sens classique et vérifie les conditions au bord au sens classique.

Ce que les gens ne comprennent pas en général c'est que les conditions aux bords sont déjà incluse dans la formulation variationnelle et qu'il ne faut pas essayer de rajouter des équations en plus au système matriciel pour les vérifier.

Et la ca va mieux !

Je ne suis pas du tout au point sur toute cette partie mathematique, et j'espere que la question suivante n'est pas trop stupide :

En se placant en regime permanent, soit les conditions aux limites telles que decrites plus haut ne nous permettent pas de resoudre , je me trompe ?
Comment est-ce possible ?

Encore merci pour ton aide !!!

[A voir en vidéo sur Futura]

[KerLannais]

Pour l'équation stationnaire tu résous juste BC+D=0

j'ai écris un petit programme sous scilab (qui est l'analogue de matlab mais en gratuit) la syntaxe est vraiment très semblable à celle de matlab , je pense donc que tu pourras comprendre

function erreur=Chaleur(henv,L,D,T,N,delta_t)
  
//On veut résoudre l'équation de la chaleur sans terme source, en une
//dimension d'espace avec des conditions aux bords qui modélisent de
//la convection (ce sont alors des conditions aux bords de type Robin
//-Fourier.

//On travaille sur un segment [0,L]

//On veut résoudre l'équation
//d_t u(t,x)-D*d_x^2 u(t,x)=0      (1)

//avec D le coefficient de diffusion strictement positif

//Pour que l'équation ait un sens classique il faut que la fonction u
//solution soit une fonction de [0,+infini[*[0,L] dans R qui soit
//dérivable par rapport à sa première variable et deux fois dérivable
//par rapport à sa deuxième variable. On supposera en plus que toutes
//ses dérivées partielles sont continues sur [0,+\infini[*[0,L] pour
//pouvoir les intégrer tranquillement. On a noté d_t la dérivée par
//rapport à la première variable de u et d_x^2 la dérivée seconde par
//rapport à la deuxième variable de u.

//On aura la condition initiale suivante (choisie exprès)
//u(0,x)=u0(x)=henv+sin(x)

//avec henv une constante donnant la température extérieure, on peut
//toujours se ramener à ce que cette température soit hext=0 quitte à
//rajouter henv à la solution finale trouvée. Ainsi on pose

function y=u0(x)
  y=sin(x);
endfunction

//On pose une condition de Dirichlet homogène en x=0 qui vaut henv,
//autrement dit l'extrémité en x=0 est à température constante égale à
//la température extérieure soit u(t,0)=henv et donc on suppose ici
//u(t,0)=0.

//On pose une condition de type convection en x=L, avec b un
//coefficient de convection que l'on définit ici volontairement par
//rapport à L

b=-cos(L)/sin(L);

//On veut imposer d_x u(t,L)=-b*(u(t,L)-henv). On imposera donc
//d_x u(t,L)=-b*u(t,L).

//On va utiliser la méthode de Galerkin en espace ce qui impose
//d'utiliser des fonctions seulement H^1 et pas forcément dérivables
//au sens classique. On va donc utiliser une formulation
//variationnelle de (1) qui aura un sens pour des fonctions u
//seulement H^1 et qui sera équivalent pour des fonctions u qui ont
//la régularité susdite. Pour cela on suppose que u est régulière et
//solution de (1), on multiplie par une fonction test v de classe C^1
//sur [0,L] et à valeurs dans R et on intègre sur le segment [0,L]
//(on notera int_0^L cette intégrale). On supposera que la fonction
//test vaut 0 en x=0.

//int_0^L d_t u(t,x)v(x)dx-D*int_0^L d_x^2u(t,x)v(x)dx=0

//On fait alors une intégration par partie sur la deuxième intégrale

//int_0^L d_t u(t,x)v(x)dx-D*d_x u(t,L)v(L)
//+D*int_0^L d_x u(t,x)*d_x v(x)dx=0

//On remplace la condition au bord en x=L

//int_0^L d_t u(t,x)v(x)dx+D*b*u(t,L)v(L)
//+D*int_0^L d_x u(t,x)*d_x v(x)dx=0

//Cette formulation a alors un sens pour u H^1 en x qui s'annule en
//x=0 et dérivable en t avec une dérivée en t continue. Si cette
//formulation variationnelle est vérifiée pour tout v H^1 qui
//s'annule en a et si u est régulière alors u est solution classique.
//La méthode de Galerkin consiste à prendre des familles libres
//v_(0,N), v_(1,N), ... , v_(N,N) dans l'ensemble des fonctions de
//H^1(]0,L[) qui s'annulent en x=0 et telles que pour toute fonction
//H^1 qui sannule en x=0 s'approche bien pour N suffisamment grand
//par une combinaison linéaire de cette famille. On prendra les
//éléments P1, il s'agit alors de prendre une subdivision de N+1
//points de l'intervalle [0,L]

x=linspace(0,L,N+1);

//et de considérer les fonctions continues v(i,.), affines par
//morceaux sur cette subdivision qui valent 1 sur un x(i) avec i>0 et
//0 sur les autres points de la subdivision. On cherche alors u(t,x)
//sous la forme u(t,x)=c_1(t)*v(1,x)+...+c_N(t)v(N,x)
//et on impose que u vérifie la formulation variationnelle seulement
//pour des fonctions test dans cet espace de dimension fini (et donc
//en fait seulement pour les N fonctions test v(i,.)). On obtiendra
//alors une suite de solutions approchées pour chaque valeur de N qui
//converge vers la solution en norme H^1 uniformément en temps. Quand
//on remplace on trouve pour 1<=k<=N
//c_1'(t)'int_0^L v(1,x)*v(k,x)dx+...+c_N'(t)int_0^L v(N,x)v(k,x)dx
//+D*b*c_N(t)*delta(k=N)+D*c_1(t)int_0^L d_x v(1,x)*d_x v(k,x) dx+...
//+D*c_N(t)int_0^L d_x v(N,x)*d_x v(k,x) dx=0
//avec delta(k=N)=1 si k=N et 0 sinon.

//Soit, en notant A la matrice (symétrique) dont les coeffients sont
//les int_0^L v(i,x)v(k,x)dx, C le vecteur vertical de composantes
//c_i(t), B la matrice de composante
//D*int_0^L d_x v(i,x)*d_x v(k,x)dx+D*b*delta(i=k=N)
//AC'+BC=0
//Il est facile de calculer explicitement la matrice A et la matrice
//B et il est clair que A est inversible dès que les v(i,.) sont
//libres ce qui est évidemment le cas.

dA=2*L/(3*N)*ones(1,N);
dA(N)=L/(3*N);
sA=L/(6*N)*ones(1,N-1);
//A=diag(dA,0)+diag(sA,1)+diag(sA,-1);

dB=2*D*N/L*ones(1,N);
dB(N)=D*N/L+D*b;
sB=-D*N/L*ones(1,N-1);
B=diag(dB,0)+diag(sB,1)+diag(sB,-1);

//La matrice A étant symétrique on peut la mettre sous la forme de
//Cholesky et vu sa forme cette décomposition est très simple elle
//donne lieu à la matrice L triangulaire inférieure avec seulement la
//diagonale et la sous-diagonale qui est non nulle:

dL=ones(1,N);
sL=ones(1,N-1);
dL(1)=sqrt(dA(1));
for i=1:(N-1)
    sL(i)=sA(i)/dL(i);
    dL(i+1)=sqrt(dA(i+1)-sL(i)^2);
end

//L=diag(dL,0)+diag(sL,-1);
//Il est facile de calculer l'inverse de Linv de L puisqu'il s'agit
//de résoudre un système triangulaire.

Linv=zeros(N,N);
for n=1:N
    if n==1
        Linv(1,n)=1/dL(1);
    else
        Linv(1,n)=0;
    end
    for m=2:N
        if m==n
            Linv(m,n)=(1-sL(m-1)*Linv(m-1,n))/dL(m);
        else
            Linv(m,n)=-sL(m-1)*Linv(m-1,n)/dL(m);
        end
    end
end

//Si on note F=LC et E=(Linv)^tBLinv (où ^t est la transposée) alors
//F'+EF=0
//On se ramène donc à la résolution d'un système d'EDO linéaire
//d'ordre 1, puis on calcule C=LinvF.

E=Linv'*B*Linv;

//Résolvons cette EDO à l'aide de la méthode D'Euler explicite. On
//prend un pas Delta_t petit. On va résoudre sur un intervalle [0,T].
//L'ensemble des valeurs de temps où l'on calcule la solution est
//alors
t=0:delta_t:T;

F=zeros(N,length(t));
//Initialement on a C(0)=u0(x) et donc

F(:,1)=(u0(x(2:N+1)).*dL+u0(x(1:N)).*[0 sL])';
for i=1:(length(t)-1)
  F(:,i+1)=F(:,i)-delta_t*E*F(:,i);
end

C=Linv*F+henv;

//Sous ces hypothèses la solution exacte est connue, c'est

Ce=zeros(N,length(t));
for i=1:N;
  Ce(i,:)=henv+exp(-D*t)*sin(x(i+1));
end

s=int(linspace(1,length(t),100));

figure(0)
clf
surf(t(s),x(2:N+1),C(:,s))

figure(1)
clf
surf(t(s),x(2:N+1),Ce(:,s))

erreur=max(abs(C-Ce));

endfunction

les lignes commençant par // sont des commentaires
pour tester il ne faut pas oublier que la condition CFL (de stabilité) dit qu'il faut que le pas de temps soit de l'ordre du carré du pas d'espace, par exemple j'ai pris

et

Par ailleurs j'avais pris




en prenant comme condition initiale

on a alors dans ce cas une solution exacte

j'obtiens une erreur de entre la solution numérique que je calcule et la solution exacte ce qui est raisonnable.

[Nicolas666666]

Désolé de ne pas avoir répondu plus tôt, juste un grand merci !
J'ai mis du temps a comprendre et a programmer mais ça fonctionne !

Encore merci pour ton temps et ton aide...

--
Nicolas