Thermoelastic axisymmetric problem

Dear All,

I still try to solve a circular mirror deformation (in steady state) heated by a laser beam in its center.

First, I solve the temperature axisymmetric problem (with radiation) similar to the freefem documentation example => OK

Then, I use same equations than the “system of elasticity” freefem example for which I added the temperature field :

image

and then, this equation is written in the axisymmetric case.

I got an analytical code which is able to plot the surface deformation of the mirror and then, I can compare with the result from freefem.
With the analytical code, the surface deformation has almost the same shape than the temperature of the surface.
unfortunately, freefem does not give the same result.

I thought it was coming from the axisymmetric problem writting for which the div(u) and strain tensor let appear some 1/r terms which are maybe not defined for r=0.
thus, I wrote the problem with new terms U=u/r to make 1/r terms disapear in the formulation… but the result is still the same.

I checked the stress tensor solution on the borders and it seems to be properly 0.
but I also compute div(stress tensor) which should be almost zero… but it is not !

thus, I think I’m doing something wrong… but I don’t see what…
thank you if you can help

best regards
ronic

here is the code :

// define macro to reduce equation writting
macro epstens(Ur,uz) [Ur+x*dx(Ur),Ur,dy(uz),(x*dy(Ur)+dx(uz))/sqrt(2.)]
//
macro grad(u) [dx(u),dy(u)]
//
macro div(Ur,uz) ( 2*Ur+x*dx(Ur)+dy(uz) )

// source function
real wz=1e-3;
sf=SE*exp(-2*x^2/wz^2);

// setting equations problem with dr=>dx and dz=>dy
problem TempSimul(dT,psi) = int2d(Th)(k*grad(psi)' *grad(dT))
-int1d(Th,B3d,B3r)(psi*sf)+int1d(Th,B1r,B1f,B2,B3d,B3r)(psi*K1*dT);

// ur=r*Ur, phir=r*PHIr						
problem DeformSimul([Ur,uz],[PHIr,phiz]) = 
int2d(Th)(2*mu*epstens(PHIr,phiz)' *epstens(Ur,uz))
+int2d(Th)(div(PHIr,phiz)*lambda*div(Ur,uz))-int2d(Th)(div(PHIr,phiz)*nu*dT);

// solving problem
dT=0;
for (int n=0;n<5;n++)
{
// Linearization : T^4-Text^4 = (T+Text)*(T^2+Text^2)*(T-Text) = K1*(T-Text)
K1=epsilon*sigma*((Text+dT)+Text)*((Text+dT)^2+Text^2); 
TempSimul;
}
DeformSimul;
ur=x*Ur;
phir=x*PHIr;