结构分析的有限元法与MATLAB程序设计 里面k=stiffnessmat...

结构分析的有限元法与MATLAB程序设计 里面k=stiffnessmat...,第1张

限元

有限元法(FEA,Finite Element Analysis)的基本概念是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成,对每一单元假定一个合适的(较简单的)近似解,然后推导求解这个域总的满足条件(如结构的平衡条件),从而得到问题的解。这个解不是准确解,而是近似解,因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解,而有限元不仅计算精度高,而且能适应各种复杂形状,因而成为行之有效的工程分析手段。

英文:Finite Element 有限单元法是随着电子计算机的发展而迅速发展起来的一种现代计算方法。它是50年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法,随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 有限元法分析计算的思路和做法可归纳如下:

编辑本段1) 物体离散化

将某个工程结构离散为由各种单元组成的计算模型,这一步称作单元剖分。离散后单元与单元之间利用单元的节点相互连接起来;单元节点的设置、性质、数目等应视问题的性质,描述变形形态的需要和计算进度而定(一般情况单元划分越细则描述变形情况越精确,即越接近实际变形,但计算量越大)。所以有限元中分析的结构已不是原有的物体或结构物,而是同新材料的由众多单元以一定方式连接成的离散物体。这样,用有限元分析计算所获得的结果只是近似的。如果划分单元数目非常多而又合理,则所获得的结果就与实际情况相符合。

编辑本段2) 单元特性分析

A、 选择位移模式 在有限单元法中,选择节点位移作为基本未知量时称为位移法;选择节点力作为基本未知量时称为力法;取一部分节点力和一部分节点位移作为基本未知量时称为混合法。位移法易于实现计算自动化,所以,在有限单元法中位移法应用范围最广。 当采用位移法时,物体或结构物离散化之后,就可把单元总的一些物理量如位移,应变和应力等由节点位移来表示。这时可以对单元中位移的分布采用一些能逼近原函数的近似函数予以描述。通常,有限元法我们就将位移表示为坐标变量的简单函数。这种函数称为位移模式或位移函数。 B、 分析单元的力学性质 根据单元的材料性质、形状、尺寸、节点数目、位置及其含义等,找出单元节点力和节点位移的关系式,这是单元分析中的关键一步。此时需要应用d性力学中的几何方程和物理方程来建立力和位移的方程式,从而导出单元刚度矩阵,这是有限元法的基本步骤之一。 C、 计算等效节点力 物体离散化后,假定力是通过节点从一个单元传递到另一个单元。但是,对于实际的连续体,力是从单元的公共边传递到另一个单元中去的。因而,这种作用在单元边界上的表面力、体积力和集中力都需要等效的移到节点上去,也就是用等效的节点力来代替所有作用在单元上的力。

编辑本段3) 单元组集

利用结构力的平衡条件和边界条件把各个单元按原来的结构重新连接起来,形成整体的有限元方程 (1-1) 式中,K是整体结构的刚度矩阵;q是节点位移列阵;f是载荷列阵。

编辑本段4) 求解未知节点位移

解有限元方程式(1-1)得出位移。这里,可以根据方程组的具体特点来选择合适的计算方法。 通过上述分析,可以看出,有限单元法的基本思想是"一分一合",分是为了就进行单元分析,合则为了对整体结构进行综合分析。 有限元的发展概况 1943年 courant在论文中取定义在三角形域上分片连续函数,利用最小势能原理研究St.Venant的扭转问题。 1960年 clough的平......

当然有了。不过现在的有限元代码多了去了,什么语言都用。

我见过用Fortran 的, c++的,C#的,Java的,一句话,基本上任何一种高级语言都可以写有限元代码。

Fortran是科学计算比较早期的语言,运行效率较高。很多商业软件核心是用它写的。

clc

clear

tic

% n 是行与列划分的格子数,对整个[0,1]*[0,1]有n^2个划分,can为方程中的参数k

%这里我们用1,2,3按逆时针来表示一个三角形的各个顶点

% s是一个n^2*10的关联矩阵,s(i,1)表示第i个三角形,s(i,2),s(i,3),s(i,4)分别表示第i个三角形的1,2,3所对应的顶点

% s(i,5),s(i,6);s(i,7),s(i,8);s(i,9),s(i,10)分别表示顶点1,2,3所代表的坐标

%生成关联矩阵s

%A是总刚矩阵

