找回密码
 注册
Simdroid-非首页
楼主: yygao56

拓扑优化经典99行程序解读

[复制链接]
发表于 2011-9-20 16:52:55 | 显示全部楼层 来自 江苏南京
shalldy 发表于 2011-9-2 11:01
想请教楼主,99行代码中的敏度推导是不是这样的。以前一直对敏度前面的一个负号耿耿于怀,现在推了一遍,还 ...

倒数第二行是矩阵运算得到应该很自然的吧
回复 不支持

使用道具 举报

发表于 2011-10-15 01:44:10 | 显示全部楼层 来自 陕西西安
Simdroid开发平台
正在看99行,多谢楼主的科普。
回复 不支持

使用道具 举报

发表于 2011-10-26 21:45:13 | 显示全部楼层 来自 江苏南京
都太强了,我只会用软件瞎算
回复 不支持

使用道具 举报

发表于 2011-11-9 20:49:12 | 显示全部楼层 来自 吉林长春
顶楼主,大爱特爱!
回复 不支持

使用道具 举报

发表于 2011-11-22 14:38:37 | 显示全部楼层 来自 北京
本帖最后由 fantasy93404 于 2011-11-22 15:01 编辑

楼主,你好,首先感谢你的分享,我看了你传得pdf,其中有篇参考文献没有找到,请问你那有吗?它是Bends?e, M.P. 1995: Optimization of structural topology,shape and material . Berlin, Heidelberg, New York: Springer
Zhou, M.; Rozvany, G.I.N. 1991: The COC algorithm, part II: Topological, geometry and generalized shape optimization. Comp. Meth. Appl. Mech. Engrng. 89, 197–224这两篇我们学校都没有pdf下载。
回复 不支持

使用道具 举报

发表于 2011-11-23 15:08:54 | 显示全部楼层 来自 北京
等不及了,花钱买了份,有需要的找我
回复 不支持

使用道具 举报

发表于 2011-12-6 18:56:03 | 显示全部楼层 来自 北京
我推了一下平面矩形单元的刚度矩阵,好复杂,结果如下:
K=4ab*
[                                             (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1)]
[   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1)]
[                         (e*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1),   (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1)]
[ - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1)]
[     (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1)]
[   (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) + x/(4*a*b))^2)/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1)]
[                       - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1),   (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1)]
[   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) + (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                         (e*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1) + (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1),   (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) - y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),     (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) - y/(4*a*b))*(1/(4*a) + y/(4*a*b)))/(h^2 - 1), - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1) - (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                       - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))*(1/(4*b) + x/(4*a*b)))/(h^2 - 1),   (e*h*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1) - (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))*(1/(4*b) - x/(4*a*b)))/(h^2 - 1),                                             (e*(h/2 - 1/2)*(1/(4*a) + y/(4*a*b))^2)/(h^2 - 1) - (e*(1/(4*b) - x/(4*a*b))^2)/(h^2 - 1)]
回复 不支持

使用道具 举报

发表于 2011-12-6 18:56:59 | 显示全部楼层 来自 北京
是用matlab推得,其中e是弹性模量,h是泊松比,ab是半个 边长
回复 不支持

使用道具 举报

 楼主| 发表于 2011-12-7 09:28:41 | 显示全部楼层 来自 湖北武汉
fantasy93404 发表于 2011-12-6 18:56
是用matlab推得,其中e是弹性模量,h是泊松比,ab是半个 边长

没这么复杂吧~~~可以参考一下《matlab有限元结构动力学分析与工程应用》
回复 不支持

使用道具 举报

发表于 2011-12-7 10:00:31 | 显示全部楼层 来自 北京
本帖最后由 fantasy93404 于 2011-12-7 13:30 编辑
yygao56 发表于 2011-12-7 09:28
没这么复杂吧~~~可以参考一下《matlab有限元结构动力学分析与工程应用》


我是根据李人宪《有限元法基础》编写的,书中有三角形的单元刚度矩阵,但是对于矩形的单元刚度矩阵,书中只是提到类似,所以。。。
我利用matlab算了一遍得到了上述结果
我看了你说的这本书219页,对比一下,我的做法更一般点,我的矩形单元是
3ab
4a-b
2-ab
1-a-b

《matlab有限元结构动力学分析与工程应用》中将其正则化了,使得a=1,b=1,这样更简单了点


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×

评分

1

查看全部评分

回复 不支持

使用道具 举报

发表于 2011-12-7 19:01:13 | 显示全部楼层 来自 湖北宜昌
本帖最后由 shalldy 于 2011-12-7 19:04 编辑
fantasy93404 发表于 2011-12-7 10:00
我是根据李人宪《有限元法基础》编写的,书中有三角形的单元刚度矩阵,但是对于矩形的单元刚度矩阵,书中 ...

