Solving elasticity Equations

[joelvuillon]

Hi,
I am working with Florian De Vuyst who is in touch with Daniel Wheeler  for 
solving Elasticity system with Fipy.

Daniel provides the following formulation 

> mesh = ...
>
>    u = CellVariable(mesh=mesh)
>    v = CellVariable(mesh=mesh)
>
> The equations will be,
>
>    Xhat = (1, 0)
>    Yhat = (0, 1)
>
>    u_eqn = DiffusionTerm(mu) \
>                  + mu * (u.getFaceGrad().dot(Xhat) * Xhat \
>                           + v.getFaceGrad().dot(Xhat) * Yhat).getDivergence() 
> \
>                  + ConvectionTerm(lambda * Xhat) \
>                  + lambda * v.getGrad().dot(Yhat) \
>                  + f
>
>    v_eqn = DiffusionTerm(mu) \
>                  + mu * (u.getFaceGrad().dot(Yhat) * Xhat \
>                           + v.getFaceGrad().dot(Yhat) * Yhat).getDivergence() 
> \
>                  + ConvectionTerm(lambda * Yhat) \
>                  + lambda * u.getGrad().dot(Xhat) \
>                  + g


for solving the system for U=(u,v)  :

div.[ mu grad(U) + mu [grad(U)]^T + lambda I tr(grad(U)) ] + F = (0,0).

However, I do not agree with that formulation for the particular term div.[ 
lambda I tr(grad(U)) ] which seems has been written as div.[ lambda (U) ].

Instead, I propose the following formulation :

       u_eqn = DiffusionTerm(mu) \
                  + mu * (u.getFaceGrad().dot(Xhat) * Xhat \
                           + v.getFaceGrad().dot(Xhat) * Yhat).getDivergence() \
                  + lambda * (u.getFaceGrad().dot(Xhat) \
                                  + v.getFaceGrad().dot(Yhat))* 
Xhat).getDivergence()  \
                  + f

       v_eqn = DiffusionTerm(mu) \
                  + mu * (u.getFaceGrad().dot(Yhat) * Xhat \
                           + v.getFaceGrad().dot(Yhat) * Yhat).getDivergence() \
                  + lambda * (u.getFaceGrad().dot(Xhat) \
                                  + 
v.getFaceGrad().dot(Yhat))*Yhat).getDivergence()  \
                  + g

Best Regards,

Joel