# Thermoelasticity - Free FEM

**URL:** <https://community.freefem.org/t/thermoelasticity-free-fem/1456>\
**Category:** General Discussion\
**Created:** [January 28, 2022, 7:58am UTC](https://community.freefem.org/t/thermoelasticity-free-fem/1456 "2022-01-28T07:58:45Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![amir](https://avatars.discourse-cdn.com/v4/letter/a/c37758/32.png) [@amir](https://community.freefem.org/u/amir)\
**Post date:** [January 28, 2022, 7:58am UTC](https://community.freefem.org/t/thermoelasticity-free-fem/1456/1 "2022-01-28T07:58:45Z")

</div>

Hi,  
How can I write the code of the following equation??

![11111111](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/48871095bb85bf2bf39e37b98ffa6f71d5477c9a.jpeg)

epsilon is strain, and T hat is weight function.

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [January 28, 2022, 3:26pm UTC](https://community.freefem.org/t/thermoelasticity-free-fem/1456/2 "2022-01-28T15:26:24Z")

</div>

A small exemple of Thermo elasticity with

- alpha \* T\*([1, 1, 0.]’\* AA\*epsilon(uu,vv)) = \alpha T\ tr( AA \*\varepsilon(uu,uu)) 

This is not exactly the same but it is simple modif.

```auto
// file lame.edp
mesh Th=square(10,10,[20*x,2*y-1]);
fespace Vh(Th,P2);
Vh u,v,uu,vv,T,TT;

macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))*0.5]// EOM
macro AA [ 
 [l2mu, lambda, 0.],
 [lambda, l2mu, 0.],
 [0., 0., mu2] ] //
macro Grad(u) [dx(u),dy(u)] //

real E = 21e5, nu = 0.28, mu= E/(2*(1+nu)); 
real lambda = E*nu/((1+nu)*(1-2*nu)), f = 0; //  
real l2mu = lambda + 2*mu;
real mu2 = mu*2;
real K=1;
real ft=1;
real alpha =1e-2;

solve lame([u,v,T],[uu,vv,TT])
  = int2d(Th)(  
               epsilon(u,v)'*AA*epsilon(uu,vv) 
            - alpha * T*([1, 1, 0.]'* AA*epsilon(uu,vv))
            + K*(Grad(T)'*Grad(TT))
             )
        - int2d(Th)(f*vv+ ft*TT) 
        + on(4,u=0,v=0,T=0);
real coef=100;
plot([u,v],T,wait=1,ps="lamevect.eps",coef=coef);
plot(T,wait=1,fill=1);

mesh th1 = movemesh(Th, [x+u*coef, y+v*coef]);
 plot(th1,wait=1,ps="lamedeform.eps");
real dxmin = u[].min; 
real dymin = v[].min; 

cout << " - dep. max x = "<< dxmin<< " y=" << dymin << endl; 
cout << " dep. (20,0) = " << u(20,0) << " " << v(20,0) << endl; 

```