%声明符号变量x y

n=20

can=20

s=zeros(2*n^2,10)

h=1/n

st=1/(2*n^2)

A=zeros((n+1)^2,(n+1)^2)

syms x y

for  k=1:1:2*n^2

        s(k,1)=k

        q=fix(k/(2*n))

        r=mod(k,(2*n))

        if (r~=0)

            r=r

        else r=2*nq=q-1

        end

        if (r<=n)

            s(k,2)=q*(n+1)+r

            s(k,3)=q*(n+1)+r+1

            s(k,4)=(q+1)*(n+1)+r+1

            s(k,5)=(r-1)*h

            s(k,6)=q*h

            s(k,7)=r*h

            %%

            % 

            %   for x = 1:10

            %       disp(x)

            %   end

            % 

            s(k,8)=q*h

            s(k,9)=r*h

            s(k,10)=(q+1)*h

        else

            s(k,2)=q*(n+1)+r-n

            s(k,3)=(q+1)*(n+1)+r-n+1

            s(k,4)=(q+1)*(n+1)+r-n

            s(k,5)=(r-n-1)*h

            s(k,6)=q*h

            s(k,7)=(r-n)*h

            s(k,8)=(q+1)*h

            s(k,9)=(r-n-1)*h

            s(k,10)=(q+1)*h

        end

end

%下面生成基函数L(i)表示第i个点顶点的基函数

%生成单刚矩阵d

%生成单刚矩阵并将其加入总纲矩阵

d=zeros(3,3)

B=zeros((n+1)^2,1)

b=zeros(3,1)

%生成A的总刚

for   k=1:1:2*n^2

     L(1)=(1/(2*st))*((s(k,7)*s(k,10)-s(k,9)*s(k,8))+(s(k,8)-s(k,10))*x+(s(k,9)-s(k,7))*y)

     L(2)=(1/(2*st))*((s(k,9)*s(k,6)-s(k,5)*s(k,10))+(s(k,10)-s(k,6))*x+(s(k,5)-s(k,9))*y)

     L(3)=(1/(2*st))*((s(k,5)*s(k,8)-s(k,7)*s(k,6))+(s(k,6)-s(k,8))*x+(s(k,7)-s(k,5))*y)

     for i=1:1:3

        for j=i:3

            d(i,j)=int(int(((((diff(L(i),x))*(diff(L(j),x)))+((diff(L(i),y))*(diff(L(j),y))))-((can^2)*L(i)*L(j))),x,0,1),y,0,1)

            d(j,i)=d(i,j)

        end

     end   

     for i=1:1:3

         for j=1:1:3

             A(s(k,(i+1)),s(k,(j+1)))=A(s(k,(i+1)),s(k,(j+1)))+d(i,j)

         end

     end

     for i=1:1:3

         b(i)=int(int((L(i)),x,0,1),y,0,1)

         B(s(k,(i+1)),1)=B(s(k,(i+1)),1)+b(i)

     end

      

end

%下面根据边界条件来求解有限元方程组Mx=B,齐次边界条件约掉了很多项

M=zeros((n+1)^2,n^2)

j=n^2

for i=(n^2+n):-1:1

    if ((mod(i,(n+1)))~=1)

        M(:,j)=A(:,i)

        j=j-1

    else continue

    end

end

%preanswer是未知点的值是(n+1)^2*(n^2)的

preanswer=M\B

%得到所有节点的值

answer=zeros((n+1)^2,1)

j=1

for i=1:1:(n^2+n)

    if ((mod(i,(n+1)))~=1)

        answer(i)=preanswer(j)

        j=j+1

    else answer(i)=0

    end

end

%%

%   for x = 1:10

%   for x = 1:10

%   for x = 1:10

%       disp(x)

%   end

%       disp(x)

%   end

%       disp(x)

%   end

%生成求解后的图

Z=zeros((n+1),(n+1))

for i=1:1:(n+1)^2

      s=fix(i/(n+1))+1

      r=mod(i,(n+1))

      if(r==0)

          r=n+1

          s=s-1

      else

      end

      Z(r,s)=answer(i)

end

[X,Y]=meshgrid(1:-h:0,0:h:1)

surf(X,Y,Z)

 

 

 

toc

t=toc


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/11048388.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-13
下一篇 2023-05-13

发表评论

登录后才能评论

评论列表(0条)

保存