View previous topic :: View next topic 
Author 
Message 
viroxa
Joined: 28 Jul 2017 Posts: 78

Posted: Fri Oct 13, 2017 11:40 pm Post subject: Detecting Semantic Errors 


Hello, Everyone,
I know that this is an extremely general question, but I'm desperate for some input because I've been turning in circles the last two weeks:
What's your approach to finding the source of any semantic errors?
I've been working on a small FEprogram and now that it's supposedly finished, the results are wrong. This is especially frustrating since I know from comparison with a friends work and literature that the element stiffness array is assembled correctly.
The bits that assemble the global stiffness array and carry out the Choleskyfactorization and forward / backwardsubstitution are taken from "Programming the Finite Element Method" by Smith, Griffiths, Margetts (fifth edition) and are therefore assumed to be correct (at least, after four editions, one should expect them to be).
The load array can't be the source either, because at the moment it only contains 3 values (vertical loads at 3 nodes) for testing.
Considering
K*f = u
K and f are correct, the solver should work fine, but u is still wrong.
Any advice is greatly appreciated!
Have a nice weekend, everyone! 

Back to top 


mecej4
Joined: 31 Oct 2006 Posts: 976

Posted: Sat Oct 14, 2017 12:43 am Post subject: 


Something is wrong with your description. If K is the stiffness matrix, u is the vector of displacements and f is the vector of loads, you should have
K.u = f
Given f, you would solve these equations for the displacements u. 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 1979 Location: Sydney

Posted: Sat Oct 14, 2017 2:16 am Post subject: 


There can be a long way between "supposedly finished" and correct.
If you have separate solver routines to:
reduce K to upper diagonal
reduce f and back substitute f to obtain u
An approach would be to provide a simple version of K and f and see if the solver routines solve a known problem. ie bypass the generation of K and f and use a simple "known" matrix with a known solution.
Also: error = f  K.u is a check of the solution, so find the maximum value in abs(error). This is an easy/quick problem to solve, either by saving a copy of K and f, before reduction or regenerating K and f, then doing K.u as error_vec = f  sum_of_vectors ( K_column_i * u(I) ).
The idea is to try and test the various stages of the FE process.
If it is any help, in my first FE program (in 1975), I set the property id of each element, outside the element DO loop, so every element had the same id. It took me days to locate this error, when I was sure I had coded the program correctly. Debugging was a bit harder when using card decks and batched runs.
John 

Back to top 


viroxa
Joined: 28 Jul 2017 Posts: 78

Posted: Sat Oct 14, 2017 11:46 am Post subject: 


Sorry, mecej4, K*u = f was what I meant to write.
Thanks for the suggestions, John. I'll give them a try. 

Back to top 


JohnSilver
Joined: 30 Jul 2013 Posts: 910 Location: Aerospace Valley

Posted: Sun Oct 15, 2017 5:00 pm Post subject: 


viroxa ... let's go back to some basics first ...
how do you know the results are wrong ?
what are you comparing against ?
have you taken the code directly from the book and if so are you 100% sure it's correctly copied ? (this can take anywhere between 2 and 100+ comparisons to find an error if it's there !)
if it's based on youtìr coding of the theory written in the book then ditto bur it wuìill probably take even longer to find !
An alternative  start with a 2x2 or 3s3 simple inìvented matrix and solve that, then compare with a solution done in excel ??? 

Back to top 


viroxa
Joined: 28 Jul 2017 Posts: 78

Posted: Mon Oct 16, 2017 12:35 am Post subject: 


John, I know that the results are wrong, because I'm using a cantilevered plate which should only bend around the yaxis. But the free end is also warped.
The author of the book offers Fortran files, so I just copy/pasted the code.
In the meantime I've figured out that the error was hidden in the assembly of the nodal freedom array. I haven't yet fixed it completely, but I'm getting there. 

Back to top 


JohnCampbell
Joined: 16 Feb 2006 Posts: 1979 Location: Sydney

Posted: Mon Oct 16, 2017 2:42 am Post subject: 


Quote:  I'm using a cantilevered plate which should only bend around the yaxis. 
It would depend on the mesh orientation and how the load is applied. If it is an edge load, with say a 2x2 plate mesh, then the loads on the 3node end edge should be 1:2:1, not 1:1:1 (which would curl the end of the plate). I am assuming 4node plate elements;, and quite a lot else!
You should test against known solutions. 

Back to top 


LitusSaxonicum
Joined: 23 Aug 2005 Posts: 1844 Location: Yateley, Hants, UK

Posted: Mon Oct 16, 2017 12:46 pm Post subject: 


Viroxa,
I'm afraid that your structural appreciation is at fault. The effect you are getting is known as 'anticlastic bending' and results from Poisson's Ratio. In a cantilevered plate (say) the top of the plate in a spanwise direction is in tension, and the bottom is in compression, i.e. it is 'hogging'. The stresses induced at right angles to this spanwise direction are tension at the top and compression at the bottom. This will result in an upward curl of the plate in a 'sagging' sense, i.e. the transverse curvature is 'sagging' with the plate edges rising.
Of course along the fixed edge (cantilever root), transverse curvature won't show as it is fixed, but it will show in the distribution of reactions.
You can demonstrate this to yourself if you flex something in your fingers with a low E and high Poisson's Ratio, like an ordinary rubber eraser. You will see the transverse curvature right away.
A correct FE formulation will show this effect, simple bending theory won't.
You will also see it in box girder bridges.
Beware too of bending of nonsymmetric sections, which rotate as well as deflecting.
Best regards
Eddie 

Back to top 


