4

I solved the steady state dynamic linear elastic model in a solid. My equation is a function of frequency and the strong form is: $$\operatorname{div}(\operatorname{stress}(\vec x, w)) + w^2 \rho u( \vec x,w)=0$$ where $\rho$ is the density.

The BC is $$ \sigma(\vec x,w) \cdot n(x)=T(x,w), \qquad u(x,w)=U_0$$

The physical problem is a plate with dimension of $1 \times 1 \times 0.1$ with a harmonic load on the all the top and also the bottom of the plate is fixed. I solved this as a test with the fixed frequency $(w=200)$.

The weak form of the this equation is: $$\int \sigma(x,w)\cdot \operatorname{grad}(v)\, dv - \rho w^2 \int v(x)\cdot u(x,w) \, dv - \int v\cdot T(x,w)\, ds(2) = 0 $$

where $v$ is test function, $ds(2)$ is the top of the plate.

I compared the results with the model with converged mesh from Abaqus, I checked almost everything in the code the order of the displacement is mostly correct but the signs of the displacements are wrong and also the error is large between the solution of the Abaqus and this code.

I checked the variation form and I can't find any error in it. So the only thing that I suspected to be wrong is the Neumann boundary condition but I can't find any bug in it.

I would be very thankful if someone can help me to figure out what is wrong in my code that I get different answers from the code. Or I am not sure if this is a bug and my code is right.

Bahram
  • 111
  • 3
  • Your code matches PDE+BCs given when assuming U0=0 on bottom, T=f on top and T=0 elsewhere. Did you set f[1] = Magnitude < 0 intentionally? – Jan Blechta May 12 '13 at 18:57
  • I just notice this comment about my post.yes, by setting Magnitude<0 I meant to make the T.ns which make a negative number since the direction of the traction is (-). I am not sure why I am receiving different answers from Abaqus since the code is not complicated. can it be because of the bug in library? – Bahram May 22 '13 at 17:30
  • I would say that direction of displacement seems correct. – Jan Blechta May 24 '13 at 10:44
  • when I changed the sign of the f[1] to a positive number then I saw that the direction of the load change in plot. but when the magnitude is negative the displacement plot make sense since the direction of the displacement in toward down. – Bahram May 24 '13 at 19:51
  • So where's the problem? You said that direction of displacement is different from one obtained by Abaqus but now we agreed that this is correct. – Jan Blechta May 24 '13 at 20:06
  • The direction and values of the displacement for some points doesn't match with Abaqus however with current version the vertical displacement in top of the plate is toward down.. but in the middle the sign for some points change which doesn't make sense. Beside that the absolute value for the displacement for different points doesn't match with Abaqus. – Bahram May 28 '13 at 01:19
  • But as you said the direction of displacement on the top of the plate is correct. – Bahram May 28 '13 at 14:13
  • Sorry, but speaking about displacement in y-direction it is negative in whole domain. – Jan Blechta May 28 '13 at 14:22
  • Thank you very much. sorry, You are right, I found out what was the problem. The first problem was related to convergence and the second problem was related to me extracting the nodal displacement from the Fenics in to a text file. the vertical displacement have less than 3% relative error but the largest error is related to the x-displacement in (0.5,:,0.1) path which in abaqus I have the X-displacment equal to zero but from the Fenics it is not zero. – Bahram May 28 '13 at 18:38
  • I would focus on boundary conditions. Also make sure that same set of elastic parameters is used and correctly converted. I discourage you also from typing int parameters where they're evidently float by purpose - it often leads to hidden programming errors, for example 2/3 = 0. – Jan Blechta May 28 '13 at 18:52
  • Probably you have too coarse mesh so that SubDomain.mark() function can't hit any entire facet. Remember that for given facet all its vertices plus midpoint must meet SubDomain.inside() condition to be marked. See Johan suggestion http://scicomp.stackexchange.com/a/7381 how to debug SubDomains. – Jan Blechta May 28 '13 at 20:15
  • Do I understand correctly that you have the z-displacement equal zero on the entire 1 x 1 bottom surface of the plate? This will be a difficult dynamics problem to solve with any finite element analysis code. Exactly how are you solving it with abaqus-- element type and mesh density? If your objective is to experiment with dynamic elasticity in Fenics, I suggest a different problem. – Bill Greene Apr 09 '15 at 13:00

0 Answers0