Why Does My Simulation Not Show Lateral Expansion in a Compressed 2D Square?

[Will Flannery]

For a while I've been working on applying the finite difference method to various problems and am now trying to get a model of stress and strain working. I wrote a sim a few weeks ago and I got one example to run and thought it was good. Now, I've tried it on another example and I see that something is seriously wrong with the sim and my understanding of stress and strain.

The problem is that if I model a 2D square with a uniform force pressing down on the top side, the thing compresses like a spring but doesn't expand laterally, even though the Poisson's ratio of the material is non-zero.

Letting u and v be the displacement functions, the computed solution is that v is linear, like the compression of a spring, from the bottom to the top edge. And, u is 0 everywhere. The thing is, that these functions satisfy the differential equation model, in which all the derivatives are 2nd order and all the 2nd order derivatives of u and v are 0, so, it satisfies the model, but it ain't right, as there is no lateral displacement.

I have a web page I created with the details !

www.berkeleyscience.com/stress.htm

The thing is that the solution has x-direction stress as Hooke's law says that that x stress = A*xstrain + Av*ystrain, and ystrain is not 0, but there is no x displacement associated with the xstress, so ... I think this is where the problem is ..

The beauty of FDM is that it's a cookbook thing more or less, and the error isn't there as the sim solution satisfies the differential equation model, so the error has to be in the differential equations, which seem OK, or the boundary conditions which also seem OK.

Any insights into the matter will be appreciated.

 

[Chestermiller]

U being a linear function will satisfy the differential equation also. The lateral boundary condition is that the lateral stress is equal to zero. This will allow you to determine the lateral strain in terms of the compressive strain.

Chet

 

[AlephZero]

I haven't tried to follow everything on your web page but the fact that it doesn't seem to say anything about boundary conditions seemed a significant omission.

As Chestermiller said, if you haven't satisfied the boundary conditions everywhere (and don't forget that applying a zero external force is still a boundary condition!) it's not surprising that the solution is wrong.

Personally I wouldn't draw any conclusions from the fact that your first test "looks right" unless you have some independent results to compare it with. The most dangerous type of errors are those that "look right" but are actually 20% wrong. (If it's 200% wrong you would probably notice, and if it's 2% wrong you might not care - but 20% wrong is in the danger zone.)

 

[Will Flannery]

The boundary conditions are u and v are 0 along the bottom, unconstrained on the sides (i.e. the derivative wrt x is 0*** the homogeneous Neumann b.c.), u is unconstrained on the top and the derivative of v on the top is -1/A so that the stress is A*(-1/A) = -1.

Of course the first result was misleading - I had to work so hard to get something that 'looked right' that I celebrated early.

The thing is not only that the solution is wrong, it is that it satisfies the differential equation model, which means the model is wrong or the boundary conditions are wrong. I know the referenced equations are right as I got them from a popular textbook. The derivation of the model was some simple substitutions, which I've rechecked 20 times. The boundary conditions seem OK. Whazzzz going on ? I'll start experimenting tomorrow.

*** Yes, this seems suspicious. What other b.c. are possible? If u is linear in x, then specifying the derivative as 0 on one boundary (not to mention two) would make the sim favor a constant solution ... which isn't good ... However the real effect of this BC is just that the boundary value is copied from it's interior neighbor (hence the derivative is 0), which is why I call it the unconstrained BC.

 

[AlephZero]

The boundary conditions are u and v are 0 along the bottom, unconstrained on the sides (i.e. the derivative wrt x is 0*** the homogeneous Neumann b.c.)
...
*** Yes, this seems suspicious.

I agree it looks suspicious, as in "wrong."

If there are no external loads on the left edge, the boundary conditions are ##\sigma_{xx} = 0## and ##\tau_{xy} = \tau_{yx} = 0##. ##\sigma_{yy}## can be anything.

You can transform those into boundary conditions on the strains, and therefore on the derivatives of the displacements. After doing that, you might be able to guess why people use finite elements (formulated using virtual work) in preference to finite differences for stress analysis

We don't know what your final objective is here, but if you want to write a finite difference program for 2-D stress analysis, an easier way (once you have understood the theory!) would be using Airy stress functions: http://en.wikiversity.org/wiki/Airy_stress_function

 

