Elasticity problem

Question

Hi
I want to solve a elasticity problem in a 2 dimensional domain. My equation is:

We can write it in this form which is the equation I am trying to solve : (eq.*)

In general

and C is a tensor of rank 4 that can be written in a 66 matrix. But in plane-stress case it reduces to a 33 matrix.

This is a part of my code:

C= Constant(((E11,E12,0),(E12,E22,0),(0,0,E33)))
domain = mshr.Rectangle(Point(0,0), Point(10,5))
mesh = mshr.generate_mesh(domain, 60)
V=VectorFunctionSpace(mesh1, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
a =.....+dt*dt*inner(C*grad(u), grad(v))*dx+.....

The problem is that the C tensor (which is a 33 matrix) in my weak form can not be multiplied by the displacement. Because I only have two dimensions and displacement is a 21 vector . That is why I get "mismatch.shape" error in the line where I have defined "a" according to (eq.*). Could you help me to resolve this issue?

Answer 1

If your medium is isotropic then the strain-stress relations might be simplified to

$$\sigma=\lambda\;tr(\varepsilon)I_d + 2\mu\varepsilon,$$

where $I_d$ is the identity matrix (2x2 in 2D case), $\mu$ and $\lambda$ the Lammé parameters, $tr(\varepsilon)$ the first invariant of the strain tensor and

$$\varepsilon=\frac{1}{2}(grad(u) + grad(u)^T).$$

Of this manner, your equation of motion can be written as

$$div(\lambda\;tr(\varepsilon)I_d + 2\mu\varepsilon) + f = \rho\ddot{u}.$$

About the implementation, maybe is recomedable for you take a look to the sections 30.3.1 and 30.6.4 of the fenics book and to the elastodynamics demo.

Regards!

Comments

I have already done that for an isotropic medium but in this case my material is anisotrop. The problem is that I do not know how to define tensor C in which I can use in eq.*
Any idea?

---

Hello,
Since your C is a fourth order compliance tensor (for 2D - (2, 2, 2, 2)), you need to define it,as like(Tensor product expressions ):

C_numpy = zeros((2, 2, 2, 2))
# fill non zeroes elements of tensor
C_numpy[..,..,..,..] = ...
# and finally
sigma = as_tensor(C_numpy[i,j,k,l]*eps[k,l],(i,j))

Or, you can compress tensor notation to Voight one and than multiply all like matrices(elasticity - Voight notation):

def eps(u):
        """ Returns a vector of strains of size (3,1) in the Voigt notation
        layout {eps_xx, eps_yy, gamma_xy} where gamma_xy = 2*eps_xy"""
        return df.as_vector([u[i].dx(i) for i in range(2)] +
                            [u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])

---

the integral part of $div(\sigma(u))$ it's ok. If v is a test function in a apropiate space then the weak form of your problem is

$$\int_{\Omega} \rho\ddot{u}\cdot v\;d\Omega + \int_{\Omega}\sigma(u):\varepsilon(v)\;d\Omega = \int_{\Gamma_N}g\cdot v \;d\Gamma + \int_{\Omega}f\cdot v \;d\Omega,$$

where g are tractions in the Newmann boundary (the second term at the left side of the equality and the first to the other side corresponds to $div(\sigma)$ part).

If $\varepsilon$ was defined previously as an user defined function (as in the second link) and using the Voight notation then $\int_{\Omega}\sigma(u):\varepsilon(v)\;d\Omega = \int_{\Omega}C\varepsilon(u):\varepsilon(v)\;d\Omega$.

---

The weak form of my last reply it's ok (see page 70 of the book Full Seismic Waveform Inversion for some references).

---

Hernan
Thank you so much for your help. It was really helpful.