Your reasoning regarding the time step size makes sense. Yep, I changed to a full-scale model instead of the 1/8th just to make sense of the strange behavior. I'll change back. The main change in the script is the calculation of stress, which I now do in the time loop:
line 387: stress = pb.evaluate('de_cauchy_stress.2.Omega(quasi.D, u)', copy_materials=False)
This is then passed to update_stiffness and repeated to get the stresses at the qp's:
stress=nm.repeat(stress,4)
stress.shape=strain.shape
szz1 = avgqp(stress[..., 2:3, :], problem)I'm doing this instead of calculating stress via: szz = ezz * smix
Not sure of the relevance of "copy_materials=False" in the stress function. I also opted to only calculate the D matrix instead of using Lame's parameters, which seems to speed things up.
thx
On Tue, Feb 8, 2011 at 2:41 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
I do not have a definite answer, only some remarks that might be related to the problem:
When I run the input srs_03.py, I get:
Loading rate (mm/s): -0.02117 warning: (almost) singular matrix! (estimated cond. number: 4.38e+14) Initing ... 6.66666666667 0.000905933418981
- the singular matrix means (provided the Lame's parameters are ok), that the EBCs are not sufficient to prevent rigid body motions. The current EBCs just fix the 'z' axis motion, but you need to fix at least one other node completely and another node to prevent rotations around the 'z' axis.
This alone might cause the problems, as who knows what the linear solver does when it solves with a singular matrix.
Otherwise to your problem - you keep the loading rate the same, but change the displacement increment size in a single time step, right?
There is always a trade-off when approximating (time) derivatives with differences - if the step is too large, the derivative is not accurate, while if the step is too small, one gets big rounding-off error due to floating point arithmetics. So there is an optimal step, but I cannot tell you what it is for your problem...
The scheme we have is sort of ad-hoc devised, maybe a mathematician should look at the equations in future...
As for the stress computation, how did you change the way of its calculation?
So you see, I have more questions than answers :)
r.
On Mon, 7 Feb 2011, Andre Smit wrote:
Changing the way the stress is calculated gives more realistic results -
my assumption regarding Hooke's law was perhaps invalid.
On Mon, Feb 7, 2011 at 3:37 PM, Andre Smit <freev...@gmail.com> wrote: Robert
I'm having problems relating the model's output to the actual test output. One issue that is confusing is that if I reduce the displacement interval from 0.001 to 0.0001, say and the time during this interval to ensure a consistent loading rate, then the results change significantly, both the maximum load at failure and the displacement at maximum load. This is unexpected. I can't think what the problem can be - any ideas on your side?On Wed, Jan 26, 2011 at 10:50 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote:
Ok, interesting. Well, I cannot help you much at this point :) You might want to try a finer mesh, just to compare with the current one...r.
On Wed, 26 Jan 2011, Andre Smit wrote:
I think the iterations are indirectly sorting this out. In fact I found that the initial stiffness(0) doesn't influence the results at all. I've attached results of compression (comp) and tension (tens) tests at rate of 1mm/sec as well as the script in its current state. Based on lab tests, comparative compression failure at this rate occurs at 6.6 N/mm2 whereas tension failure occurs at 3.0 N/mm2. So the script is accurately predicting compression failure but under-estimating tension failure. Still exploring ... a On Wed, Jan 26, 2011 at 10:19 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote: Well, this flag is used essentially only in the generic solver used by simple.py. In the strain rate solver script we solve every time step regardless this flag. The first time step is different, however, as strain equals strain0 (no previous strain defined), and so initial dstrain is zero. Maybe this causes some problems? r. On Wed, 26 Jan 2011, Andre Smit wrote: Alright, that fixes it! Now I understand the reason for quasistatic! In the non-linear analyses, the time stepper is set programatically - is is possible to force solving for the first time step here as well as in: for ii, disp in ds: #my_output(dformat % (ii + 1, ds.n_step, disp)) force0 = 0.0 pb.ebcs['Load'].dofs['u.2'] = disp pb.ts.is_quasistatic = True # <<<=========== Force quasistatic pb.time_update(ts) On Wed, Jan 26, 2011 at 9:57 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote: So the problem is the single dot outside the linear curve, right? What if you add 'quasistatic' : True, to the time-stepper options? As it is now, the first time step is not solved but taken from the initial conditions (unspecified = zero) and boundary conditions (nonzero!), so it is not in equlibrium! r. On Wed, 26 Jan 2011, Andre Smit wrote: Sorry - forgot to attach the figure - your point is taken re the nodal residual. Forces are calculated and summed on each of the nodes in the Top region. I checked this a while back and plotted in Paraview - the forces are zero at the other nodes within the model but equal and opposite in direction to the Top at the corresponding Bottom nodes of the model (as you'd expect). a On Wed, Jan 26, 2011 at 9:26 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote: Where can I see the line? maybe you forgot to attach a figure? As for the forces, it might be better to compute them using a (inner) surface integral of stress instead of the nodal rezidual. Btw. you look at the nodal force in the first node of the Top region only, right? What are the values in the other nodes? r. On Wed, 26 Jan 2011, Andre Smit wrote: Thanks Robert! I've attached a modification that shows the strain/stress output for the elastic case of the cylinder under compression. The slope of the line is Young's modulus as you'd expect. As with our non-linear analyses, the forces calculated after the first time step appear to be too high - not sure what the reason for this is. On Wed, Jan 26, 2011 at 2:25 AM, Robert Cimrman <cimr...@ntc.zcu.cz> wrote: There was some old code in the stress_strain() function. The exception was caused by putting directly the traction function into the EBC definition, instead of its name. The attached file should work. r. -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en. -- Andre -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en. -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en. -- Andre -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en. -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en. -- Andre -- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
-- Andre
-- Andre
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.
-- Andre