[Chestermiller]

I really don't understand the problem you are solving yet. Are you saying that you have a square sheet of side s and thickness h, and you are compressing it in the thickness direction? Or are you saying you have an infinitely long square cylinder of side s, and you are compressing between two of the opposite sides?

Also, I don't understand how the stress can be -1, with no units, and how, if the strain is -1/A, the stress is obtained by multiplying by A.

Chet

 

[Will Flannery]

Yep, the boundary conditions are the problem. Trying to code them both presents another problem, sigmaxx = 0 determines one value for u on the boundary, and txy=0 another. I've encoded just the first, sigmaxx = 0 and things look better, the thing is expanding laterally, I'm amazed at the difference. But, how to code the txy = 0. I'm 'solving' the discrete problem two ways, by iteration and using a linear equation solver; you only get one new value of u per iteration and I've have two equations if I include the eq for txy = 0. I'd have the same problem trying to include an extra equation (for txy = 0) in the linear equation solver.

Chet, the problem I'm trying to solve is a thin square 2-D sheet, there is a 2D plane stress form of Hooke's law that is applicable. The A is my own shorthand notation for E/(1-v*v) given on the web page.

Edit - I put the boundary value for v into play and was able to code both boundary conditions, and the new output 'looks good'. So, I'm celebrating again ! I put the new output on the web page.

My motivation is pedagogical - engineering and mathematical physics is to a large extent applied differential equations. But I didn't learn this till I was a post doc (math and comp sci)! The basic paradigm is to model a process using differential equations and then to solve the DEs. Solving the DEs has been the rub, as most are unsolvable, so the educational system relegates them to the outback. Enter 'computational calculus', I (actually one of my critics) coined the phrase and I discovered one other fellow, Claes Johnson in Norway, who has been saying the same thing for years, teach computational calculus first. So, he's written a few books, and his thing is the FEM so maybe 6 months ago I started studying it. I wanted to do some interesting projects, and did a few, but then latched onto computational fluid dynamics and flow over a backward step as my holy grail. I discovered that computational fluid dynamicists prefer the FVM, but pressed onward. Finally I discovered The 12 Step Program to Navier-Stokes on YouTube. Fantastic. Using the FDM presented in this course I was able to get a sim for the backward step. And I came to the opinion that the FEM introduces all sorts of machinery that is necessary for real problems with something other than rectangular flat sheets, but that obscures simple problems, and that FDM is a 100% transparent almost cookbook procedure start to finish. So, pedagogically, I think the FDM is the way to begin, and my goal is to get some interesting demonstrations of how processes are modeled and analyzed using FDM.

 

[Chestermiller]

Chet, the problem I'm trying to solve is a thin square 2-D sheet, there is a 2D plane stress form of Hooke's law that is applicable. The A is my own shorthand notation for E/(1-v*v) given on the web page.

I see what your difficulty is now.

With the top and bottom boundary conditions you use, this is the same problem (by symmetry) as if the sheet were twice as thick and u and v were zero at both the top and bottom of the sheet.

Are you considering this as a true 3D problem (with u and v varying with z at the edges, or are you treating it as a pure 2D problem (with u and v not required to satisfy the zero traction boundary condition at the edges)? If you are solving it as a 2D approximation, then you don't need to consider the zero traction boundary condition, and u can be taken as zero at the edges (the tractions at the edges will not be zero). However, if you are trying to determine the full 3D solution to the problem, then you are dealing with a "boundary layer" situation near the edges. In this situation, the solution is virtually indistinguishable from the 2D solution throughout nearly the entire area of the square, and the stresses and compressive force you get are very good approximations. However, near the edges, there is a region on the order of one or two sheet thicknesses from the edges in which the displacements and stresses are varying very rapidly, from the 2D solution inboard of the edge to the zero traction boundary condition at the edge.

If you are solving the problem in 3D, you need to include these boundary layer regions in your analysis. To do this with a finite difference code might be computationally expensive because you will have to use a fine grid throughout the region just to get accuracy on the solution at the edges. This is where a FEM code might be better because you have control over the grid.

If you are solving it as a 2D problem, then you don't need to satisfy the zero traction BC at the edges, and you already have your solution (with u zero at the edge).

Chet