塑性返回屈服面的程序，悬赏查错

吃草的鸟

hillyuan

吃草的鸟 发表于 2013-5-2 14:15
对应eml
对应ebl

1. Equation (2) wrong. Third terms in your second part does not exists.
>>You are right. He forgot to add geometry nonlinear terms.

2. Consistent with your equation (2)~(4), you must modify your contitutive equations let stress(zz)=0
>>I think his formulation already took into the consideration that the stress normal to the surface is zero.

3. Equation (5) wrong! You do consider strain(23), strain(31). Therefore stress(23) and stress(31) are not zero.
>>He doesn't consider the transverse shear strain. I think his shell element formulation is based on Kirchoff-Love hypothesis instead of Mindlin-Reissner hypothesis. In his formulation, in addition to the membrane part, he also includes the bending part due to the change in curvature. That's why he probably made a mistake in comparing with Abaqus since most likely Abaqus considers transverse shear components.

4. Equation (9)-(11) define strain and stress inconsistent to that of equation (2)~(6). It is used in text book in material in stength but not in continuum mechanics. The stress defined cannot be used to calculate Misess stress! Therefore
!) Forget {\epsilon}_g and {\sigma}_g
!) Use the strain defined in Equation (7)
!) Use the elastic matrix [G]=[C]
>>He is trying to solve a hard problem without the full understanding of everything. I would suggest he solves the problem in a step-by-step way.
1. Get elastic-only case working. It seems like he already got it done.
2. Get elastic+large displacement working.
3. Get plastic-only case working.
4. Get plastic + large displacement working.

To understand the nonlinear problem, it might be a better idea to set up a simple 1D problem, for example, a bar subject to uniaxial tension. Just use linear element to manually solve a two or three element model. Use this 1D model to play with large displacement and plasticity. After you are comfortable with this 1D model. You are ready to focus on shell element.

hillyuan

Program structure well defined. Good! But don't use GOTO. It is dangerous!

I cannot read your program in detail, but at least following errors
1. Your strain calculation is probably wrong. I am not sure what your deml mean but I cannot find strain(33), which along the thickness of shell, in your six strain component.
2. Do you consider the condition strain(33)=0 when calculate your elastic matrix?
3. where is your simgl_eq2?
and so on.

A suggestion,  if(fp2.gt.1e-2) your calculation wrong, because 理想弹塑性 you can obtain the right answer in one step. Just delete relevent code, it is erro-prone.

吃草的鸟

hillyuan 发表于 2013-5-2 09:51
Program structure well defined. Good! But don't use GOTO. It is dangerous!

I cannot read your progr ...

Thank you very much!

hillyuan

吃草的鸟 发表于 2013-5-2 10:20
Thank you very much!
1、关于strain（33）：这里直接设stress(33)=0. stain（33）=0？or strain(11)+str ...

If you suppose strain(33)=0, then you have 5 strain components, not six. The elastic matrix should be 5*5, not 6*6, right? By the way, you cannot expect strain(33)=0.d0 in your programs (it is 0 in linear elastic analysis if you rightfully programed but not in nonlinear analysis)

Is your sigmlg0 3*3?  Then how about matmul(transpose(sigmlg0),HH), this is 3*3 matrix multiple a 6*6 matrix. what you are doing here?

You need reconsider your formulus from the very beginning.

吃草的鸟

hillyuan 发表于 2013-5-2 11:43
If you suppose strain(33)=0, then you have 5 strain components, not six. The elastic matrix should ...

现在程序进行的是非线性分析，不考虑板厚变化，这里哪不妥呢？               

吃草的鸟

hillyuan 发表于 2013-5-2 11:43
If you suppose strain(33)=0, then you have 5 strain components, not six. The elastic matrix should ...

有问题吗？

hillyuan

吃草的鸟 发表于 2013-5-2 12:00
这里是将应变和应力分解为膜应变和弯曲应变Δε=Δεm+Δεb，这里只是将不为0的独立分量拿出来，因为不 ...

Well. I am wrong, the supposition of shell is that stress(33)=0, not strain(33). In 非线性分析, you cannot obtain stress(33)=0 directly. You need another NR iteration to do it.

We have shear strain but there is no 弯曲应变 in continuum mechanics. By the way, your σbxx,σbyy,σbxy is meaningless. What your strain and stress tensor should be strain(exx, eyy, ezz, exy, 0, 0 ) and sigmlg0={σxx,σyy,0, σxy,0, 0} as you supposed.

hillyuan

吃草的鸟 发表于 2013-5-2 12:05
平面应力下σ（MISES）=transpose{σxx,σyy,σxy} [AB] {σxx,σyy,σxy} ,不对吗？
这里只是将σxx,σy ...

what matmul(matmul(transpose(sigmlg0),HH),sigmlg0) do is just calcualte equilvent stress, write down it directly. It is stupid to do such complicated calculation.

吃草的鸟

hillyuan 发表于 2013-5-2 12:10
Well. I am wrong, the supposition of shell is that stress(33)=0, not strain(33). In 非线性分析, you ...

这里并不是说弯曲应变，只是将应变分成在板厚方向不变和线性变化两部分，这个是为后续计算应用的！这里的塑性应力修正就是修正膜应力和弯曲应力部分！横向剪力不考虑！总的和应变当然就是3*3，应力也是3*3！

吃草的鸟

hillyuan 发表于 2013-5-2 12:20
what matmul(matmul(transpose(sigmlg0),HH),sigmlg0) do is just calcualte equilvent stress, write do ...

这是将应变分开后，比较直接的计算方法吧！

hillyuan

吃草的鸟 发表于 2013-5-2 13:00
这是将应变分开后，比较直接的计算方法吧！

It is just 0.5*(s11*s11+s22*s12+s12*s12) for your problem, I think.

吃草的鸟

hillyuan 发表于 2013-5-2 13:24
It is just 0.5*(s11*s11+s22*s12+s12*s12) for your problem, I think.

对，就是不用求S11，S22，S12，直接用矩阵相乘得到了！这个肯定没问题，弹性时也是这么算的！

hillyuan

吃草的鸟 发表于 2013-5-2 12:59
这里并不是说弯曲应变，只是将应变分成在板厚方向不变和线性变化两部分，这个是为后续计算应用的！这里的 ...

You should not have six non-zero strain component because you don't consider 横向剪力.

Write down your formulus from the beggining!

吃草的鸟

已经解决了，谢谢大家

吃草的鸟

hillyuan 发表于 2013-5-2 14:56
1. Equation (2) wrong. Third terms in your second part does not exists.
2. Consistent with your e ...

Thank you very much!

吃草的鸟

hillyuan 发表于 2013-5-2 14:56
1. Equation (2) wrong. Third terms in your second part does not exists.
>>You are right. He forgot ...

非常感谢二位的帮忙！