1// non linear elasticity model 2// 3// ------------------------------- 4// with huge utilisation of macro new version 5// more simple 6// --------------------------- 7// optimize version 8// ------------ 9// problem is find $(uu,vn)$ minimizing $J$ 10// $ min J(un,vn) = \int 1/2 (F2) - int Pa * un $ 11// $ dJ(u,u,uu,vv) = int dF2(u,v,uu,vv) df(F2(u,v)) $ 12// where $F2 = (^t {E} A {E} ) $, 13// $E(U) = 1/2 (\nabla U + \nabla U^t + \nabla U^t \nabla U) $ 14// ($u_1$) 15// with U=( ) 16// ($u_2$) 17// so: 18//$$ E_ij = 0.5 ( d_i u_j + d_j u_i ) + \sum_k d_i u_k * d_j*u_k \leqno(1)$$ 19// the symetric tensor $t_{ij}$ are a vector [t11,2*t12,t22] 20// date from 21/* 22Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow 23 24Stefan Turek and Jaroslav Hron⋆ 25Institute for Applied Mathematics and Numerics, University of Dortmund, 26Vogelpothsweg 87, 44227 Dortmund, Germany 27 28http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.550.1689&rep=rep1&type=pdf 29 30section 4.2 CSM test: 31*/ 32 33macro EL(u,v) [dx(u),(dx(v)+dy(u)),dy(v)] // is $[\epsilon_{11},2\epsilon_{12},\epsilon_{22}]$ 34 35macro ENL(u,v) [ 36(dx(u)*dx(u)+dx(v)*dx(v))*0.5, 37(dx(u)*dy(u)+dx(v)*dy(v)) , 38(dy(u)*dy(u)+dy(v)*dy(v))*0.5 ] // EOM ENL 39 40macro dENL(u,v,uu,vv) [(dx(u)*dx(uu)+dx(v)*dx(vv)), 41 (dx(u)*dy(uu)+dx(v)*dy(vv)+dx(uu)*dy(u)+dx(vv)*dy(v)), 42 (dy(u)*dy(uu)+dy(v)*dy(vv)) ] // 43 44 45macro E(u,v) (EL(u,v)+ENL(u,v)) // is $[\E{11},\sqrt2\E{12},\E{22}] 46macro dE(u,v,uu,vv) (EL(uu,vv)+dENL(u,v,uu,vv)) // 47macro ddE(u,v,uu,vv,uuu,vvv) dENL(uuu,vvv,uu,vv) // 48macro F2(u,v) (E(u,v)'*A*E(u,v)) // 49macro dF2(u,v,uu,vv) (E(u,v)'*A*dE(u,v,uu,vv) ) // 50macro ddF2(u,v,uu,vv,uuu,vvv) ( 51 (dE(u,v,uu,vv)'*A*dE(u,v,uuu,vvv)) 52 + (E(u,v)'*A*ddE(u,v,uu,vv,uuu,vvv)) )// EOM 53 54// for hyper elasticity problem 55// ----------------------------- 56macro f(u) (u) // end of macro 57macro df(u) (1) // end of macro 58macro ddf(u) (0) // end of macro 59 60// ------------------------------- 61// 62// $ \sigma = 2 \mu E + \lambda tr(E) Id $ 63// $ A(u,v)= \sigma(u):\E(v) $ 64// 65// ( a b ) 66// ( b c ) 67// 68// tr*Id : (a,b,c) -> (a+c,0,a+c) 69// so the associed matrix is: 70// ( 1 0 1 ) 71// ( 0 0 0 ) 72// ( 1 0 1 ) 73// ------------------ 74verbosity=0; 75 76// data CMS 1 77real xp=0.35,yp=0.0;// 78real upbench= -0.007187, vpbench= -0.066; 79 80real EE = 1.4e6; 81real rho = 1000; 82real sigma = 0.4; 83real mu = EE/(2*(1+sigma)); 84real lambda = EE*sigma/((1+sigma)*(1-2*sigma)); 85real gravity = -2; 86int nnn = 10; 87 88mesh Th = square(17*nnn,nnn,[0.35*x,0.02*(y-0.5)]); 89cout<<" Th nt : " << Th.nt<<endl; 90int bottom=1, right=2,upper=3,left=4; 91 92plot(Th); 93 94real a11= 2*mu + lambda ; 95real a22= mu ; // because [0,2*t12,0]' A [0,2*s12,0] = 2*mu*(t12*s12+t21*s21) = 4*mu*t12*s12 96real a33= 2*mu + lambda ; 97real a12= 0 ; 98real a13= lambda ; 99real a23= 0 ; 100// symetric part 101real a21= a12 ; 102real a31= a13 ; 103real a32= a23 ; 104func A = [ [ a11,a12,a13],[ a21,a22,a23],[ a31,a32,a33] ]; 105 106 107 108 109 110 111fespace Vh(Th,[P2,P2]); 112 113 114 115 116 117 118Vh [uu,vv], [w,s],[un,vn]; 119[un,vn]=[0,0];// intialisation 120[uu,vv]=[0,0]; 121 122 123 124real errb=0; 125// Newton's method 126// --------------- 127for (int i=0;i<10;i++) 128{ 129 cout << "Loop " << i << endl; 130 solve NonLin([uu,vv],[w,s])= 131 int2d(Th)( ddF2(un,vn,uu,vv,w,s) ) 132 - int2d(Th)( dF2(un,vn,w,s) - rho*gravity*s) 133 + on(left,uu=0,vv=0); 134 ; 135 real res = uu[].linfty ; // norme L^2 of [uu,vv] 136 errb = max(abs(un(xp,yp)-upbench) , abs(vn(xp,yp)-vpbench)); 137 138 un[] -= uu[]; 139 cout << " Linfty residual = " << res << endl; 140 cout << " u of A =" << un(xp,yp) << " " << vn(xp,yp) << " bench " << upbench << " "<< vpbench << " err bench =" << errb << endl; 141 plot(movemesh(Th, [x+un, y+vn]),wait=1,cmm=" i =" +i + " err =" + res); 142 143 if (res<1e-10) break; 144} 145 146cout << " err bench = " << errb << endl; 147//plot([un,vn],wait=1); 148mesh th1 = movemesh(Th, [x+un, y+vn]); 149plot(th1,wait=1,ps="nl-elas.eps"); 150assert( errb < 1e-2);