function KE = MechanicalRecEleStiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4)
E=1;NU=0.3;h=1;x1=0;y1=0;x2=1;y2=0;x3=1;y3=1;x4=0;y4=1;
syms s t;
N1 = (1-s)*(1-t)/4;%interpolation functions for 4-node rectangular elements
N2 = (1+s)*(1-t)/4;
N3 = (1+s)*(1+t)/4;
N4 = (1-s)*(1+t)/4;
x = N1*x1 + N2*x2 + N3*x3 + N4*x4;
y = N1*y1 + N2*y2 + N3*y3 + N4*y4;
xs = (-x1+x2+x3-x4)/4+((x1-x2+x3-x4)/4)*t;%differentiate x with repect to s
xt = (-x1-x2+x3+x4)/4+((x1-x2+x3-x4)/4)*s;
ys = (-y1+y2+y3-y4)/4+((y1-y2+y3-y4)/4)*t;
yt = (-y1-y2+y3+y4)/4+((y1-y2+y3-y4)/4)*s;
N1s = -1/4+1/4*t;%differentiate N1 with repect to s
N2s =  1/4-1/4*t;
N3s =  1/4+1/4*t;
N4s = -1/4-1/4*t;
N1t = -1/4+1/4*s;
N2t = -1/4-1/4*s;
N3t =  1/4+1/4*s;
N4t =  1/4-1/4*s;
N1x = yt*N1s - ys*N1t;%differentiate N1 with repect to X
N2x = yt*N2s - ys*N2t;
N3x = yt*N3s - ys*N3t;
N4x = yt*N4s - ys*N4t;
N1y = xs*N1t - xt*N1s;
N2y = xs*N2t - xt*N2s;
N3y = xs*N3t - xt*N3s;
N4y = xs*N4t - xt*N4s;
J = xs*yt - ys*xt;%determinant of Jacobi matrix
B=[ N1x 0   N2x 0   N3x 0   N4x 0;
    0   N1y 0   N2y 0   N3y 0   N4y;
    N1y N1x N2y N2x N3y N3x N4y N4x]/J;
D = (E/(1-NU*NU))*[1, NU, 0 ; NU, 1, 0 ; 0, 0, (1-NU)/2];
BDB = h*B'*D*B*J;
KE = double(int(int(BDB, t, -1, 1), s, -1, 1));

我按照公式推导编的一个刚度矩阵求解的程序,算起来感觉很方便的,矩形单元的四个点的坐标可以任意赋值,不过没用高斯积分,你试试看。

评分

1

查看全部评分

回复 不支持

使用道具 举报

 楼主| 发表于 2011-12-8 09:31:38 | 显示全部楼层 来自 湖北武汉
shalldy 发表于 2011-12-7 19:01
function KE = MechanicalRecEleStiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4)
E=1;NU=0.3;h=1;x1=0;y1=0;x2 ...

挺给力的,给加点分
回复 不支持

使用道具 举报

发表于 2011-12-8 16:22:21 | 显示全部楼层 来自 北京
本帖最后由 fantasy93404 于 2011-12-9 10:51 编辑
shalldy 发表于 2011-12-7 19:01
function KE = MechanicalRecEleStiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4)
E=1;NU=0.3;h=1;x1=0;y1=0;x2 ...


很给力啊,感谢
你给了我一个新思路,原来matlab不只是能进行矩阵计算,编程这么给力!
我按照这个思路编程了,求取点G(0,0)处的矩形单元刚度矩阵,结果不一样啊。。。。
function [ke]=elementstiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4)
E=1;NU=0.3;x1=0;y1=0;x2=1;y2=0;x3=1;y3=1;x4=0;y4=1;
A=[1 x1 y1 x1*y1 0 0 0 0
    0 0 0 0 1 x1 y1 x1*y1
    1 x2 y2 x2*y2 0 0 0 0
    0 0 0 0 1 x2 y2 x2*y2
     1 x3 y3 x3*y3 0 0 0 0
    0 0 0 0 1 x3 y3 x3*y3
     1 x4 y4 x4*y4 0 0 0 0
    0 0 0 0 1 x4 y4 x4*y4]
x=(x1+x2)/2
y=(y2+y3)/2
C=[0 1 0 y 0 0 0 0
    0 0 0 0 0 0 1 x
    0 0 1 x 0 1 0 y]
B=C/A
D=E*(1-NU)/((1+NU)*(1-2*NU))*[1 NU/(1-NU) 0
  NU/(1-NU) 1 0
  0 0 (1-2*NU)/(2*(1-NU))]
