| Author | Message | |||||
     Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 24 Registered: 03-2004 |
For displacements that are larger than the thickness of the plate, the middle surface gets streched. According to Timoshenko, "Theory of Elastic Stability", edition 1, 1936, deflection W can be approximated by the following 2-dimensional equations. 1. del2(del2(W)) = (h/D)*[ dyy(F)*dxx(W) + dxx(F)*dyy(W) - 2*dx(dy(F))*dx(dy(W)) ] 2. del2(del2(F)) = E[ dx(dy(W))^2 - dxx(W)^2*dyy(W)^2 ] where F is a scalar stress function of x and y h is thickness D = E*h^3/(12*(1-nu^2)) E is Young's Modulus nu is Poisson's ratio and h is the plate thickness. After solving (1) and (2): The components of stress in the middle plane , Nx, Ny and Nxy are calculated from Nx = h*dyy(F) Ny = h*dxx(F) Nxy = -h*dx(dy(F)) The strains in the middle plane, ex, ey and exy are calculated from ex = (1/(h*E))*( Nx - nu*Ny) ey = (1/(h*E))*(Ny - nu*Nx) exy = (2*(1 + nu)/E)*Nxy The strains are defined in terms of derivatives as follows ex = dx(U) + .5*dx(W)^2 ey = dy(V) + .5*dy(W)^2 exy = dy(U) + dx(V) + dx(W)*dy(W) The boundary conditions are especially perplexing At two opposing ends W = f(y), dx(U) = constant and dx(V) = 0 The sides are free. Any suggestions? I've looked at freeplate and fixed plate. But, I'm still completely baffled and confess that I can't even understand the justification for natural(V) = L*U in fixplate. | |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 25 Registered: 03-2004 |
I've made some progress. It turns out that equation (1) is nothing more than a tortured form of the equations of equilbrium. Timoshenko developed it by manipulating the strain and stress-strain equations. So, I can replace (1) with: (1a) U: div(vector(Nx, Nxy)) = 0 (1b) V: div(vector(Nxy, Ny)) = 0 and define Nx, Ny and Nxy in the definitions section as I would for an ordinary plane stress problem. Equation (2) is still needed because the z-axis displacement, W, appears in the definitions of strain used in (1). So, now I only need to do something with (2) that will get around the high order derivative problem and the stuff on the right side. I could use your trick from "Fixplate" and introduce an intermediate variable with del2(W) = something . But, I don't see how to use "something" to create meaningful boundary conditions, such as value(W) = f(y). | |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 26 Registered: 03-2004 |
Beg your pardon. I should have said equation (2) where I said equation (1) and vice versa. Furthermore (1) can be written in terms of Nx, Ny and Nxy. So, the equations are (1) W: del2(del2(W)) = (1/D)*[ Nx*dxx(W) + Ny*dyy(W) + 2*Nxy*dx(dy(W)) ] (2a) U: div(vector(Nx, Nxy)) = 0 (2b) V: div(vector(Nxy, Ny)) = 0 | |||||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 694 Registered: 06-2003 |
The condition Natural(V)=L*U that we used in the fixplate example was given to us by Said Doss at Lawrence Livermore Lab. I thought it was rather cute. The PDE is delsq(V)+Q=0. At a straight boundary, the tangential curvature is zero, so delsq(V) is the normal curvature. In the absence of Q, this means the slope is constant near the boundary. The BC says that the normal derivative is a very large multiple of U; ergo, a very steep slope if U does not vanish. But the solution cannot rise rapidly; it must curl over and meet the other boundary. The only way to satisfy all these conditions is if U goes to (near) zero, implying a modest normal derivative at the boundary. The bigger L is, the closer U approaches zero. As far as the formulation of your system, I have never done this particular model before, so I don't know how much useful information I can contribute. Clearly, you need two auxilliary variables, one for delsq(W) and one for delsq(F). This means four equations and eight boundary conditions, just as the numbers require. My guess is that a condition Natural(M)=L*(W-f(y)) would be the analog of the fixplate BC, driving W to f(y). | |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 27 Registered: 03-2004 |
Robert, Thanks very much for the suggestions. I finally came to understand the Natural(V) = L*U earlier today. Another way to look at it is that dV/dn is the shear strain and physical considerations require that it cannot become unreasonably large. Your explanation is better. I agree that it is a clever trick and I'm going to use it. The W-f(Y) notion had occurred to me too. I think I've finally got the model working. It's attached.
| |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 28 Registered: 03-2004 |
I see that the attachment didn't make it. Here it is as text. I have found that it's necessary to decrease errlim as low as 0.00001 in some cases. This can make solutions take a bit of time. Maybe you can see a way of making an improvement in accuracy vs speed. ******************** { FIXPLATE_LARGE_DEFLECTION.PDE } { *********************************************** This example considers the bending of a thin rectangular plate under a distributed transverse load producing deflections that are larger than the plate thickness. When a plate undergoes large deflections, the mid-plane is stretched. These stresses interact with the bending stress to limit the deflection. Timoshenko in "Theory of Elastic Stability" describes a way to handle such problems without resorting to a full 3D treatment. The method was originally developed by Theodore von Karman. It uses the 2D plane stress equations in combination with a modied biharmonic plate equation. The equations are as follows div(vector(Nx, Txy)) = 0 div(vector(Txy, Ny)) = 0 del2(del2(W)) = (1/D)*[ Q + Nx*dxx(W) + Ny*dyy(W) + 2*Nxy*dx(dy(W)) ] where Q is the load distribution, D = E*h^3/(12*(1-nu^2)) (Modulus of rigidity) E is Young's Modulus nu is Poisson's ratio Nx, Ny and Nxy are the mid-plane stresses of the plate and h is the plate thickness. The boundary conditions to be imposed depend on the way in which the plate is mounted. Here we consider the case of a clamped boundary, for which U = 0 (x displacement) V = 0 (y displacement) W = 0 (z displacement) dW/dn = 0 FlexPDE cannot directly solve the fourth order equation, but if we define rho = del2(W), then the deflection equation becomes U: del2(W) = rho (rho is equivalent to the sum of the x and y plate curvatures) rho: del2(rho) = Q with the boundary conditions dW/dn = 0 d(rho)/dn = L*W where L is a very large number. Natural(W) = 0 makes the spatial derivative of displacement (slope of plate mid-plane) in the direction normal to the boundary equal to zero. Natural(V) = L*W makes the spatial derivative of total curvature in the direction normal to the boundary into a large number unless U = W. This forces U->0 because d(rho)/dn is the shear in the direction normal to the boundary In this system, d(rho)/dn can only remain bounded if W -> 0, satisfying the value condition on W. The particular problem addressed here is a plate of 16-gauge steel, 8 x 11.2 inches, covering a vacuum chamber, with a pressure. A comparison with FIXPLATE is as follows: Max displacement, edges clamped in all x and y in FIXPLATE_LARGE_DEFLECTION Pressure (Q)------------------FIXPLATE-----------------FIXPLATE_LARGE_DEFLECTION---Maximu m in-plane stress ----0.01 psi----------------------0.000149-----------------0.000149---------------------- --------------<0.025 psi ----0.1-----------------------------0.00149------------------0.00148------------ ---------------------------2.93 ----1.0-----------------------------0.0149--------------------0.00898----------- -----------------------------106 ----14.7---------------------------0.219----------------------0.0394------------ ----------------------------2013 Max displacement, edges clamped only in z plane in FIXPLATE_LARGE_DEFLECTION Pressure (Q) FIXPLATE FIXPLATE_LARGE_DEFLECTION Maximum in-plane stress ----0.01 psi----------------------0.000149-----------------0.000149---------------------- --------------<0.009 psi ----0.1-----------------------------0.00149-------------------0.00149----------- ---------------------------1.34 ----1.0-----------------------------0.0149--------------------0.0136------------ -----------------------------107 ----14.7---------------------------0.219---------------------- 0.055------------------------------------------1351 FIXPLATE results diverge from FIXPLATE_LARGE_DEFLECTION for deflections above a few thousandths of an inch. This is to be expected since the equations used in FIXPLATE are good only for deflections that are significantly less than the plate thickness. Note that the maximum displacements in FIXPLATE_LARGE_DEFLECTION diverge more from FIXPLATE results when the plate edges are clamped in all three directions. This makes sense since FIXPLATE puts no constraints on edge displacements in x or y directions. } title " Plate bending with stretching and large deflection " Select cubic { Use Cubic Basis } autostage = off Variables U(0.1) V(0.1) W(0.1) rho(0.1) Definitions k = staged(1, .1, .01) ! The purpose of this is to see the effect of decreasing errlim Select errlim = k*.001 Definitions xslab = 11.2 yslab = 8 h = 0.0598 {16 ga} L = 1.0e4 E = 29e6 Q = 1 nu = .3 D = E*h**3/(12*(1-nu**2)) {stress coefficients } C11 = E/(1-nu^2) C12 = C11*nu C33 = E/(2*(1+nu)) {Strains} ex = dx(U) ey = dy(V) uy = dy(U) vx = dx(V) wx = dx(W) wy = dy(W) exy = uy + vx + [wx*wy] ux = ex + .5*( wx^2) vy = ey + .5*(wy^2) {Stresses} Nx = C11*ux+C12*(vy) { Stresses } Ny = C11*vy+C12*(ux) Nxy = C33*(exy) {Vectors for divergence} Pgeneral = vector(Nx, Nxy) Qgeneral = vector(Nxy, Ny) gap = 0.2 Initial values U = 0 V = 0 W = 0 rho = 0 Constraints {use to prevent rigid body motion when U and V are not clamped} surf_integral(U) = 0 surf_integral(V) = 0 Equations U: div(Pgeneral) = 0 V: div(Qgeneral) = 0 W: del2(W) = rho ! rho is the total curvature [dxx(W) + dyy(W)] of the mid-plane of the plate rho: del2(rho) = (1/D)*[ Q + Nx*dxx(W) + Ny*dyy(W) + 2*Nxy*dx(dy(W)) ] Boundaries Region 1 start (0,0) natural(U) = 0 ! Use natural(U) = 0 if not clamped in x direction natural(V) = 0 ! Use natural(V) = 0 if not clamped in y direction natural(W) = 0 natural(rho) = L*W ! Clamps edge in z direction line to (xslab,0) to (xslab,yslab) to (0,yslab) to close monitors contour(W) plots contour (W) as "Displacement" elevation(W) from (0,yslab/2) to (xslab,yslab/2) as "Displacement" surface(W) as "Displacement" surface(rho) as "rho" contour(Nx) as "x-axis stress" contour(Ny) as "y_axis stress" contour((Nx*dxx(W) + Ny*dyy(W) + 2*Nxy*dx(dy(W)))/D, Q/d) elevation(W) from (0, 0) to (xslab, 0) elevation(U) from (0, 0) to (xslab, 0) elevation(V) from (0, 0) to (xslab, 0) elevation(Ny) from (xslab/2, 0) to (xslab/2, yslab) elevation(Nx) from (0, yslab/2) to (xslab, yslab/2) summary("") report(val(W, xslab/2, yslab/2)) report(val(Nx, xslab/2, yslab/2)) report(val(Ny, xslab/2, yslab/2)) report(errlim) end 220693428 | |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 29 Registered: 03-2004 |
I would also like to offer some praise for FlexPDE. During this summer I've gone down lots of roads looking for what you see above. A big part of the work was learning enough about plate and shell theory to know what to look for. FlexPDE has not only given me a very cost effective tool for implementing the solution. Its ability to let me experiment with equations and see immediate feedback has been an indispensable part of my learning process. Thanks for a great product. Jerry Brown | |||||
| Peter Cendula (cendo) Junior Member Username: cendo Post Number: 3 Registered: 06-2007 |
Hi, I think you have an error in your script in equation rho: del2(rho) = (1/D)*[ Q + Nx*dxx(W) + Ny*dyy(W) + 2*Nxy*dx(dy(W)) ] there should be rho: del2(rho) = (1/D)*[ Q + h*Nx*dxx(W) + h*Ny*dyy(W) + 2*h*Nxy*dx(dy(W)) ] because you omitted "h". The comparisom with fixplate.pde is not so drastically different for larger Q, e.g. for Q=14.7 the Max displacement, edges clamped only in z plane in FIXPLATE_LARGE_DEFLECTION is 0.14 instead of your 0.055. Please respond, thanks peter | |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 37 Registered: 03-2004 |
Peter, Thanks for the note. I had not noticed it until today. Your're quite right. I discovered the mistake a while ago and should have come back here to alert anyone who might try to use my script. I'm still trying to find a way to use this model for a problem where two edges are clamped ond the other two are free. If the free edge is parallel with the y axis the boundary conditions are dxx(W) + nu*dyy(W) = 0 and dxxx(W) + (2 - nu)*dx(dyy(W)) = 0. Of course, I have to use the "Fixplate" trick for the third derivatives. But, I can't see a way to do it so that both kinds of edge conditions can be defined. Thanks, Jerry | |||||
| Peter Cendula (cendo) Member Username: cendo Post Number: 5 Registered: 06-2007 |
Hi Jerry, thank you for comment. I'm actually also solving problem with two free ends and I first started with the change your posted code, see attached. I enforce free ends by natural conditions on U,V and it looks working. I think in your BC, dxx(W) + nu*dyy(W) = 0 and dxxx(W) + (2 - nu)*dx(dyy(W)) = 0, you mistyped W instead of F (Airy stress function) and +nu instead of -nu. Hopefully it helps you Peter
| |||||
| Jerry Brown (jerrybrown11743) Member Username: jerrybrown11743 Post Number: 38 Registered: 03-2004 |
Peter, Thanks for the feedback. Sorry I haven't responded sooner. Since we're working on the similar problems, I will start checking this thread daily. Two heads are always better than one. Your solution works only if the two free ends are parallel with the y axis. Here is why. The free edge boundary conditions described in my previous post come from Timoshenko's "Theory of Elastic Stability". For a free edge that is parallel with the y-axis the conditions are: Mx = 0 Mxy = 0 Qx = 0 where Mx is the bending moment per unit length of the edge parallel to the y-axis, Qx is the shearing force per unit length and Mxy is the twisting moment per unit length about the x-axis. Mx = -D*(dxx(W) + nu*dyy(W)) Mxy = -Myx = D*(1 - nu)*dx(dy(W)) Qx = dx(Mx) + dy(Myx) Timoshenko's sign conventions are used. Timoshenko explains that the Qx = 0 and Mxy = 0 conditions can be reduced to one: Qx - dy(Mxy) = 0 This and Mx = 0 reduce to the two BCs in my previous post when expressed in terms of derivatives of W. In your problem the free edge is parallel to the x-axis. So, the BCs would be. My = D*(dyy(W) + nu*Dxx(W)) = 0 Qy + dx(Myx)= dyyy(W) + (2 - nu)*dy(dxx(W)) = 0 Your BC of natural(W) = 0 is equivalent to dy(W) = 0. That works provided there are no forces that are trying to twist the plate about the x-axis, because it is consistent with the condition that My = 0. The BC of natural(rho) = 0 is equivalent to Qy = 0 which is correct for a free edge. I think you will get correct answers so long as the applied load, Q, is uniform in the y direction. I'm trying to solve a problem in which the two clamped ends are twisted relative to each other (they are not parallel to the y-axis or to each other). So, dy(W) will not generally be equal to 0. |
