Hi Adnan Shahriar ,

Usually, if your algorithm is correct, the reasons for this error are mesh generation, boundary conditions, load loading methods, and material parameter model problems. You can check them one by one.

Zhengwei Fan Thanks!

I have checked all except loading method (still haven'e done surface integration, so go with concentrated load).

2 things I found are,

1. If I increase my mesh density (keeping Wedge elements angle proper), the result actually is converging to Abaqus result!

2. Also, for Abaqus, the fixed BC gives negligible deformation (in order of -34). This may be negligible and round off error, but signals they use different technique to solve the system of linear equations. In my code, it is impossible to get non 0 for fixed boundary, no matter how small it is!

3. Also, I have done 1 element test, the error is in same order!

Idk, may be they are doing something to get more accurate result with less mesh density or something else! But I will try more to find out if there is any error in my code.

Hello Adnan Shahriar,

I had the same problem when comparing my FEM Matlab code with Abaqus.

When you apply concentrated forces on 3D solid element the shear locking control method adopted by Abaqus will have a strong impact on the displacement response (I presume that you didn't implement that in your code or not the same of Abaqus). If you apply a uniform surfacique force the displacement will be less impacted. I personally got a 0.1% discrepancy on such load cases.

I suggest you to validate your code using a uniform surfacique force.

If you wish to have the exact same displacements on the concentrated force (Good luck! XD) you can give a look to some literature techniques for shear locking control and then implement one of them but the result is not guaranteed being Abaqus implementation not described in the documentation (for obvious reasons).

Best regards

Simone

Thanks a lot Simone Coniglio! I thought using the full integration negate shear locking correction in Abaqus. I will adopt the surface integration within days and match the results.

My model does not have any shear locking control. But one way or another I need to adopt higher order elements as it seems to face significant bending.

Adnan Shahriar What is the exact test case with which you get a discrepancy of 20% in the 3D case and 3% in the plane strain case?

Hi Philipp Gampe

For 2D case, 4 node quad gives some discrepancy but 8 node matches perfectly!

For 3D, its a hollow structure with fixed end condition at base and concentrated load at top.

Hi Adnan Shahriar,

if you don't use a fixed end condition you can get an analytical solution for this case, assuming an elastic material and a static problem. This allows you to compare not only your results to the ABAQUS result but both FEA results to an analytical result.

In the 2D case you could use the the plane strain case of an infinite thick-walled hollow cylinder under internal pressure. This allows you to test your FEA code with 2D plane strain elements, and 3D elements with the appropriate boundary conditions to get the plane strain. It's a great way to test your code because you get a variation of the strain in the radial direction.

For both cases you can find the analytical solution in the book "Applied Mechanics of Solids" by Allen F. Bower. Luckily, the contents of his book can be found on his website: http://solidmechanics.org/Text/Chapter4_1/Chapter4_1.php

When you want to compare the solutions for linear 2D elements you need to look into the manual of your FE program. Commercial programs like to use a selectively reduced integration to avoid locking phenomena. Bower actually mentions briefly mentions that in his book. http://solidmechanics.org/Text/Chapter8_6/Chapter8_6.php

"In many commercial codes, the `fully integrated’ elements actually use selective reduced integration." This can lead to a discrepancy in the results.

The picture of your test case looks nice. What tool did you use for post processing?

Thanks a lot Philipp Gampe. I will go through the links and look for the analytical solution for 3D and already have good result for 2D models with analytical.

I have used Matlab surface plot and I don't think its really a good idea! I actually have crossing or overlapping surfaces colored based on y coordinate (can be displacement or stress but require putting data into 3D Meshgrid format! I don't like doing that! Interpolation or triangulation takes huge amount of time).

Adnan Shahriar I used patch matlab function you can provide color fields on faces or on node and interpolation is carrried out internally.

best regards

Simone

Thanks Simone Coniglio . Patch is a great tool!

Meanwhile, did the Surface traction last week, the force vector matched exactly with ABAQUS, but the deformation error is still in same order.

Also, I found that, ABAQUS is probably use one integration point (even if Reduced Integration unchecked) when determining Nodal force from temperature distribution hence get same magnitude vector (Reaction Force) in 4 nodes, while 4 IP change the vectors drastically. The results are below.

1. Uniaxial temperature distribution

2. Reaction Forces for only one 4 noded full integration element

3. Force vector for same type element in Matlab

Hi Adnan Shahriar,

I don't think that ABAQUS uses a reduced integration for to determine thermal stresses and thus, the resulting forces. That just doesn't make sense. But I would need to test it myself to check that.

You wrote about a "uniaxial temperature distribution", but the vector for the reaction forces in the screenshot from ABAQUS to me look more like you have a uniform temperature distribution. The reaction forces in the ABAQUS solution and the forces in your solution looks different.

Hello again,

You can select the integration type in abaqus...

Check your input file you should have S4 elements if you have S4R means you didn't selected a full integration scheme.

Best regards

Simone

Thanks Philipp Gampe .

The temperature is linearly increasing along y axis, but constant along x following expression T=20*Y. I have used Fixed boundary condition for the whole domain to extract Reaction force, which should be (-) Nodal force vectors. Hence although the vector is different, the nodal force is just negative of that.

My code generates the nodal force and I can directly extract that. So, it looks different. But the idea is the same.

Hi Simone Coniglio . Yes the element type is S4. So it should be full integration.

Hi Adnan Shahriar,

thanks for clarifying that. The idea is the same, but how about the values?

From the ABAQUS screenshot it seems that the values for the reaction forces for each node are the same which is not the case for the your MATLAB solution. Of course, they should not be the same. How do the stress and strain field look for your case?

Hi Philipp Gampe . Sorry for this delay!

You can see the results in the attached files. The temperature has been set as linearly decreasing from 10 to 0 from left to right and lower part has set as Fixed. The Abaqus result for 4 node case showed constant temperature and expansion. But the Matlab result with 4 node and Abaqus with 8 node showed different displacement!

Dear Adnan Shahriar

in solid mechanics (Not thermal), the response of your FEM code depends entirely on how the mass matrix was constructed, the ways like lump-method or ... also the amplitude of the applied force/Displacement is important, And these wrongs usually cause discrepancies between in-houses FEM codes and some commercial software like ABAQUS solutions. finally, if you want, you can send your code for debugging.