* fichier : ther9.dgibi ************************************************************************ * Section : Thermique Conduction ************************************************************************ ************************************************************************* * NOM : THER9 * DESCRIPTION : Cas-test du p-laplacien (p>=1) inspiré par \cite[chapitre * 8]{saramito}. * * Teste la thermique en massif 2D conduction anisotrope. * * * Il s'agit de l'équation non-linéaire suivante : * * - div k(T) grad T = f * * avec terme source : f = 1 * et coefficient de conductivité : k(T) = |grad T|^(p-2) * condition aux limites T = 0 sur le bord du domaine. * * Cette équation est intéressante à plusieurs titres : * - la régularité de la solution varie avec p ; * - pour p=2, on a le Laplacien normal * - pour p=1, on a l'opérateur courbure moyenne * - pour p-> +\inf, la solution T tend vers la fonction * distance au bord \cite{belyaev}. * - suivant les valeurs de p, le comportement des * algorithmes de résolution non-linéaire (point fixe, * Newton, accélération de convergence) varie. * * Pour la méthode de Newton, la matrice tangente du * p-Laplacien s'exprime sous la forme d'un Laplacien avec * un tenseur de diffusion anisotrope : * * k2(T) = k(T) I + (p-2) |grad T|^(p-4) gradT \ptensoriel gradT * * * * Sur la d-sphère unité, pour p>1, on a une solution * radiale analytique de l'équation : * T = C (1-r^a) * * avec : a = p / (p-1) * C = d^(-1/(p-1)) / a * * En effet, pour une solution T radiale ne dépendant que * de r, le p-Laplacien s'exprime comme \cite{franzina} : * -\delta_p T(r) = -(p-1) |T'(r)|^(p-2) T''(r) * -(d-1)/r |T'(r)|^(p-2) T'(r) * * * Bibliographie (bibtex) : *@unpublished{saramito, * TITLE = {{Efficient C++ finite element computing with Rheolef}}, * AUTHOR = {Saramito, Pierre}, * URL = {https://cel.archives-ouvertes.fr/cel-00573970}, * NOTE = {Notes de cours de DEA accessibles depuis : \url{https://cel.archives-ouvertes.fr/cel-00573970} }, * TYPE = {DEA}, * ADDRESS = {Grenoble, France, France}, * PAGES = {176}, * YEAR = {2016}, * MONTH = Jun, * KEYWORDS = {Finite Element Method FEM ; C++ LIBRARY ; Navier Stokes Equations ; Characteristic Method ; Elasticity equations}, * PDF = {https://cel.archives-ouvertes.fr/cel-00573970/file/rheolef.pdf}, * HAL_ID = {cel-00573970}, * HAL_VERSION = {v13}, *} * *@article{franzina2010existence, * title={Existence and uniqueness for a p-laplacian nonlinear eigenvalue problem}, * author={Franzina, Giovanni and Lamberti, Pier Domenico}, * journal={Electron. J. Differential Equations}, * volume={26}, * number={10}, * year={2010} *} * * *@article{belyaev, *author = {Belyaev, Alexander G. and Fayolle, Pierre-Alain}, *title = {On Variational and PDE-Based Distance Function Approximations}, *journal = {Computer Graphics Forum}, *volume = {34}, *number = {8}, *issn = {1467-8659}, *url = {http://dx.doi.org/10.1111/cgf.12611}, *doi = {10.1111/cgf.12611}, *pages = {104--118}, *year = {2015}, *} * * * * LANGAGE : GIBIANE-CAST3M * AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA) * mél : stephane.gounand@cea.fr ********************************************************************** * VERSION : v1, 02/05/2017, version initiale * HISTORIQUE : v1, 02/05/2017, création * HISTORIQUE : * HISTORIQUE : ************************************************************************ 'OPTI' 'ECHO' 1 ; 'SAUT' 'PAGE' ; interact = faux ; graph = faux ; * * Procédure d'accélération de la convergence. * 'DEBP' ACCCVG ; 'ARGU' res*'CHPOINT' ; 'ARGU' tabacc*'TABLE' ; 'SI' ('NON' ('EXIS' tabacc 'iter')) ; tabacc . 'iter' = 0 ; tabacc . 'znacce' = 4 ; tabacc . 'itdep' = 3 ; tabacc . 'acfp1' = 'COPI' ('*' res 0.) ; tabacc . 'acfp2' = tabacc . 'acfp1' ; tabacc . 'acfp3' = tabacc . 'acfp1' ; tabacc . 'acfep1' = tabacc . 'acfp1' ; tabacc . 'acfep2' = tabacc . 'acfp1' ; tabacc . 'correc' = 0. ; tabacc . 'freap' = 0. ; 'FINS' ; it = '+' (tabacc . 'iter') 1 ; znacce = tabacc . 'znacce' ; itdep = tabacc . 'itdep' ; tabacc . 'iter' = it ; tabacc . 'correcp' = tabacc . 'correc' ; tabacc . 'correc' = 0. ; tabacc . 'acfp0' = 'ENLEVER' ('-' res (tabacc . 'freap')) 'FLX' ; tabacc . 'acfep0' = '-' (tabacc . 'acfp0') (tabacc . 'correcp') ; 'SI' ('ET' ('MULT' it znacce) ('>' it itdep)) ; 'MESS' 'acceleration' ; correc = 'ACT3' (tabacc . 'acfep2') (tabacc . 'acfep1') (tabacc . 'acfep0') (tabacc . 'acfp3') (tabacc . 'acfp2') (tabacc . 'acfp1') (tabacc . 'acfp0') ; res = '-' res correc ; 'FINS' ; 'SI' ('>' it 3) ; 'DETR' (tabacc . 'acfp3') ; 'DETR' (tabacc . 'acfep2') ; 'FINSI' ; tabacc . 'acfp3' = tabacc . 'acfp2 ' ; tabacc . 'acfp2' = tabacc . 'acfp1 ' ; tabacc . 'acfp1' = tabacc . 'acfp0 ' ; tabacc . 'acfep2' = tabacc . 'acfep1 ' ; tabacc . 'acfep1' = tabacc . 'acfep0 ' ; 'RESP' res tabacc ; 'FINP' ; * * Options globales * 'OPTI' 'DIME' 2 ; 'OPTI' 'ELEM' 'TRI3' ; * * Maillage * * nx : nombre de maille par côté * idom : type de domaine 0 (carré) ou 1 (quart de cercle) * Sur 1, il y a une solution analytique pour p>1 nx = 40 ; idom = 1 ; * Physique * vp : exposant p du p-laplacien vp = 3. ; * Numérique * nitmax : nombre d'itérations max. * tol : tolérance pour la convergence de la boucle non-linéaire * omega : facteur de relaxation * iktan : type de matrice tangente : * 0 laplacien ; 1 point fixe ou ; 2 newton nitmax = 100 ; tol = 1.D-10 ; *omega = '/' 2. vp ; *iktan = 1 ; omega = 1. ; iktan = 2 ; lacc = 0 ; tabacc = 'TABL' ; chpar = 'CHAI' 'p=' vp ' omega=' omega ' iktan=' iktan ' lacc=' lacc ; 'MESS' '*************************' ; 'MESS' '* P laplacien ' chpar ; 'MESS' '*************************' ; * Points 'OPTI' 'DENS' ('/' 1. ('FLOT' nx)) ; pA = 0. 0. ; pB = 1. 0. ; pC = 1. 1. ; pD = 0. 1. ; * * Maillage carré * 'SI' ('EGA' idom 0) ; * Lignes lAB = 'DROI' nx pA pB ; lBC = 'DROI' nx pB pC ; lCD = 'DROI' nx pC pD ; lDA = 'DROI' nx pD pA ; * Surface *mt = 'DALL' lAB lBC lCD lDA ; mt = 'SURF' (lAB 'ET' lBC 'ET' lCD 'ET' lDA) ; * Blocages mailblo = 'CONT' mt ; 'FINS' ; * * Maillage quart de cercle * 'SI' ('EGA' idom 1) ; * Lignes lAB = 'DROI' pA pB ; lBD = 'CERC' pB pA pD ; lDA = 'DROI' pD pA ; * Surface *mt = 'DALL' lAB lBC lCD lDA ; mt = 'SURF' (lAB 'ET' lBD 'ET' lDA) ; mailblo = lBD ; 'FINS' ; * 'SI' graph ; 'TRAC' mt 'TITR' 'MAILLAGE' ; 'FINS' ; * * Modele et materiau initial * modt = 'MODE' mt 'THERMIQUE' 'ISOTROPE' ; modt2 = 'MODE' mt 'THERMIQUE' 'ANISOTROPE' 'CONS' 'T2'; matt = 'MATE' modt 'K' 1. ; * mlap0 = 'COND' modt matt ; * terme source tsou = 'SOUR' modt 1. mt ; * conditions aux limites mcli = 'BLOQ' 'T' mailblo ; * résolution mat = mlap0 'ET' mcli ; smb = tsou ; sol = 'RESO' mat smb ; solt = 'EXCO' 'T' sol 'T' ; *ch0 = 'Solution au pas 0' ; *'TRAC' solt mt 'TITR' ch0 ; ipas = 0 ; soli = sol ; lnres = 'PROG' ; cgtx = 'MOTS' 'T,X' ; cgty = 'MOTS' 'T,Y' ; cgt = cgtx 'ET' cgty ; ck11 = 'MOTS' 'K11' ; ck21 = 'MOTS' 'K21' ; ck22 = 'MOTS' 'K22' ; * * Boucle pour la non-linearite * 'REPE' bcl nitmax ; gsoli = 'GRAD' modt soli ; ngsoli = 'PSCA' gsoli gsoli cgt cgt ; * Nouveau matériau avec la nouvelle conductivité chk = ngsoli '**' ('/' ('-' vp 2.) 2.) ; matt = 'MATE' modt 'K' chk ; * Matrice et calcul du résidu mlap = 'COND' modt matt ; mat = mlap 'ET' mcli ; res = tsou '-' ('*' mat soli) ; * nres = 'MAXI' res 'ABS' 'AVEC' ('MOTS' 'Q') ; lnres = 'ET' lnres ('PROG' nres) ; 'SI' ('<' nres tol) ; 'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres ' < tol=' tol) ; 'QUIT' bcl ; 'SINO' ; 'MESS' ('CHAI' 'Pas i=' ipas ' max|Ri|=' nres) ; 'FINS' ; * Matrice tangente = rigidite au pas 0 'SI' ('EGA' iktan 0) ; mattan = mlap0 ; 'FINS' ; * Matrice tangente = Matrice de point fixe 'SI' ('EGA' iktan 1) ; mattan = mlap ; 'FINS' ; 'SI' ('EGA' iktan 2) ; * Matrice tangente = Matrice de Newton * 2 \eta' (|grad T|^2) (grad T \ptens grad T) * avec \eta'(z) = \frac{p-2}{2} z^{\frac{p-4}{2}} coeff = '*' (ngsoli '**' ('/' ('-' vp 4.) 2.)) ('-' vp 2.) ; * On change le type de gsoli sinon il veut une syntaxe avec un modele gsoli = 'CHAN' 'TYPE' gsoli ' ' ; vk11 = '*' gsoli gsoli cgtx cgtx ck11 ; vk21 = '*' gsoli gsoli cgty cgtx ck21 ; vk22 = '*' gsoli gsoli cgty cgty ck22 ; matt2 = 'MATE' modt2 'DIRECTION' (1. 0.) 'PARALLELE' 'K11' ('CHAN' ('*' vk11 coeff) 'CONS' 'T2') 'K21' ('CHAN' ('*' vk21 coeff) 'CONS' 'T2') 'K22' ('CHAN' ('*' vk22 coeff) 'CONS' 'T2') ; matnew2 = 'COND' modt2 matt2 ; mattan = mlap 'ET' matnew2 ; 'FINS' ; mattot = ('*' mattan ('/' 1. omega)) 'ET' mcli ; smb = res ; 'SI' ('EGA' lacc 1) ; smb tabacc = ACCCVG smb tabacc ; 'FINS' ; dsol = 'RESO' mattot smb ; 'SI' ('EGA' lacc 1) ; tabacc . 'freap' = '*' mcli dsol ; 'FINS' ; solip = '+' soli dsol ; * solipt = 'EXCO' 'T' solip 'T' ; * chip = 'CHAI' 'Solution vp=' vp ' au pas ' ('+' ipas 1) ; * 'TRAC' solipt mt 'TITR' chip 'NCLK' ; soli = solip ; ipas = '+' ipas 1 ; 'FIN' bcl ; * * Courbe du résidu * npas = 'DIME' lnres ; 'SI' ('>' npas 1) ; lpas = 'PROG' 0. 'PAS' 1. 'NPAS' ('-' npas 1) ; llnres = '/' ('LOG' lnres) ('LOG' 10.) ; evl = 'EVOL' 'MANU' lpas llnres ; 'SI' graph ; titx = 'iter' ; tity = 'Log. |res|' ; tit = 'CHAI' tity '(' titx ')' ' ' chpar ; 'DESS' evl 'TITR' tit 'TITX' titx 'TITY' tity ; 'FINS' ; 'FINS' ; * * Tracé en élévation * 'SI' graph ; 'OPTI' 'DIME' 3 ; solit = 'EXCO' 'T' soli ; depz = 'NOMC' solit 'UZ' ; 'FORM' ('/' depz ('MAXI' depz 'ABS')) ; 'OPTI' 'OEIL' (20. 30. 10.) ; 'OPTI' 'ISOV' 'SULI' ; chi = 'CHAI' 'Solution vp=' vp ' au pas ' ipas ; 'TRAC' solit mt ('CONT' mt) 30 'TITR' chi ; 'OPTI' 'DIME' 2 ; 'FINS' ; * * Comparaison à la solution analytique * 'SI' ('EGA' idom 1) ; x y = 'COOR' mt ; r2 = (x '*' x) '+' (y '*' y) ; expo = vp '/' ('-' vp 1.) ; vd = 'VALE' 'DIME' ; Const = '/' (vd '**' ('/' -1 ('-' vp 1.))) expo ; solana = '*' ('-' 1. ('**' r2 ('/' expo 2.))) const ; solnum = 'EXCO' 'T' soli ; dsol = '-' solana solnum ; 'SI' graph ; 'TRAC' dsol mt ('CONT' mt) 30 'TITR' 'Difference analytique-numerique' ; 'FINS' ; mdsol = 'MAXI' dsol 'ABS' ; ch = 'CHAI' 'Maxi. |analytique - numerique|=' mdsol ' ' chpar; 'MESS' ch ; 'FINS' ; * * Tests : * 1) La méthode de Newton converge en 7-8 itérations pour p=3 (cf. Saramito * p.137) * 2) L'erreur par rapport à la solution analytique vaut environ 4.e-4 * avec 40 mailles * lok = vrai ; lok = 'ET' lok (npas '