Modélisation numérique sous MATLAB
Par Rafic YOUNES
Enseignant - Chercheur à la Faculté de Génie – Université Libanaise Responsable du DEA Mécanique E2M
R.Y. 02.04
1. Grandes lignes de la méthode des éléments finis
1.1 Exposé de la démarche
1.2 Un exemple détaillé
2. Utilisation du pdetool
3. Applications
3.1 Application 1
3.2 Application 2
3.3 Application 3
3.4 Application 4
3.5 Application 5
Quelle place occupe le calcul dans l'industrie ? Quels sont les principaux champs d'application du calcul ?
Le calcul numérique permet à l'ingénieur d'effectuer des simulations numériques de phénomènes physiques. Le calcul occupe une place stratégique avec la CAO et les autres technologies de simulation (essais) dans le développement d'un produit complexe qui touche à différents domaines de la physique. Cela concerne les industries automobiles, navales, aéronautiques, ferroviaires, mais aussi les industries lourdes: centrales électriques, platesformes pétrolières, et le génie civil.
Le calcul est indispensable lorsque l'on cherche à obtenir une solution optimisée pour réduire les coûts et les délais de fabrication. Grâce au calcul, même simplifié, l'ingénieur peut tester plusieurs configurations pour optimiser le comportement d'un modèle à une prestation donnée. Cela évite de multiplier les prototypes et les essais tests réels, les supports physiques ne servent plus à chercher une solution, ils permettent de la valider.
Le calcul s'applique aussi dans les domaines du « process ». Les procédés de fabrication tels que l'emboutissage, l'usinage grande vitesse, les dépôts de peinture, l'assemblage de tôlerie, la mise en forme des plastiques, peuvent être modélisés par éléments finis. Ici c'est une bonne représentation du comportement du phénomène physique qui sera recherchée pour pouvoir vérifier et valider un procédé de fabrication d'une pièce.
Enfin le calcul de conception dans les bureaux d'études, c'est sans doute le plus répandu car grâce aux outils de calcul simplifié dont disposent les logiciels de CAO modernes, la simulation numérique fait partie des outils de conception pour obtenir un comportement défini à priori qui détermine le dimensionnement, donc le dessin, des pièces mécaniques.
La méthode des éléments finis est de toutes les méthodes de discrétisation la plus utilisée car :
A) elle peut traiter des problèmes de géométrie complexe. B) elle couvre de nombreux domaines de la physique. C) les moyens informatiques actuels (puissance des calculateurs, outils de visualisation) la rende facile de mise en oeuvre. D) de nombreux logiciels généraux ou dédiés sont disponibles sur le marché, le tableau ci-dessous vous permettra de visiter les sites Internet proposés par les principaux fournisseurs de logiciel éléments finis.
Pdetool de MatLab, Abaqus, I-deas, Adina, Marc, Algor, Radioss, Ansys, Samcef, Catia, Systus
L'objectif de ce cours est de vous donner les bases essentielles sur lesquelles repose la méthode des éléments finis. Ensuite sont présentés quelques applications utilisant le pdetool de MatLab.
Beyrouth – 01/03/2005
1 Exposé de la démarche :
La méthode consiste à rechercher une solution approchée de la solution exacte sous la forme d'un champ F(M,t) défini par morceaux sur des sous domaines de ?. Les n sous domaines ?i doivent être tels que
où ?~i désigne l'intérieur de ?i. Autrement dit, les ?i sont une partition de .
~
Les champs f (M,t), définis sur chaque sous domaines sont des champs choisis parmi une famille arbitraire de champs (généralement polynômiaux). La famille de champs locaux est appelée espace des fonctions d'interpolation de l'élément. ~
La famille de champs globaux F(M,t), obtenus par juxtaposition des champs locaux est appelée espace des fonctions d'interpolation du domaine.
Le champ dans chaque sous domaine ?i est déterminé par un nombre fini de valeurs du champ (ou de valeurs de ses dérivées) en des points choisis arbitrairement dans le sous domaine, et appelés Q°uds. Le champ local est une interpolation entre les valeurs aux noeuds. Le sous domaine muni de son interpolation est appelé élément.
Chercher une solution par éléments finis consiste donc à déterminer quel champ local on attribue à chaque sous domaine pour que le champ global ~
F(M,t) obtenu par juxtaposition de ces champs locaux soit proche de la solution du problème.
Parmi les contraintes qu'on impose à la solution approchée cherchée, il y a souvent au moins une continuité simple (C0) à la frontière entre les sous domaines.
La figure 1 montre une solution approchée discontinue d'un champ scalaire sur un domaine de dimension 1. La famille de champs locaux est la famille des champs constants par morceaux.
Figure 1:Solution approchée discontinue
La figure 2 montre une solution approchée continue C0 d'un champ scalaire sur un domaine de dimension 1. La famille de champs locaux est la famille des champs polynômiaux de degré 1.
La figure 3 montre une solution approchée continue C1 d'un champ scalaire sur un domaine de dimension 1. La famille de champs locaux est la famille des champs polynômiaux de degré 3.
Figure 2:Solution approchée continue C0
Figure 3:Solution approchée continue C1
La qualité de la solution approchée dépend de la division en sous domaines (nombre et dimensions des sous domaines), du choix de la famille de champs locaux dans chaque sous domaine, et des conditions de continuité qu'on impose aux frontières des sous domaines (C0, C1,...). Une fois ces choix faits, il reste à rechercher, une combinaison de champs locaux qui satisfait approximativement les équations.
Pour résoudre un problème par la méthode des éléments finis, on procède donc par étapes successives :
1 - On se pose un problème physique sous la formed'une équation différentielle ou aux dérivés partielles à satisfaire en tout point d'un domaine , avec des conditions aux limites sur le bord nécessaires et suffisantes pour que la solution soit unique.
2 - On construit une formulation intégrale du systèmedifférentiel à résoudre et de ses conditions aux limites: C'est la formulation variationnelle du problème.
3 - On divise en sous domaines : C'est le maillage. Les sous domaines sont appelés mailles.
4 - On choisit la famille de champs locaux, c'est à direà la fois la position des noeuds dans les sous domaines et les polynômes (ou autres fonctions) qui définissent le champ local en fonction des valeurs aux noeuds (et éventuellement des dérivées). La maille complétée par ces informations est alors appelée élément.
5 - On ramène le problème à un problème discret : C'est la discrétisation. En effet, toute solution approchée est complètement déterminée par les valeurs aux noeuds des éléments. Il suffit donc de trouver les valeurs à attribuer aux noeuds pour décrire une solution approchée. Le problème fondamental de la méthode des éléments finis peut se résumer en deux questions :
(a) Comment choisir le problème discret dont la solution est <<proche>> de la solution exacte? (b) Quelle signification donner au mot <<proche>> ?
6 - On résout le problème discret: C'est la résolution
7 - On peut alors construire la solution approchée àpartir des valeurs trouvées aux noeuds et en déduire d'autres grandeurs : C'est le post-traitement.
8 - On visualise et on exploite la solution pour juger desa qualité numérique et juger si elle satisfait les critères du cahier des charges : C'est l'exploitation des résultats.
Les étapes 1, 2, 3, 4 et 5 sont souvent rassemblées sous le nom de prétraitement.
Le travail de ces différentes étapes est assisté par les logiciels. Pour maîtriser leur utilisation, il est indispensable de comprendre les fondements de la méthode, notamment les phases 3 et 4, ne serait-ce que pour comprendre et choisir intelligemment parmi les options qu'ils proposent.
2 UN EXEMPLE DETAILLE:
On se propose de rechercher une solution approchée du problème suivant :
Trouver f(x) dans le domaine ? = [0, 1] satisfaisant l'équation différentielle :
:
est un domaine de dimension 1
et le temps n'apparaît pas. Pour décrire le domaine, on n'a besoin que d'une seule variable x. L'équation à résoudre est donc une équation différentielle ordinaire. La frontière se réduit à deux points.
La solution exacte de ce problème est
Son graphe est donné figure 6 page 8.
On divise arbitrairement en trois mailles de même taille (voir figure 4):
Figure 4:Maillage du problème
On décide de prendre des éléments à 3 noeuds, et pour la famille de champs locaux des polynômes de degré 2. Les noeuds sont choisis aux extrémités et au milieu de chaque maille.
On peut alors déterminer chaque champ local en fonction des valeurs aux 3 noeuds.
Remarquer que le fait d'avoir utilisé des noeuds aux extrémités de chaque élément présente deux avantages :
1 - Le nombre de noeuds est réduit, car il y a desnoeuds communs à deux éléments.
2 - On assure ainsi une continuité C0 de la solution approchée : les champs locaux de deux éléments voisins auront la même valeur à leur nud commun.
Remarquer encore qu'il n'est pas nécessaire de prendre les éléments identiques : on aurait pu prendre certains à deux noeuds et d'autres à trois noeuds.
Chaque champ local peut donc s'exprimer en fonction des valeurs aux noeuds. En effet, il n'existe qu'un seul polynôme du second degré qui satisfasse les conditions suivantes :
où x1, x2 et x3 sont les coordonnées des 3 noeuds de l'élément, et f1, f2 et f3 sont des valeurs aux 3 noeuds de l'élément.
On obtient :
Le polynôme avec les coefficients a1, b1 et c1 cidessus est appelé fonction d'interpolation de l'élément 1.
Les polynômes d'interpolation des 2 autres sont de la même forme : Pour le second, on remplace les indices 1, 2, 3 par les indices 3, 4, 5 ; et pour le troisième, on remplace les indices 1, 2, 3 par les indices 5, 6, 7.
Dans cet exemple, on a 7 noeuds. Si on donne une valeur arbitraire à chaque noeud, les fonctions d'interpolation sont déterminées et définissent par
~
morceaux une fonction f continue C0 sur le domaine. L'espace des fonctions , définies en 3 morceaux, est de dimension 7.
Les trois fonctions d'interpolation sur les 3 éléments sont :
On montre en analyse fonctionnelle que
conditions de régularité qu'on n'explicitera pas ici.
Autrement dit, résoudre l'équation différentielle est équivalent à chercher f(x) tel que :
Cette formulation est une formulation variationnelle (on dit aussi formulation intégrale ou encore formulation faible) du problème.
On obtient une autre formulation variationnelle en utilisant l'intégration par parties :
Compte tenu des conditions aux limites, f(0)=0 et ???dfdx(x)? = 0 , une autre formulation variationnelle
??x=1
du problème est donc : Trouver f(x) tel que
Il est important de noter que dans les termes de bord de cette dernière formulation, les conditions aux limites sont prises en compte. En comparant les deux formulations variationnelles (2.7) et (2.9), on constate que, dans la seconde forme, la fonction inconnue f(x)intervient avec des dérivées d'ordre inférieur, et que par contre les dérivées des fonctions ?(x) apparaissent. L'équivalence entre ces deux formulations variationnelles et le problème initial est soumise à des conditions sur f(x), ?(x) et précisées par l'analyse fonctionnelle et qu'on supposera satisfaites.
Remarque :
La démarche suivie ici pour construire une formulation variationnelle du problème est purement mathématique. Il arrive souvent que la formulation variationnelle ait une interprétation physique. On peut donc aussi les établir par des raisonnements physiques.
Le principe de l'approximation par éléments finis est le suivant :
On choisit une formulation variationnelle (par exemple (2.9) et on cherche à satisfaire l'équation (2.9) avec
~
les f déterminés par les valeurs aux noeuds définis précédemment. Il va de soit que ce problème, n'a généralement pas de solution si on garde la condition
~
<<? ?i(x) >> (sinon f serait une solution exacte !). On ne vérifiera donc cette équation que pour certains ?i(x)seulement.
À chaque ?i(x) correspondra une équation scalaire fonction des valeurs aux noeuds. Il suffit de les choisir indépendants et en nombre égal au nombre d'inconnues. Dans notre exemple, il en faut 7 pour déterminer les valeurs aux noeuds f1,...,f7. Les fonctions ?(x) sont appelées fonctions test ou encore fonctions de pondération. On est ainsi ramené à un système algébrique de 7 équations à 7 inconnues.
Il y a une infinité de manières de choisir les 7 ?i(x). Elles engendrent différentes variantes de la méthode des éléments finis. Chacune aboutit à une solution approchée différente. La seule condition impérative est que le choix des ?i doit conduire à un système algébrique à solution unique. Des théorèmes d'analyse numérique suggèrent certains choix, pour ~
lesquels la solution approchée f est garantie unique et convergente vers la solution exacte 2.9 quand le maillage se raffine. Dans notre exemple, si on choisit la formulation (2.9), les fonctions ?i(x) doivent avoir une dérivée non nulle partout, pour que la formulation ne soit pas triviale. Par contre, dans cette formulation, la dérivée seconde de la fonction f(x), qui sera
~
remplacée par l'approximation f (x), n'apparaît plus.
On pouvait donc choisir une famille d'approximation
~
f (x)plus simple.
On choisit par exemple les 7 fonctions test ?1, ?2 ,…, ?7 définies par la figure 5.
Figure 5:Les sept fonctions ?i(x).
Leurs dérivées sont constantes (+6 ou -6) ou nulles et on a ?i(xj) = ?ij
Pour chaque fonction test ?i. on écrit une équation scalaire:
~
où les f (x) sont des fonctions de x et des valeurs aux noeuds f1,f2,...,f7.
Les intégrales de ces 7 équations sont faciles à calculer : elles se réduisent à des intégrales sur le segment où les dérivées de ?i sont non nulles.
La construction de ce système d'équations s'appelle l'assemblage. Pour un choix différent de fonctions ?i(x), les intégrales sur chaque élément sont moins simples à calculer.
Par exemple, la première équation (avec ?1) est :
? ?
avec a, b et c définis en (2.1). On procède de même pour les 6 autres équations.
On obtient ainsi un système de 7 équations à 7 inconnues (les valeurs aux noeuds f1,f2,...,f7) :
La solution de ce système donne les valeurs aux
~
noeuds de l'approximation f (x) cherchée :
f 1 = f 2 =? f 3 =? f 4 =? f 5 =? f 6 =? f7 =?
Les choix faits précédemment (les fonctions tests et les fonctions d'interpolation) conduisent à un système d'équations non symétrique. On verra plus loin comment obtenir des systèmes symétriques, pour lesquels on dispose d'algorithmes numériques plus efficaces.
Les fi étant déterminées, on connaît alors les interpolations dans les trois sous domaines, et donc aussi la solution approchée en juxtaposant ces interpolations.
La figure 6 montre une comparaison entre la solution exacte et la solution approchée.
Figure 6:Comparaison des solutions exacte et approchée
On peut noter que puisqu'on n'a pas imposé aux
~ fonctions f (x)de respecter exactement les conditions aux limites, celles-ci ne le sont qu'approximativement:
Bien que cela n'apparaisse pas clairement sur la figure 6, les dérivées sont bien discontinues :
À titre de comparaison on donne figure 7 le graphe de la solution approchée avec le même maillage mais avec des éléments à interpolation linéaire à deux
~
noeuds. L'espace des F(x)n'est que de dimension 4 et la solution est beaucoup plus grossière.
Figure 7:Approximation par 3 éléments linéaires
Figure 2.8: Approximation par 3 éléments
linéaires,avec C.L. imposées exactes
Figure 9:Approximation par 3 éléments d'ordre 2, avec C.L. imposées exactes
On donne aussi figure 8 le graphe avec 3 éléments à interpolation linéaire à deux noeuds où on a astreint les interpolations des éléments 1 et 3 à respecter exactement les conditions aux limites. L'espace des
~
F(x)n'est plus que de dimension 2.
On donne enfin figure 9 le graphe avec 3 éléments à interpolation du second degré à 3 noeuds où on a astreint les interpolations des éléments 1 et 3 à respecter exactement les conditions aux limites. ~
L'espace des F(x) est de dimension 5.
Le pdetool sous MatLab est un outil informatique permettant le traitement des équations aux dérivées partielles en deux dimensions en utilisant la méthode des éléments finis. Ce document constitue un guide simple permettant l'utilisation efficace de cette interface graphique.
Tapez sous la fenêtre des commandes de MatLab la commande pdetool, vous devez recevoir l'interface graphique suivante.
Nous allons traiter l'équation de Laplace avec des conditions aux limites de type Dirichlet. Ce problème décrit généralement un équilibre des forces, un bilan d'énergie thermique ou électrostatique. L'exemple suivant peut être résolu à la main en faisant appel au développement par séries de Fourier. Ceci permet de vérifier l'exactitude des résultats numériques obtenus par pdetool. Cependant, vous pouvez compliquer par la suite votre problème à des cas non évidents pour un traitement analytique.
1. Dessin du domaine: supposons qu'il s'agit durectangle [0, 3]3[0, 1].
Dans le menu Options, sélectionner Grid and Snap; puis choisir Axes' Limits entre [-1 4] et [-1 2]. Faites Apply et fermer la boite Axes' Limits. Le résultat sera affiché à travers la figure suivante.
Dans le menu Draw, sélectionner Rectangle/Square (non centré).
Dans la zone du dessin, placer le curseur sur le point (0,0), cliquer en permanence avec le bouton gauche de la souris en déplaçant le curseur jusqu'au point (3,1). La région marquée sera désignée par R1.
2. Passer aux conditions aux limites: donner à lasolution u recherché la valeur -1 sur les cotes horizontales, et +1 sur les cotes verticales.
Dans le menu Boundary, select Boundary Mode. La figure change comme suit :
Dans la région du dessin, double cliquer sur la cote gauche vertical du rectangle. La boite de dialogue boundary condition apparaît. La condition du type Dirichlet h*u=r est prise par défaut. La valeur h=1 est correcte; changer la valeur de r vers +1.
Faite la même chose pour la cote verticale droite du triangle, c.a.d. donner à r la valeur -1. Refaire aussi le même travail sur les deux cotes horizontales en donnant à r la valeur +1. Fermer la boite de dialogue.
3. Spécification des paramètres de l'équation PDE :
Dans le menu PDE, sélectionner PDE specification. La boite de dialogue pde specification apparaît; le problème elliptique de Laplace -div(c*grad(u)) + a*u
= f, avec a=0, f=0, et c=1 est donné par défaut.
Dans le menu Mesh, sélectionner Initialize Mesh. Vous obtenez un maillage régulier et standard.
Dans le menu Solve, sélectionner Solve PDE. La figure dans la section du dessin change vers la solution du problème.
6. Génération des autres présentations pour la solution:
Dans le menu Plot, sélectionner Parameters. La boite de dialogue plot selection apparaît; cliquer sur les boutons de votre choix pour le type de présentation, les paramètres des abscisses et des ordonnées, la couleur et terminer en appuyant sur le bouton plot.
Donner pour la cote horizontale inférieure la condition
Si vous adoptez les paramètres de la figure ci-dessus, -.05 * x.^4 .* (3 - x).^2 et pour la cote supérieur la nouvelle figure suivante apparaît. horizontale la condition -.05*x.^2.*(3 - x).^4.
Pour sauvegarder votre travail, cliquer dans le menu File sur save as… Donner à votre programme le nom "monpde1.m" et quitter MatLab. Plus tard, vous pouvez reprendre l'exécution de toutes les phases de votre programme en tapant sous la fenêtre des commandes la simple commande monpde1.
Il existe dans la barre du menu une série des boutons pouvant remplacer un nombre assez important des commandes contenues dans les sous rubriques de chaque menu. Il est important de les remarquer et les retenir.
4.1. Ecoulement dans un convergent :
C'est un exemple d'écoulement potentiel, une application classique dans la mécanique des fluides incompressibles. Il s'agit d'un fluide newtonien qui traverse une tube conique.
Si le flux de fluide est irrotationnel, alors le vecteur v de vitesse est exprimé comme gradient d'un potentiel inconnu u, c.-à-d. v=grad(u). Avec la densité constante, la conservation de la masse donne div(v)=0, et nous devons ainsi résoudre l'équation div(grad(u))=0.
Le fluide entre de la gauche avec une vitesse constante v=v0. Il sort à la droite avec une vitesse inconnue, mais avec un potentiel constant u=0. le débit est nul par les autres segments de frontière.
Le domaine numérique:
Le tube convergent dans un repère cylindrique.
En raison de la symétrie, nous pouvons ramener les coordonnées de cylindrique (x, r, thêta) du problème à seulement deux coordonnées, (x,r). À la frontière de création récente, nous présentons l'écoulement parallèle à l'axe des abscisses, ce qui induit une autre condition aux limites Neumann.
Sur les frontières pleines, la vitesse dans la direction normale disparaît, qui implique des états de frontière de Neumann, c.-à-d., vous indiquez le composant normal du grad(u) à l'apport, la vitesse est perpendiculaire à la frontière, aussi un état de Neumann. La vitesse sur toutes les frontières ne peut pas déterminer le potentiel, parce qu'ajouter des constantes à u ne change pas v. L'écoulement est régi par la différence potentielle plutôt puis le potentiel lui-même. À la sortie x=1, nous ne savons pas la vitesse, mais nous pouvons assumer un écoulement laminaire loin du goulot. Là nous pouvons fixer le potentiel à la constante, par exemple, u=0, qui est un état de frontière de Dirichlet.
Employer des symétries et les coordonnées de cylindrique, les changements d'équation à : y*grad(u))=0 de div(. Note : À l'apport devient vous devez indiquer le composant normal du y*grad(u), qui est 10y.
L'équation est linéaire même si les coefficients sont variables. Appuyer sur le bouton de solution produit la solution, qui est le potentiel u. Le champ v de vitesse est le gradient du potentiel et est visualisé en utilisant la parcelle de terrain de flèche pour le grad(u).
Les lignes de niveau donne la taille de v, alors que les flèches donnent la distribution de vitesse.
Le flux est inverse proportionnel au secteur de la section transversale. Puisque le secteur de la sortie est 4 fois plus petit puis que de l'apport, la vitesse augmente par un facteur 4, clairement vu en comparant la longueur des flèches au dans et à la sortie. Notez la singularité au coin de réentrée, où la vitesse devient infinie. Ceci indique la validité limitée du modèle potentiel d'écoulement pour des problèmes réalistes d'écoulement.
A soap film bounded by a wire frame has the shape determined by the minimization of its internal energy. The shape of the film is described by the heightu(x,y) over a given domain. The height is fixed at the domain boundary, by the height of the wire frame g If the film is homogeneous, minimizing the internal energy is equivalent to minimizing its area, the integral of s(u) = sqrt(1+|grad(u)|^2).
This in its turn is equivalent to finding u such that
• u(x,y) = g(x,y), on the domain boundary In this example you look at a wire frame with two circular components. The inner frame is at height 0 while the outer frame is given by g(x,y) = x^2.
The wire frame for a soap bubble. The computational domain Geometry
First you draw the domain which is an annulus (a twodimensional doughnut), i.e. the region between two concentric circles:
An annulus obtained by subtracting the small disk from the large one. This is expressed by the set formula C1-C2
Next you set the boundary conditions, u = 0 on the inner frame, and u(x,y) = x^2 on the outer frame:
Boundary conditions dialog box.
Setting the equation coefficients amounts to typing in MATLAB statements. These can be constants, expressions containing the space variables x and , y , the unknown u , its partial derivatives ux and uy , or even user-defined functions. The general form of the time-independent, elliptic, scalar PDE is
In the soap bubble problem, you have a = f = 0 and c = 1./sqrt(1+ux.^2+uy.^2).
A triangularization of the domain
The mesh generator initializes a coarse mesh. Further refinements of the mesh are necessary to obtain a reasonable solution. Mesh generation is done automatically, but it you can also steer it by setting parameters. In this example the default mesh is initialized and then uniformly refined. The fine mesh contains 590 nodes and 1096 triangles.
The Finite Element Method discretization is performed. Linear problems are solved by MATLAB's internal solver. However, the soap bubble equation is nonlinear, and you have to tell that to the solver. After setting the nonlinear option among the Solve parameters, the equations are discretized and solved. A damped Gauss-Newton iterative is employed. Set the Jacobian option to Full. While the solver works, you can see the reduction of the residual in your workspace window.
The shape of the soap film
This is an example of a non-linear system of two state-variables, electric potential u and temperature v, in equilibrium.
An electric current passing through a metal will produce heat. It is known that the electric conductivity depends on the temperature of the conductor. Thus the heat and current distributions are related to each other.
The electric potential u satisfies the equation -div( s*grad(u) )=0, where s=s(v) is a function of the temperature v. The current density -s*grad(u) produces heat at the rate s|u|^2. This is a source term in the heat equation. Over a wide range of temperatures, the electric conductivity is inverse proportional to the temperature. The system of PDEs modeling the electric and heat conductance is:
• -div( s(v)*grad(u) ) = 0
• -K*div( grad(v) ) = s(v)/|grad(u)|^2 Imagine a rectangular metal plate with a potential difference of unit size applied on two opposite faces. Those faces are at ``room'' temperature, 300 degrees Kelvin. The other two faces are insulated.
A metal plate with an electric current heating it.
The geometry is very simple, a unit square. Select the Generic system mode where you define the boundary conditions end the equations.
There are two kind of boundary conditions: On the insulated sides, we select zero (default) Neumann boundary conditions.
On the other two sides, we select Dirichlet boundary conditions, setting the potentials r1 to 0 and 1, and the temperature r2 to 300. The Equation
The generic system is:
-div( c11*grad(u) ) - div( c12*grad(v) ) + a11*u +a12*v = f1
In this case, c11 = 1e+8./u(2,:), c22 = 400, and f2 = 1e+8./u(2,:).*(ux(1,:).^2+uy(2,:).^2).
Specification of a nonlinear system. u, ux, uy contain two rows, one for each unknown.
One uniform refinement of the default original mesh produce enough mesh-points (689 points in the unit square). Next you have to set the options for the solver. You have to tell the solver that the problem is nonlinear and you can ``help'' by giving a good initial guess of the solution. There are other parameters that are used for fine-tuning the solver in specific cases, and these should be used when the default values fail. The default initial guess u=v=0 will fail because the current density becomes infinite. You know the temperature on two sides, so set u(2,:)=300. We also know that the potential varies, e.g., from 0 to 1 in the x-direction, so set u(1,:)=x.
Options in the nonlinear solver. Click to enlarge.
Pressing the solve-button starts the nonlinear iterations (a damped Gauss-Newton iteration). The nonlinear solver performs 17 steps to reduce the residual to the specified tolerance. Using the full approximation of the Jacobian, the solution is found (up to the same tolerance) in only 4 steps. Since the computing time for the full Jacobian is long, this strategy should be used only for extremely tight tolerances.
Heating due to electric currents.
Using the plot preferences, you can put different properties on the same image. The color scale is a thermometer showing the heat distribution in the metal plate. Note that the change in electric conductance is not means that the equation for the electric potential $u$ cannot be uncoupled from the heat equation.
Note that the middle section of the plate becomes very hot. You can cool it by circulating ``room temperature'' air through a hole. Carve out a circular hole with diameter 0.2 in the center of the plate. At the new boundary, the temperature u(2,:) is fixed again to 300 degrees Kelvin. Since this boundary is electrically insulated, the current passing through it is zero, i.e. a Neumann condition is needed. This is an example of mixed Neumann-Dirichlet conditions. The Dirichlet condition is expressed by the constraint h*[0 1;1 0]*[u;v]=[300;0]
Mixed Neumann-Dirichlet boundary conditions. The electric potential satisfies a Neumann condition, while the temperature is constrained by the Dirichlet condition v=300.
The cooling effect is substantial. The maximum temperature drops from 117 to 83 Celsius:
Cooling the heated plate.
The air cooling through the middle of the plate produces a significant drop in the temperature. Note the strong variations in the current density around the cooling hole, which makes the problem genuinely nonlinear.
4.4 . Stresses in a Wrench
In draw mode you create a Constructive Solid Geometry model of the geometry. Different building blocks can be combined using a set formula. The picture on the right shows the Constructive Solid Geometry model of a wrench.
Define the boundary conditions by double clicking on the boundary segments. You can select Dirichlet, Neumann or mixed boundary conditions. Using the Dirichlet boundary condition you clamp the "claw" part of the wrench and with the Neumann boundary condition you apply a downward force to the handle.
In mesh mode the geometry is automatically meshed. The picture to the right shows the result of the adaptive mesh refinement algorithm. The triangular mesh has been refined the most where the variation in stress is large.
In plot mode you have an almost unlimited number of plot options. The picture on the left shows the x direction stresses visualized using contour and color plots. The plot clearly shows tension (red) and compression (dark blue) in different parts of the wrench. The displacements are shown by the distorted figure.