ke=B.'*D*B
你的结果是
ans =
  Columns 1 through 7
    0.4945    0.1786   -0.3022   -0.0137   -0.2473   -0.1786    0.0549
    0.1786    0.4945    0.0137    0.0549   -0.1786   -0.2473   -0.0137
   -0.3022    0.0137    0.4945   -0.1786    0.0549   -0.0137   -0.2473
   -0.0137    0.0549   -0.1786    0.4945    0.0137   -0.3022    0.1786
   -0.2473   -0.1786    0.0549    0.0137    0.4945    0.1786   -0.3022
   -0.1786   -0.2473   -0.0137   -0.3022    0.1786    0.4945    0.0137
    0.0549   -0.0137   -0.2473    0.1786   -0.3022    0.0137    0.4945
    0.0137   -0.3022    0.1786   -0.2473   -0.0137    0.0549   -0.1786
  Column 8
    0.0137
   -0.3022
    0.1786
   -0.2473
   -0.0137
    0.0549
   -0.1786
    0.4945
我的结果是
ans =
  Columns 1 through 7
    0.4327    0.2404   -0.2404    0.0481   -0.4327   -0.2404    0.2404
    0.2404    0.4327   -0.0481    0.2404   -0.2404   -0.4327    0.0481
   -0.2404   -0.0481    0.4327   -0.2404    0.2404    0.0481   -0.4327
    0.0481    0.2404   -0.2404    0.4327   -0.0481   -0.2404    0.2404
   -0.4327   -0.2404    0.2404   -0.0481    0.4327    0.2404   -0.2404
   -0.2404   -0.4327    0.0481   -0.2404    0.2404    0.4327   -0.0481
    0.2404    0.0481   -0.4327    0.2404   -0.2404   -0.0481    0.4327
   -0.0481   -0.2404    0.2404   -0.4327    0.0481    0.2404   -0.2404
  Column 8
   -0.0481
   -0.2404
    0.2404
   -0.4327
    0.0481
    0.2404
   -0.2404
    0.4327
回复 不支持

使用道具 举报

发表于 2011-12-8 16:23:44 | 显示全部楼层 来自 北京
yygao56 发表于 2011-12-8 09:31
挺给力的,给加点分

感谢版主的评分加了20仿真币
回复 不支持

使用道具 举报

发表于 2011-12-8 19:10:23 | 显示全部楼层 来自 湖北宜昌
fantasy93404 发表于 2011-12-8 16:22
很给力啊,感谢
你给了我一个新思路,原来matlab不只是能进行矩阵计算,编程这么给力!
我按照这个思路编 ...

显然是你的有问题,我的跟99行的结果是一样的。
不过我现在没时间看你的,有个代码要看,你自己先研究一下。
要不你加我QQ181802707讨论
回复 不支持

使用道具 举报

发表于 2011-12-9 10:28:48 | 显示全部楼层 来自 北京
shalldy 发表于 2011-12-8 19:10
显然是你的有问题,我的跟99行的结果是一样的。
不过我现在没时间看你的,有个代码要看,你自己先研究一 ...

很高兴有人一起来讨论问题,呵呵
QQ加上
回复 不支持

使用道具 举报

发表于 2011-12-9 10:38:48 | 显示全部楼层 来自 北京
正如楼主所说“这里不解释,自己可查有限元理论。个人认为这里的单元刚度求解不是很好”,这里,你的“结果和99行程序的结果是一样”反而说明了。。。。。,恕我直言,最起码你的结果楼主认为“求解不是很好”
还有,我的程序中单元刚度矩阵与单元所在的位置有关,不同的x,y处的单元具有不同的单刚值,所以,咱们的结果存在差异也是正常的,因为我只是计算了x=0,y=0的情形
回复 不支持

使用道具 举报

发表于 2011-12-17 05:21:08 | 显示全部楼层 来自 加拿大
top(nelx,nely,volfrac,penal,rmin) 这里的参数定义很重要,不同的参数会对拓扑结构有很大的影响。

一版penal取3,为了尽可能多得消除中值密度,因为在0和1之间的密度在数学角度上是不合理的。
rmin的取值是为了避免checkerboard而采用的sensitivity filtering. 另外拓扑优化还有个难点,就是mesh-denpendency,单元大小不同对结果也有很大影响。

附件中是我调了参数的对比。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×

点评

感谢eett的辛勤劳动,又是一位拓扑优化高手啊  发表于 2011-12-17 09:33

评分

1

查看全部评分

回复 不支持

使用道具 举报

发表于 2012-6-6 13:43:33 | 显示全部楼层 来自 北京
eett 发表于 2011-12-17 05:21
top(nelx,nely,volfrac,penal,rmin) 这里的参数定义很重要,不同的参数会对拓扑结构有很大的影响。

一版pe ...

真不错,左孔天的博士论文里好像也提到这点了
回复 不支持

使用道具 举报

发表于 2012-12-6 14:49:49 | 显示全部楼层 来自 广东广州
楼主还有研究这个不?毕设做这个想问些问题。
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-29 20:45 , Processed in 0.065188 second(s), 13 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表