有限元方法[Matlab]-笔记

合集 - 数值计算(17)

1.【数学公式】mathtype和word2016集成2023-04-19 2.【MathType】word2016数学公式编号2023-04-19 3.【Matlab】基于KDtree的最近邻搜索和范围搜索2023-07-30 4.【Matlab】判断点和多面体位置关系的两种方法实现2023-07-29 5.《空间三角面片对相交判断算法》的matlab实现_ 0.2微秒2023-11-11 6.Mpmath库-学习笔记2024-07-02

7.有限元方法[Matlab]-笔记2024-04-13

<-- 访问笔记代码仓库 -->

<< MATLAB Codes for Finite Element Analysis - Solids and Structures (Ferreira) >> 笔记

第一章:matlab基础
第二章: 离散系统
- 笔记、例题
- Matlab代码
  - problem1.m
  - outputDisplacementsReactions.m
第三章:1D 杆单元
第四章 2D桁架分析
第五章 3D 桁架分析
第六章 2D 欧拉-伯努利梁分析
第7章 2d Frame单元静态分析
- 单元刚度矩阵推导
- 例题
  - Matlab 代码
第八章 3d Frame单元静态分析
第九章 analysis of grids
第十章 Timoshenko beams 静态分析
- 10.Timoshenko梁的静力分析-Formulation

第一章:matlab基础

略

第二章: 离散系统

笔记、例题

Matlab代码

problem1.m

 % MATLAB codes for Finite Element Analysis
 
%  Problem 1: 3 springs problem
 
% clear memory
clear ;
 
% ? define vars
P=10; % N
ki=[1;1;1;1];
 
% elementNodes: connections at elements
elementNodes=[1 2;
            2 3;
            2 4];
 
% numberElements: number of Elements
numberElements=size(elementNodes,1);
 
% numberNodes: number of nodes
numberNodes=4;
 
% for structure:
% displacements: displacement vector N x 1
% force : force vector N x 1
% stiffness: stiffness matrix N x N
displacements=zeros(numberNodes,1); %#ok<PREALL>
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes);
 
% applied load at node 2
force(2)=P;
% computation of the system stiffness matrix
for e=1:numberElements
    % elementDof: element degrees of freedom (Dof)
    elementDof=elementNodes(e,:) ;
    stiffness(elementDof,elementDof)=...
    stiffness(elementDof,elementDof)+ki(e)*[1 -1;-1 1];
end
% boundary conditions and solution
% prescribed dofs
prescribedDof=[1;3;4];
% free Dof : activeDof
% setdiff() 函数-->求两个数组的差集
% C = setdiff(A,B) 返回 A 中存在但 B 中不存在的数据,不包含重复项.C 是有序的.
activeDof=setdiff((1:numberNodes)',prescribedDof);
 
% solution
% !KU=F --> U=K\F
displacements=stiffness(activeDof,activeDof)\force(activeDof);
 
% positioning all displacements
displacements1=zeros(numberNodes,1);
displacements1(activeDof)=displacements;
 
% output displacements/reactions
outputDisplacementsReactions(displacements1,stiffness,numberNodes,prescribedDof)

outputDisplacementsReactions.m

 function outputDisplacementsReactions(displacements,stiffness,GDof,prescribedDof)
    % output of displacements and reactions in
    % tabular form
    % GDof: total number of degrees of freedom of the problem
    % displacements
    disp('Displacements')
    %displacements=displacements1;
    jj=1:GDof; % format
    [jj' displacements]
    % reactions
    F=stiffness*displacements;
    reactions=F(prescribedDof);
    disp('reactions')
    [prescribedDof reactions]
end

第三章:1D 杆单元

笔记

例题1

matlab 代码

problem2.m

 % antonio ferreira 2008
% clear memory
clear all
% E; modulus of elasticity
% A: area of cross section
% L: length of bar
E = 30e6;A=1;EA=E*A; L = 90;
% generation of coordinates and connectivities
% numberElements: number of elements
numberElements=3;
% generation equal spaced coordinates
nodeCoordinates=linspace(0,L,numberElements+1);
xx=nodeCoordinates;
% numberNodes: number of nodes
numberNodes=size(nodeCoordinates,2);
% elementNodes: connections at elements
ii=1:numberElements;
elementNodes(:,1)=ii;
elementNodes(:,2)=ii+1;
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
displacements=zeros(numberNodes,1); %#ok<PREALL>
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
 
% applied load at node 2
force(2)=3000.0;
% computation of the system stiffness matrix
for e=1:numberElements;
    % elementDof: element degrees of freedom (Dof)
    elementDof=elementNodes(e,:) ;
    nn=length(elementDof);
    length_element=nodeCoordinates(elementDof(2))...
    -nodeCoordinates(elementDof(1));
    detJacobian=length_element/2;invJacobian=1/detJacobian;
 
    % central Gauss point (xi=0, weight W=2)
    [shape,naturalDerivatives]=shapeFunctionL2(0.0);
    Xderivatives=naturalDerivatives*invJacobian;
    % B matrix
    B=zeros(1,nn); B(1:nn) = Xderivatives(:);
    stiffness(elementDof,elementDof)=stiffness(elementDof,elementDof)+B'*B*2*detJacobian*EA;
end
 
% boundary conditions and solution
% prescribed dofs
fixedDof=find(xx==min(nodeCoordinates(:)) | xx==max(nodeCoordinates(:)))';
prescribedDof=[fixedDof]
% free Dof : activeDof
activeDof=setdiff([1:numberNodes]',[prescribedDof]);
 
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
 
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness,numberNodes,prescribedDof)
 
disp(stiffness)
disp(force)
disp(displacements)

solution.m

 function displacements=solution(GDof,prescribedDof,stiffness,force)
    % function to find solution in terms of global displacements
    activeDof=setdiff([1:GDof]', ...
    [prescribedDof]);
    U=stiffness(activeDof,activeDof)\force(activeDof);
    displacements=zeros(GDof,1);
    displacements(activeDof)=U;
end

shapeFunctionL2.m

 function [shape,naturalDerivatives]=shapeFunctionL2(xi)
    % shape function and derivatives for L2 elements
    % shape : Shape functions
    % naturalDerivatives: derivatives w.r.t. xi
    % xi: natural coordinates (-1 ... +1)
    shape=([1-xi,1+xi]/2)';
    naturalDerivatives=[-1;1]/2;
end

例题2

matlab 代码

problem3Structure.m

 % MATLAB codes for Finite Element Analysis
 
% antonio ferreira 2008
% clear memory
clear all
% p1 : structure
p1=struct();
 
% E; modulus of elasticity
% A: area of cross section
% L: length of bar
% k: spring stiffness
E=70000;
A=200;
k=2000;
 
% generation of coordinates and connectivities
% numberElements: number of elements
p1.numberElements=3;
p1.numberNodes=4;
p1.elementNodes=[1 2; 2 3; 3 4];
p1.nodeCoordinates=[0 2000 4000 4000];
p1.xx=p1.nodeCoordinates;
% GDof: total degrees of freedom
p1.GDof=p1.numberNodes;
 
% for structure:
    % displacements: displacement vector
    % force : force vector
    % stiffness: stiffness matrix
p1.displacements=zeros(p1.GDof,1);
p1.force=zeros(p1.GDof,1);
p1.stiffness=zeros(p1.GDof);
 
% applied load at node 2
p1.force(2)=8000.0;
% computation of the system stiffness matrix
for e=1:p1.numberElements;
    % elementDof: element degrees of freedom (Dof)
    elementDof=p1.elementNodes(e,:) ;
    L=p1.nodeCoordinates(elementDof(2))-p1.nodeCoordinates(elementDof(1));
    if e<3
        ea(e)=E*A/L;
    else
        ea(e)=k;
    end
    p1.stiffness(elementDof,elementDof)=p1.stiffness(elementDof,elementDof)+ea(e)*[1 -1;-1 1];
end
 
% prescribed dofs
p1.prescribedDof=[1;4];
 
% solution
p1.displacements=solutionStructure(p1);
 
% output displacements/reactions
outputDisplacementsReactionsStructure(p1)
 
% Displacements
 
% ans =
 
%     1.0000         0
%     2.0000    3.3333
%     3.0000         0
%     4.0000         0
 
 
% reactions
 
% ans =
 
%     1.0000   -3.3333
%     3.0000   -3.3333
%     4.0000   -3.3333

outputDisplacementsReactionsStructure.m

 function outputDisplacementsReactionsStructure(p)
    % output of displacements and reactions in
    % tabular form
    % GDof: total number of degrees of freedom of
    % the problem
    % displacements
    disp('Displacements')
    jj=1:p.GDof; 
    [jj' p.displacements]
    % reactions
    F=p.stiffness*p.displacements;
    reactions=F(p.prescribedDof);
    disp('reactions')
    [p.prescribedDof reactions]

shapeFunctionL2Structure.m

 function shapeL2=shapeFunctionL2Structure(xi)
    % shape function and derivatives for L2 elements
    % shape : Shape functions
    % naturalDerivatives: derivatives w.r.t. xi
    % xi: natural coordinates (-1 ... +1)
    shapeL2=struct()
    shapeL2.shape=([1-xi,1+xi]/2)';
    shapeL2.naturalDerivatives=[-1;1]/2;
 
end % end function shapeFunctionL2

solutionStructure.m

 function displacements = solutionStructure(p)
    % function to find solution in terms of global displacements
    activeDof=setdiff([1:p.GDof]', [p.prescribedDof]);
    U=p.stiffness(activeDof,activeDof)\p.force(activeDof);
    displacements=zeros(p.GDof,1);
    displacements(activeDof)=U;
end

第四章 2D桁架分析

笔记:单元刚度矩阵推导

例题1

模型参数:

杨氏模量E=30e6 Pa ,泊松比 \(\mu\) =0.3
A=2 m^2
L=10 m
P=10000 N

手算过程

Abaqus 模拟结果

inp文件在码云repo中

U and RF result :

  
********************************************************************************
Probe Values Report, written on Tue Apr 16 21:43:02 2024
 
Source
-------
 
   ODB: C:/Users/Administrator/Documents/learnFEM/Matlb Fea Codes repo/fea-codes-of-matlab/chapter04/problem4/abaqus/Job-2.odb
   Step: Step-1
   Frame: Increment      1: Step Time =    1.000
   Variable: U, U1
 
Probe values reported at nodes
 
  Part Instance     Node ID                   Orig. Coords                                     Def. Coords                 
                                    X               Y               Z               X               Y               Z      
-------------------------------------------------------------------------------------------------------------------------
    PART-1-1           1              0.              0.              0.     345.178E-06    -1.32149E-03              0.
    PART-1-1           2              0.             10.              0.              0.             10.              0.
    PART-1-1           3             10.             10.              0.             10.             10.              0.
    PART-1-1           4             10.              0.              0.             10.              0.              0.
 
 
  Part Instance     Node ID  Attached elem      U, U1        U, U2           RF, RF1         RF, RF2  
-------------------------------------------------------------------------------------------------------------- 
    PART-1-1           1         3          345.178E-06    -1.32149E-03              0.              0.        
    PART-1-1           2         1                   0.    -7.92893E-33              0.     7.92893E+03        
    PART-1-1           3         2         -2.07107E-33    -2.07107E-33     2.07107E+03     2.07107E+03        
    PART-1-1           4         3          2.07107E-33              0.    -2.07107E+03             -0.

可以看出,手算和abaqus的U结果完全一致.

matlab 代码

problem4.m

 % ! problem4.m --> based on <<MATLAB codes for Finite Element Analysis>>
 
clear ;
clc
format default
p=struct();
 
% define model parameters
 
p.E = 30e6;
p.A = 2;
p.L = 10;
p.P = -10000;
 
% define elems and nodes
 
p.nodes=[0 0 ;  % node's coordinates (x,y)
         0 p.L ;
         p.L p.L ;
         p.L 0 ];
p.elems=[1 2 ; 
        1 3 ; 
        1 4];
 
p.node_num=size(p.nodes,1);
p.elem_num=size(p.elems,1);
p.node_Coord_x=p.nodes(:,1);
p.node_Coord_y=p.nodes(:,2);
 
% define problem's dimension : 1D/2D/3D;
%        = each node's dof  
p.problem_dimension=size(p.nodes,2);
 
% global degree of freedom number
p.global_dof_num=2*p.node_num;
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);
 
% define boundary conditions
% fixed dof 
p.fix_dof=[3 4 5 6 7 8]';
% load dof and amplitude
p.load_dof=2;
 
% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);
 
% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=2*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);
 
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 2d truss node has 2 node,each node have 2 dof,
    % i-th node --> U_(2*i-1) and U_(2*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2 ...
             connectivity(2)*2-1 connectivity(2)*2];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
    l=(node2_x-node1_x)/elem_length;
    m=(node2_y-node1_y)/elem_length;
 
    k_e=(p.E*p.A/elem_length)*...
        [l*l l*m -l*l -l*m;
         l*m m*m -l*m -m*m;
        -l*l -l*m l*l l*m;
        -l*m -m*m l*m m*m];
    p.elemStiffs(:,:,i)=k_e;
 
    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
 
end
 
% apply boundary condition 
p.displacements(p.fix_dof)=0;
p.node_forces(p.load_dof)=p.P;
 
% solution KU=F
p = solutionStruct(p);
 
% show compute result
disp('displacements of dofs :')
disp([(1:p.global_dof_num)', p.displacements])
 
disp('forces of dofs :')
disp([(1:p.global_dof_num)',p.node_forces])
 
figure
axis equal
 
% plot undeformed mesh
draw_2Dtruss_mesh(p.nodes,p)
 
% plot deformed mesh
deformed_node_coordinate=cal_deformed_nodes(p,500.0);
draw_2Dtruss_mesh(deformed_node_coordinate,p)
 
title('位移变形')
grid on
 
 
function new_p=solutionStruct(p)
    % free dof in structure
    free_dof=setdiff((1:p.global_dof_num)',p.fix_dof);
    
    % solute the equation :KU=F --> U=K\F
    U=p.global_stiffness_matrix(free_dof,free_dof)\p.node_forces(free_dof);
    
    p.displacements(free_dof)=U;
    % displacements=zeros(p.global_dof_num,1);
    
    p.node_forces=p.global_stiffness_matrix*p.displacements;
    
    new_p=p;
end
 
function deformed_node_coordinate=cal_deformed_nodes(p,scale)
    % initial 
    deformed_node_coordinate=p.nodes;
    if p.problem_dimension==2
        index1=1:2:p.global_dof_num;
        index2=2:2:p.global_dof_num;
        deformed_node_coordinate(:,1)=deformed_node_coordinate(:,1)+p.displacements(index1)*scale;
        deformed_node_coordinate(:,2)=deformed_node_coordinate(:,2)+p.displacements(index2)*scale;
    elseif p.problem_dimension==3
        index1=1:3:p.global_dof_num;
        index2=2:3:p.global_dof_num;
        index3=3:3:p.global_dof_num;
        deformed_node_coordinate(:,1)=deformed_node_coordinate(:,1)+p.displacements(index1);
        deformed_node_coordinate(:,2)=deformed_node_coordinate(:,2)+p.displacements(index2);
        deformed_node_coordinate(:,3)=deformed_node_coordinate(:,3)+p.displacements(index3);
    end
end
 
function draw_2Dtruss_mesh(nodes,p)
    xy=zeros(size(p.elems,2)*p.elem_num,p.problem_dimension);
    cur1=1;
    cur2=0;
    for i =1:p.elem_num
        
        connectivity=p.elems(i,:);
        attach_node_num=size(connectivity,2);
        
        cur1=cur2+1;
        cur2=cur2+attach_node_num;
        
        coord=zeros(attach_node_num,p.problem_dimension);
        for j=1:size(connectivity,2)
            coord(j,:)=nodes(connectivity(j),:);
        end
    
        xy(cur1:cur2,:)=coord;
    end
    line(xy(:,1),xy(:,2),"Marker","o","LineWidth",1,"Color",rand(1,3))
end

三种方式,计算结果一致

例题2

模型参数:

杨氏模量E = 70 GPa,泊松比 \(\mu\) = 0.3
桁架横截面积 A = \(3e^{-4} m^2\)
L = 3 m

Abaqus 结果

待补充...

matlab 计算

基于problem4.m 代码,作如下修改即可:

  
% define model parameters
% unit:SI(m)
p.E = 70e9;
p.A = 3e-4;
p.L = 3;
p.P = 1000*[-50;-100;-50];
disp('单位制:SI(m)')
% define elems and nodes
 
p.nodes=[0 0 ;  % node's coordinates (x,y)
        0 p.L ;
        p.L 0 ;
        p.L p.L ;
        2*p.L 0 ;
        2*p.L p.L];
p.elems=[1 2 ; 
        1 3; 
        2 3;
        2 4;
        1 4;
        3 4;
        3 6;
        4 5;
        4 6;
        3 5;
        5 6];
 
......
 
% define boundary conditions
% fixed dof 
p.fix_dof=[1 2 10]';
% load dof and amplitude
p.load_dof=[4 8 12]';
 
......

 单位制:SI(m)
displacements of dofs :
    1.0000         0
    2.0000         0
    3.0000    0.0071
    4.0000   -0.0090
    5.0000    0.0052
    6.0000   -0.0163
    7.0000    0.0052
    8.0000   -0.0201
    9.0000    0.0105
   10.0000         0
   11.0000    0.0034
   12.0000   -0.0090
 
forces of dofs :
   1.0e+05 *
 
    0.0000    0.0000
    0.0000    1.0000
    0.0000    0.0000
    0.0000   -0.5000
    0.0001   -0.0000
    0.0001   -0.0000
    0.0001   -0.0000
    0.0001   -1.0000
    0.0001    0.0000
    0.0001    1.0000
    0.0001         0
    0.0001   -0.5000

第五章 3D 桁架分析

笔记 - 3d truss 单元刚度矩阵

例题 1

matlab 计算

problem7.m 代码

 % ! problem7.m --> based on <<MATLAB codes for Finite Element Analysis>>
 
clear ;
clc
format default
p=struct();
 
% define model parameters
% unit:SI(m)
p.E = 1.2e6;
p.A = [0.302 0.729 0.187]';
p.P = -1000;
disp('单位制:SI(m)')
% define elems and nodes
 
p.nodes=[72 0 0 ;  % node's coordinates (x,y,z)
        0 36 0;
        0 36 72;
        0 0 -48];
 
p.elems=[1 2 ; 
        1 3;
        1 4];
 
p.node_num=size(p.nodes,1);
p.elem_num=size(p.elems,1);
p.node_Coord_x=p.nodes(:,1);
p.node_Coord_y=p.nodes(:,2);
p.node_Coord_z=p.nodes(:,3);
 
% define problem's dimension : 1D/2D/3D;
%        = each node's dof  
p.problem_dimension=size(p.nodes,2);
 
% global degree of freedom number
p.global_dof_num=p.problem_dimension*p.node_num;
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);
 
% define boundary conditions
% fixed dof 
p.fix_dof=[2 4 5 6 7 8 9 10 11 12]';
% load dof and amplitude
p.load_dof=3;
 
% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);
 
% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=p.problem_dimension*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);
 
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 3 dof,
    % i-th node --> U_(3*i-2) / U_(3*i-1) / U_(3*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*3-2 connectivity(1)*3-1 connectivity(1)*3 ...
        connectivity(2)*3-2 connectivity(2)*3-1 connectivity(2)*3];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node1_z=p.node_Coord_z(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    node2_z=p.node_Coord_z(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2+(node1_z-node2_z).^2);
    Cx=(node2_x-node1_x)/elem_length;
    Cy=(node2_y-node1_y)/elem_length;
    Cz=(node2_z-node1_z)/elem_length;
 
    T=[Cx*Cx Cx*Cy Cx*Cz;
        Cx*Cy Cy*Cy Cy*Cz;
        Cx*Cz Cy*Cz Cz*Cz];
 
    k_e=(p.E*p.A(i)/elem_length)*[T -T;-T T];
    p.elemStiffs(:,:,i)=k_e;
 
    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
 
end
 
% apply boundary condition 
p.displacements(p.fix_dof)=0;
p.node_forces(p.load_dof)=p.P;
 
% solution KU=F
p = solutionStruct(p);
 
% calculate stress
p.stress=zeros(p.elem_num,1);
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 3 dof,
    % i-th node --> U_(3*i-2) / U_(3*i-1) / U_(3*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*3-2 connectivity(1)*3-1 connectivity(1)*3 ...
        connectivity(2)*3-2 connectivity(2)*3-1 connectivity(2)*3];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node1_z=p.node_Coord_z(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    node2_z=p.node_Coord_z(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2+(node1_z-node2_z).^2);
    Cx=(node2_x-node1_x)/elem_length;
    Cy=(node2_y-node1_y)/elem_length;
    Cz=(node2_z-node1_z)/elem_length;
 
    p.stress(i)=(p.E/elem_length)*[-Cx -Cy -Cz Cx Cy Cz]*p.displacements(elem_dof);
end
 
% show compute result
disp('displacements of dofs :')
disp([(1:p.global_dof_num)', p.displacements])
 
disp('forces of dofs :')
disp([(1:p.global_dof_num)',p.node_forces])
 
disp('3d truss elem stress :')
disp([(1:p.elem_num)',p.stress])
 
figure
axis equal
 
% plot undeformed mesh
draw_3Dtruss_mesh(p.nodes,p)
 
% plot deformed mesh
deformed_node_coordinate=cal_deformed_nodes(p,50.0);
draw_3Dtruss_mesh(deformed_node_coordinate,p)
 
title('位移变形')
grid on
 
 
function draw_3Dtruss_mesh(nodes,p)
    xy=zeros(size(p.elems,2)*p.elem_num,p.problem_dimension);
    cur1=1;
    cur2=0;
    for i =1:p.elem_num
        
        connectivity=p.elems(i,:);
        attach_node_num=size(connectivity,2);
        
        cur1=cur2+1;
        cur2=cur2+attach_node_num;
        
        coord=zeros(attach_node_num,p.problem_dimension);
        for j=1:size(connectivity,2)
            coord(j,:)=nodes(connectivity(j),:);
        end
    
        xy(cur1:cur2,:)=coord;
    end
    line(xy(:,1),xy(:,2),xy(:,3),"Marker","o","LineWidth",1,"Color",rand(1,3))
end

计算结果

 单位制:SI(m)
displacements of dofs :
    1.0000   -0.0711
    2.0000         0
    3.0000   -0.2662
    4.0000         0
    5.0000         0
    6.0000         0
    7.0000         0
    8.0000         0
    9.0000         0
   10.0000         0
   11.0000         0
   12.0000         0
 
 
forces of dofs :
   1.0e+03 *
 
    0.0010         0
    0.0020   -0.2232
    0.0030   -1.0000
    0.0040    0.2561
    0.0050   -0.1281
    0.0060         0
    0.0070   -0.7024
    0.0080    0.3512
    0.0090    0.7024
    0.0100    0.4463
    0.0110         0
    0.0120    0.2976
 
 
3d truss elem stress :
   1.0e+03 *
 
    0.0010   -0.9482
    0.0020    1.4454
    0.0030   -2.8685

例题2

模型参数:

E = 210000 Mpa,\(\mu\)=0.3
A = 100 mm^2
P =-10000 N

Abaqus 结果

waiting...

Matlab 计算

problem8.m

基于problem7.m,修改参数,即可.

  
% define model parameters
% unit:SI(mm)
p.E = 210000;
p.A = 100;
p.P = -10000;
disp('单位制:SI(mm)')
% define elems and nodes
 
p.nodes=[4 4 3 ;  % node's coordinates (x,y,z)
        0 4 0;
        0 4 6;
        4 0 3;
        8 -1 1];
p.elems=[1 2 ; 
        1 3;
        1 4;
        1 5];
 
......
 
% define boundary conditions
% fixed dof 
p.fix_dof=[4:15]';
% load dof and amplitude
p.load_dof=2;
 
......

计算结果

 单位制:SI(mm)
displacements of dofs :
    1.0000   -0.0003
    2.0000   -0.0015
    3.0000    0.0003
 
forces of dofs :
   1.0e+04 *
 
    0.0004    0.0271
    0.0005         0
    0.0006    0.0203
    0.0007    0.1355
    0.0008         0
    0.0009   -0.1016
    0.0010         0
    0.0011    0.7968
    0.0012         0
    0.0013   -0.1626
    0.0014    0.2032
    0.0015    0.0813
 
3d truss elem stress :
    1.0000   -3.3865
    2.0000  -16.9326
    3.0000  -79.6809
    4.0000  -27.2610
 
deformed nodes:
    3.8791    3.3929    3.1075
         0    4.0000         0
         0    4.0000    6.0000
    4.0000         0    3.0000
    8.0000   -1.0000    1.0000

第六章 2D 欧拉-伯努利梁分析

单元刚度矩阵推导

例题 1

consider a simply-supported and clamped Bernoulli beam in bending, under uniform load

Abaqus 计算

Matlab 计算

 % ! problem9.m --> based on <<MATLAB codes for Finite Element Analysis>>
 
clear ;
clc
format default
p=struct();
Case=1;
 
% define model parameters
% unit:SI(mm)
p.EI = 1;
p.L = 1;
% * distribute force
p.P = -1;
% define elems num of bornoulli beam
p.elem_num=80;
% define problem's dimension : 1D/2D/3D;
%        = each node's dof  
p.problem_dimension=2;
 
% ? -------------------------------
 
% node's coordinates (x)
p.nodes=linspace(0,p.L,p.elem_num+1)';
p.elems=zeros(p.elem_num,2);
for i=1:p.elem_num
    p.elems(i,1)=i;
    p.elems(i,2)=i+1;
end
p.node_num=size(p.nodes,1);
p.node_Coord_x=p.nodes(:,1);
 
% global degree of freedom number
p.global_dof_num=2*p.node_num;
 
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);
 
% define boundary conditions
if Case==1
    % Case 01: clamped at x=0
    % fixed dof 
    p.fix_dof=[1 2]';
end
if Case==2
    % Case 02: clamped-clamped
    % fixed dof 
    p.fix_dof=[1 p.elem_num*2+1 2 p.elem_num*2+2]';
end
if Case==3
    % Case 03: simply supported-simply supported 
    p.fix_dof=[1 2*p.elem_num+1]';
end
 
% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);
 
% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=p.problem_dimension*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);
 
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
    % i-th node --> w_(2*i-1) / dw_(2*i-2) /  dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2  ...
        connectivity(2)*2-1 connectivity(2)*2 ];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2)); 
    
    elem_length=node2_x-node1_x;
    a=0.5*elem_length;
    k_e=(p.EI/(a^3))*0.5*[3 3*a -3 3*a;
                        3*a 4*a*a -3*a 2*a*a;
                        -3 -3*a 3 -3*a;
                        3*a 2*a*a -3*a 4*a*a];
    p.elemStiffs(:,:,i)=k_e;
    p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.P*a/3)*[3 a 3 -a])';
 
    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
end
 
% apply boundary condition 
p.displacements(p.fix_dof)=0;
 
% solution KU=F
p = solutionStruct_bornoulli_beam(p);
 
plot(p.node_Coord_x,p.displacements(1:2:2*p.node_num),'.')
minDispalce=min(p.displacements(1:2:2*p.node_num));
text(0.35,0.5*minDispalce,"挠度最小值:"+minDispalce)
 
if Case==1
    title("Case 01: clamped at x=0")
end
if Case==2
    title("Case 02: clamped-clamped")
end
if Case==3
    title("Case 03: simply supported-simply supported") 
end
 
 
function new_p=solutionStruct_bornoulli_beam(p)
    % free dof in structure
    free_dof=setdiff((1:p.global_dof_num)',p.fix_dof);
    
    % solute the equation :KU=F --> U=K\F
    U=p.global_stiffness_matrix(free_dof,free_dof)\p.node_forces(free_dof);
    
    p.displacements(free_dof)=U;
    % displacements=zeros(p.global_dof_num,1);
    
    p.node_forces_fu=p.global_stiffness_matrix*p.displacements;
    
    new_p=p;
end

例题2

matlab 计算

problem10.m

 % ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>
 
clear;
clc;
format compact;
 
p=struct();
 
% define model paraeters
p.E=1000000;
p.L=10;
t=p.L/1000;
p.Iz=(1*t^3)/12;
p.P=-1000;
% define elems of bornoulli beam
p.elem_num=20;  
% dof of each nornoulli beam nodes
p.bornoulli_node_dof_num=2;
p.spring_k=100;
 
% -------------------------------
 
% node's coordinate
p.nodes=linspace(0,p.L,p.elem_num+1)';
p.elems=zeros(p.elem_num,2);
for i=1:p.elem_num
    p.elems(i,1)=i;
    p.elems(i,2)=i+1;
end
p.node_num=size(p.nodes,1);
p.node_Coord_x=p.nodes(:,1);
 
% global degree of freedom number: 
%   ! bornoulli_beam node's dof + 1 spring node's dof
p.global_dof_num=2*p.node_num;
 
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num+1,1);
p.node_forces=zeros(p.global_dof_num+1,1);
 
% define boundary conditions
%   ! BC: clamped at x=0
p.fix_dof=[1,2,p.global_dof_num+1]';
 
% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num+1);
 
% compute all ElemStiffnessMatrix and assembly stiff matrix
beam_dof_num=p.bornoulli_node_dof_num*size(p.elems,2);
p.elemStiffs=zeros(beam_dof_num,beam_dof_num,p.elem_num);
 
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
    % i-th node --> w_(2*i-1) / dw_(2*i-2) /  dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2  ...
        connectivity(2)*2-1 connectivity(2)*2 ];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2)); 
    
    elem_length=node2_x-node1_x;
    a=0.5*elem_length;
    k_e=(p.E*p.Iz/(a^3))*0.5*[3 3*a -3 3*a;
                        3*a 4*a*a -3*a 2*a*a;
                        -3 -3*a 3 -3*a;
                        3*a 2*a*a -3*a 4*a*a];
    p.elemStiffs(:,:,i)=k_e;
    p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.P*a/3)*[3 a 3 -a])';
 
    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
end
 
% add spring nodes into global_stiffness_matrix
p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])=...
    p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])+p.spring_k*[1,-1;-1,1];
 
% apply boundary condition
p.node_forces(p.fix_dof)=0;
 
% solution KU=F
p=solutionStruct(p);
 
% ! display result
disp('Bornoulli beam result !')
 
disp('beam挠度:')
omega_index=(1:2:p.global_dof_num)';
disp([omega_index,p.displacements(omega_index)])
 
disp('beam转角:')
theta_index=(2:2:p.global_dof_num)';
disp([theta_index,p.displacements(theta_index)])
 
disp('beam支反力:')
disp([p.fix_dof,p.node_forces(p.fix_dof)])
 
plot(p.node_Coord_x,p.displacements(omega_index),'-')
title('Problem10.m')

计算结果

 Bornoulli beam result !
beam挠度:
   1.0e+05 *
    0.0000         0
    0.0000   -0.1722
    0.0001   -0.6300
    0.0001   -1.2909
    0.0001   -2.0800
    0.0001   -2.9297
    0.0001   -3.7800
    0.0001   -4.5785
    0.0002   -5.2801
    0.0002   -5.8473
    0.0002   -6.2501
    0.0002   -6.4661
    0.0003   -6.4802
    0.0003   -6.2849
    0.0003   -5.8802
    0.0003   -5.2737
    0.0003   -4.4803
    0.0003   -3.5225
    0.0004   -2.4303
    0.0004   -1.2413
    0.0004   -0.0004
 
beam转角:
   1.0e+05 *
    0.0000         0
    0.0000   -0.6588
    0.0001   -1.1450
    0.0001   -1.4738
    0.0001   -1.6600
    0.0001   -1.7188
    0.0001   -1.6650
    0.0002   -1.5138
    0.0002   -1.2800
    0.0002   -0.9788
    0.0002   -0.6250
    0.0002   -0.2338
    0.0003    0.1800
    0.0003    0.6012
    0.0003    1.0149
    0.0003    1.4062
    0.0003    1.7599
    0.0004    2.0612
    0.0004    2.2949
    0.0004    2.4462
    0.0004    2.4999
 
beam支反力:
   1.0e+04 *
    0.0001    0.6000
    0.0002    1.2479
    0.0043    0.3750

第7章 2d Frame单元静态分析

单元刚度矩阵推导

参见博文:【刚度矩阵推导】2d frame 单元

例题

alt text

Matlab 代码

 % ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>
% ! chpater 7  using 2d-frame element
 
clear ;
clc
format default
p=struct();
 
% define model parameters
 
p.E = 210000;   % SI(mm)
p.A = 200;
p.P = [15000 10e6]';
p.Iz=2e8;
 
% define elems and nodes
 
p.nodes=[0 0 ;  % node's coordinates (x,y)
         0 6000 ;
         6000 6000 ;
         6000 0;];
p.elems=[1 2; 2 3; 3 4;];
 
p.node_num=size(p.nodes,1);
p.elem_num=size(p.elems,1);
p.node_Coord_x=p.nodes(:,1);
p.node_Coord_y=p.nodes(:,2);
 
% define problem's dimension = each node's dof  
p.problem_dimension=3;
 
% global degree of freedom number
p.global_dof_num=3*p.node_num;
% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num,1);
p.node_forces=zeros(p.global_dof_num,1);
 
% define boundary conditions
% fixed dof 
p.fix_dof=[1 2 3  10 11 12]';
% load dof and amplitude
p.load_dof=[4 6]';
 
% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num);
 
% compute all ElemStiffnessMatrix and assembly stiff matrix
elem_global_dof_num=p.problem_dimension*size(p.elems,2);
p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);
 
for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 2d frame elem has 2 node,each node have 3 dof,
    % i-th node --> U_(3*i-2) and U_(3*i-1) and  U_(3*i) dof(displacement of node) 
    elem_dof=[connectivity(1)*3-2 connectivity(1)*3-1 connectivity(1)*3 connectivity(2)*3-2 connectivity(2)*3-1 connectivity(2)*3];
 
    node1_x=p.node_Coord_x(connectivity(1));
    node1_y=p.node_Coord_y(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2));
    node2_y=p.node_Coord_y(connectivity(2));
    
    elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
    l=(node2_x-node1_x)/elem_length;
    m=(node2_y-node1_y)/elem_length;
 
    % coordinate transform matrix: U_l=L*U_g
    L=[l m 0 0 0 0;
        -m l 0 0 0 0;
        0 0 1 0 0 0;
        0 0 0 l m 0;
        0 0 0 -m l 0;
        0 0 0 0 0 1;];
 
    % in local coordinates, the stiffness matrix of the frame element is obtained by com-
    % bination of the stiffness of the bar element and the Bernoulli beam element
    E=p.E;
    A=p.A;
    Iz=p.Iz;
 
    k_e_local=zeros(elem_global_dof_num,elem_global_dof_num);
 
    k_e_local(1,1)=(E*A)/elem_length;
    k_e_local(4,4)=(E*A)/elem_length;
    
    k_e_local(1,4)=-(E*A)/elem_length;
    k_e_local(4,1)=-(E*A)/elem_length;
    
    k_e_local(2,2)=12*E*Iz/(elem_length^3);
    k_e_local(2,3)=6*E*Iz/(elem_length^2);
    k_e_local(2,5)=12*E*Iz/(elem_length^3);
    k_e_local(2,6)=6*E*Iz/(elem_length^2);
    
    k_e_local(3,2)=6*E*Iz/(elem_length^2);
    k_e_local(5,2)=12*E*Iz/(elem_length^3);
    k_e_local(6,2)=6*E*Iz/(elem_length^2);
    
    k_e_local(3,3)=4*E*Iz/(elem_length);
    k_e_local(3,5)=-6*E*Iz/(elem_length*elem_length);
    k_e_local(3,6)=2*E*Iz/elem_length;
    
    k_e_local(6,3)=2*E*Iz/elem_length;
    k_e_local(5,3)=-6*E*Iz/(elem_length*elem_length);
    
    k_e_local(5,5)=12*E*Iz/(elem_length*elem_length*elem_length);
    k_e_local(5,6)=-6*E*Iz/(elem_length*elem_length);
    k_e_local(6,5)=-6*E*Iz/(elem_length*elem_length);
    
    k_e_local(6,6)=4*E*Iz/elem_length;
 
    k_e_global=L' *k_e_local *L;
    p.elemStiffs(:,:,i)=k_e_global;
 
    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e_global;
 
end
 
% apply boundary condition 
p.displacements(p.fix_dof)=0;
p.node_forces(p.load_dof)=p.P;
 
% solution KU=F
p = solutionStruct(p);
 
% show compute result
disp('displacements of dofs :')
disp([(1:p.global_dof_num)', p.displacements])
 
disp('forces of dofs :')
disp([(1:p.global_dof_num)',p.node_forces])
 
figure
axis equal
 
% plot undeformed mesh
draw_2Dtruss_mesh(p.nodes,p,'black')
 
deform_scale=500;
deformed_node_coordinate=[p.node_Coord_x+deform_scale*p.displacements(3*(1:p.node_num)-2), ...
                            p.node_Coord_y+deform_scale*p.displacements(3*(1:p.node_num)-1)];
% plot deformed mesh
draw_2Dtruss_mesh(deformed_node_coordinate,p,'red')
 
title('位移变形')
grid on

alt text

第八章 3d Frame单元静态分析

略

第九章 analysis of grids

略

第十章 Timoshenko beams 静态分析

与欧拉-伯努利梁公式不同,Timoshenko梁公式考虑了横向剪切变形.因此,它能够模拟较薄或较厚的梁.本章对Timoshenko梁的静力弯曲、自由振动和屈曲进行了分析.

10.Timoshenko梁的静力分析-Formulation

posted @ 2024-04-13 20:26 FE-有限元鹰阅读(392) 评论(1) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 【刚度矩阵推导】2d frame 单元

· 【ABAQUS文档笔记】ABAQUS刚体单元和可变形单元的review

· 《有限元分析与应用》

· 【旧文重发】MATLAB 通过函数封装一劳永逸地解决线性规划与运输问题的linprog的标准化操作（附MATLAB代码）

· 有限元代码

公告

昵称： FE-有限元鹰
园龄： 4年10个月
粉丝： 28
关注： 8

+加关注

2025年3月

日

一

二

三

四

五

六

	% MATLAB codes for Finite Element Analysis

	% Problem 1: 3 springs problem

	% clear memory
	clear ;

	% ? define vars
	P=10; % N
	ki=[1;1;1;1];

	% elementNodes: connections at elements
	elementNodes=[1 2;
	2 3;
	2 4];

	% numberElements: number of Elements
	numberElements=size(elementNodes,1);

	% numberNodes: number of nodes
	numberNodes=4;

	% for structure:
	% displacements: displacement vector N x 1
	% force : force vector N x 1
	% stiffness: stiffness matrix N x N
	displacements=zeros(numberNodes,1); %#ok<PREALL>
	force=zeros(numberNodes,1);
	stiffness=zeros(numberNodes);

	% applied load at node 2
	force(2)=P;
	% computation of the system stiffness matrix
	for e=1:numberElements
	% elementDof: element degrees of freedom (Dof)
	elementDof=elementNodes(e,:) ;
	stiffness(elementDof,elementDof)=...
	stiffness(elementDof,elementDof)+ki(e)*[1 -1;-1 1];
	end
	% boundary conditions and solution
	% prescribed dofs
	prescribedDof=[1;3;4];
	% free Dof : activeDof
	% setdiff() 函数-->求两个数组的差集
	% C = setdiff(A,B) 返回 A 中存在但 B 中不存在的数据,不包含重复项.C 是有序的.
	activeDof=setdiff((1:numberNodes)',prescribedDof);

	% solution
	% !KU=F --> U=K\F
	displacements=stiffness(activeDof,activeDof)\force(activeDof);

	% positioning all displacements
	displacements1=zeros(numberNodes,1);
	displacements1(activeDof)=displacements;

	% output displacements/reactions
	outputDisplacementsReactions(displacements1,stiffness,numberNodes,prescribedDof)

	function outputDisplacementsReactions(displacements,stiffness,GDof,prescribedDof)
	% output of displacements and reactions in
	% tabular form
	% GDof: total number of degrees of freedom of the problem
	% displacements
	disp('Displacements')
	%displacements=displacements1;
	jj=1:GDof; % format
	[jj' displacements]
	% reactions
	F=stiffness*displacements;
	reactions=F(prescribedDof);
	disp('reactions')
	[prescribedDof reactions]
	end

	% antonio ferreira 2008
	% clear memory
	clear all
	% E; modulus of elasticity
	% A: area of cross section
	% L: length of bar
	E = 30e6;A=1;EA=E*A; L = 90;
	% generation of coordinates and connectivities
	% numberElements: number of elements
	numberElements=3;
	% generation equal spaced coordinates
	nodeCoordinates=linspace(0,L,numberElements+1);
	xx=nodeCoordinates;
	% numberNodes: number of nodes
	numberNodes=size(nodeCoordinates,2);
	% elementNodes: connections at elements
	ii=1:numberElements;
	elementNodes(:,1)=ii;
	elementNodes(:,2)=ii+1;
	% for structure:
	% displacements: displacement vector
	% force : force vector
	% stiffness: stiffness matrix
	displacements=zeros(numberNodes,1); %#ok<PREALL>
	force=zeros(numberNodes,1);
	stiffness=zeros(numberNodes,numberNodes);

	% applied load at node 2
	force(2)=3000.0;
	% computation of the system stiffness matrix
	for e=1:numberElements;
	% elementDof: element degrees of freedom (Dof)
	elementDof=elementNodes(e,:) ;
	nn=length(elementDof);
	length_element=nodeCoordinates(elementDof(2))...
	-nodeCoordinates(elementDof(1));
	detJacobian=length_element/2;invJacobian=1/detJacobian;

	% central Gauss point (xi=0, weight W=2)
	[shape,naturalDerivatives]=shapeFunctionL2(0.0);
	Xderivatives=naturalDerivatives*invJacobian;
	% B matrix
	B=zeros(1,nn); B(1:nn) = Xderivatives(:);
	stiffness(elementDof,elementDof)=stiffness(elementDof,elementDof)+B'B2detJacobianEA;
	end

	% boundary conditions and solution
	% prescribed dofs
	fixedDof=find(xx==min(nodeCoordinates(:)) \| xx==max(nodeCoordinates(:)))';
	prescribedDof=[fixedDof]
	% free Dof : activeDof
	activeDof=setdiff([1:numberNodes]',[prescribedDof]);

	% solution
	GDof=numberNodes;
	displacements=solution(GDof,prescribedDof,stiffness,force);

	% output displacements/reactions
	outputDisplacementsReactions(displacements,stiffness,numberNodes,prescribedDof)

	disp(stiffness)
	disp(force)
	disp(displacements)

	function displacements=solution(GDof,prescribedDof,stiffness,force)
	% function to find solution in terms of global displacements
	activeDof=setdiff([1:GDof]', ...
	[prescribedDof]);
	U=stiffness(activeDof,activeDof)\force(activeDof);
	displacements=zeros(GDof,1);
	displacements(activeDof)=U;
	end

	function [shape,naturalDerivatives]=shapeFunctionL2(xi)
	% shape function and derivatives for L2 elements
	% shape : Shape functions
	% naturalDerivatives: derivatives w.r.t. xi
	% xi: natural coordinates (-1 ... +1)
	shape=([1-xi,1+xi]/2)';
	naturalDerivatives=[-1;1]/2;
	end

	function outputDisplacementsReactionsStructure(p)
	% output of displacements and reactions in
	% tabular form
	% GDof: total number of degrees of freedom of
	% the problem
	% displacements
	disp('Displacements')
	jj=1:p.GDof;
	[jj' p.displacements]
	% reactions
	F=p.stiffness*p.displacements;
	reactions=F(p.prescribedDof);
	disp('reactions')
	[p.prescribedDof reactions]

	function shapeL2=shapeFunctionL2Structure(xi)
	% shape function and derivatives for L2 elements
	% shape : Shape functions
	% naturalDerivatives: derivatives w.r.t. xi
	% xi: natural coordinates (-1 ... +1)
	shapeL2=struct()
	shapeL2.shape=([1-xi,1+xi]/2)';
	shapeL2.naturalDerivatives=[-1;1]/2;

	end % end function shapeFunctionL2

	function displacements = solutionStructure(p)
	% function to find solution in terms of global displacements
	activeDof=setdiff([1:p.GDof]', [p.prescribedDof]);
	U=p.stiffness(activeDof,activeDof)\p.force(activeDof);
	displacements=zeros(p.GDof,1);
	displacements(activeDof)=U;
	end


	********************************************************************************
	Probe Values Report, written on Tue Apr 16 21:43:02 2024

	Source
	-------

	ODB: C:/Users/Administrator/Documents/learnFEM/Matlb Fea Codes repo/fea-codes-of-matlab/chapter04/problem4/abaqus/Job-2.odb
	Step: Step-1
	Frame: Increment 1: Step Time = 1.000
	Variable: U, U1

	Probe values reported at nodes

	Part Instance Node ID Orig. Coords Def. Coords
	X Y Z X Y Z
	-------------------------------------------------------------------------------------------------------------------------
	PART-1-1 1 0. 0. 0. 345.178E-06 -1.32149E-03 0.
	PART-1-1 2 0. 10. 0. 0. 10. 0.
	PART-1-1 3 10. 10. 0. 10. 10. 0.
	PART-1-1 4 10. 0. 0. 10. 0. 0.


	Part Instance Node ID Attached elem U, U1 U, U2 RF, RF1 RF, RF2
	--------------------------------------------------------------------------------------------------------------
	PART-1-1 1 3 345.178E-06 -1.32149E-03 0. 0.
	PART-1-1 2 1 0. -7.92893E-33 0. 7.92893E+03
	PART-1-1 3 2 -2.07107E-33 -2.07107E-33 2.07107E+03 2.07107E+03
	PART-1-1 4 3 2.07107E-33 0. -2.07107E+03 -0.

	% ! problem4.m --> based on <<MATLAB codes for Finite Element Analysis>>

	clear ;
	clc
	format default
	p=struct();

	% define model parameters

	p.E = 30e6;
	p.A = 2;
	p.L = 10;
	p.P = -10000;

	% define elems and nodes

	p.nodes=[0 0 ; % node's coordinates (x,y)
	0 p.L ;
	p.L p.L ;
	p.L 0 ];
	p.elems=[1 2 ;
	1 3 ;
	1 4];

	p.node_num=size(p.nodes,1);
	p.elem_num=size(p.elems,1);
	p.node_Coord_x=p.nodes(:,1);
	p.node_Coord_y=p.nodes(:,2);

	% define problem's dimension : 1D/2D/3D;
	% = each node's dof
	p.problem_dimension=size(p.nodes,2);

	% global degree of freedom number
	p.global_dof_num=2*p.node_num;
	% define global Nodal displacement colum vector
	p.displacements=zeros(p.global_dof_num,1);
	p.node_forces=zeros(p.global_dof_num,1);

	% define boundary conditions
	% fixed dof
	p.fix_dof=[3 4 5 6 7 8]';
	% load dof and amplitude
	p.load_dof=2;

	% initial global stiffness matrix
	p.global_stiffness_matrix=zeros(p.global_dof_num);

	% compute all ElemStiffnessMatrix and assembly stiff matrix
	elem_global_dof_num=2*size(p.elems,2);
	p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 2d truss node has 2 node,each node have 2 dof,
	% i-th node --> U_(2i-1) and U_(2i) dof(displacement of node)
	elem_dof=[connectivity(1)2-1 connectivity(1)2 ...
	connectivity(2)2-1 connectivity(2)2];

	node1_x=p.node_Coord_x(connectivity(1));
	node1_y=p.node_Coord_y(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));
	node2_y=p.node_Coord_y(connectivity(2));

	elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
	l=(node2_x-node1_x)/elem_length;
	m=(node2_y-node1_y)/elem_length;

	k_e=(p.Ep.A/elem_length)...
	[ll lm -ll -lm;
	lm mm -lm -mm;
	-ll -lm ll lm;
	-lm -mm lm mm];
	p.elemStiffs(:,:,i)=k_e;

	% assemble stiffness matrix
	p.global_stiffness_matrix(elem_dof,elem_dof)=...
	p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;

	end

	% apply boundary condition
	p.displacements(p.fix_dof)=0;
	p.node_forces(p.load_dof)=p.P;

	% solution KU=F
	p = solutionStruct(p);

	% show compute result
	disp('displacements of dofs :')
	disp([(1:p.global_dof_num)', p.displacements])

	disp('forces of dofs :')
	disp([(1:p.global_dof_num)',p.node_forces])

	figure
	axis equal

	% plot undeformed mesh
	draw_2Dtruss_mesh(p.nodes,p)

	% plot deformed mesh
	deformed_node_coordinate=cal_deformed_nodes(p,500.0);
	draw_2Dtruss_mesh(deformed_node_coordinate,p)

	title('位移变形')
	grid on


	function new_p=solutionStruct(p)
	% free dof in structure
	free_dof=setdiff((1:p.global_dof_num)',p.fix_dof);

	% solute the equation :KU=F --> U=K\F
	U=p.global_stiffness_matrix(free_dof,free_dof)\p.node_forces(free_dof);

	p.displacements(free_dof)=U;
	% displacements=zeros(p.global_dof_num,1);

	p.node_forces=p.global_stiffness_matrix*p.displacements;

	new_p=p;
	end

	function deformed_node_coordinate=cal_deformed_nodes(p,scale)
	% initial
	deformed_node_coordinate=p.nodes;
	if p.problem_dimension==2
	index1=1:2:p.global_dof_num;
	index2=2:2:p.global_dof_num;
	deformed_node_coordinate(:,1)=deformed_node_coordinate(:,1)+p.displacements(index1)*scale;
	deformed_node_coordinate(:,2)=deformed_node_coordinate(:,2)+p.displacements(index2)*scale;
	elseif p.problem_dimension==3
	index1=1:3:p.global_dof_num;
	index2=2:3:p.global_dof_num;
	index3=3:3:p.global_dof_num;
	deformed_node_coordinate(:,1)=deformed_node_coordinate(:,1)+p.displacements(index1);
	deformed_node_coordinate(:,2)=deformed_node_coordinate(:,2)+p.displacements(index2);
	deformed_node_coordinate(:,3)=deformed_node_coordinate(:,3)+p.displacements(index3);
	end
	end

	function draw_2Dtruss_mesh(nodes,p)
	xy=zeros(size(p.elems,2)*p.elem_num,p.problem_dimension);
	cur1=1;
	cur2=0;
	for i =1:p.elem_num

	connectivity=p.elems(i,:);
	attach_node_num=size(connectivity,2);

	cur1=cur2+1;
	cur2=cur2+attach_node_num;

	coord=zeros(attach_node_num,p.problem_dimension);
	for j=1:size(connectivity,2)
	coord(j,:)=nodes(connectivity(j),:);
	end

	xy(cur1:cur2,:)=coord;
	end
	line(xy(:,1),xy(:,2),"Marker","o","LineWidth",1,"Color",rand(1,3))
	end


	% define model parameters
	% unit:SI(m)
	p.E = 70e9;
	p.A = 3e-4;
	p.L = 3;
	p.P = 1000*[-50;-100;-50];
	disp('单位制:SI(m)')
	% define elems and nodes

	p.nodes=[0 0 ; % node's coordinates (x,y)
	0 p.L ;
	p.L 0 ;
	p.L p.L ;
	2*p.L 0 ;
	2*p.L p.L];
	p.elems=[1 2 ;
	1 3;
	2 3;
	2 4;
	1 4;
	3 4;
	3 6;
	4 5;
	4 6;
	3 5;
	5 6];

	......

	% define boundary conditions
	% fixed dof
	p.fix_dof=[1 2 10]';
	% load dof and amplitude
	p.load_dof=[4 8 12]';

	......

	单位制:SI(m)
	displacements of dofs :
	1.0000 0
	2.0000 0
	3.0000 0.0071
	4.0000 -0.0090
	5.0000 0.0052
	6.0000 -0.0163
	7.0000 0.0052
	8.0000 -0.0201
	9.0000 0.0105
	10.0000 0
	11.0000 0.0034
	12.0000 -0.0090

	forces of dofs :
	1.0e+05 *

	0.0000 0.0000
	0.0000 1.0000
	0.0000 0.0000
	0.0000 -0.5000
	0.0001 -0.0000
	0.0001 -0.0000
	0.0001 -0.0000
	0.0001 -1.0000
	0.0001 0.0000
	0.0001 1.0000
	0.0001 0
	0.0001 -0.5000

	% ! problem7.m --> based on <<MATLAB codes for Finite Element Analysis>>

	clear ;
	clc
	format default
	p=struct();

	% define model parameters
	% unit:SI(m)
	p.E = 1.2e6;
	p.A = [0.302 0.729 0.187]';
	p.P = -1000;
	disp('单位制:SI(m)')
	% define elems and nodes

	p.nodes=[72 0 0 ; % node's coordinates (x,y,z)
	0 36 0;
	0 36 72;
	0 0 -48];

	p.elems=[1 2 ;
	1 3;
	1 4];

	p.node_num=size(p.nodes,1);
	p.elem_num=size(p.elems,1);
	p.node_Coord_x=p.nodes(:,1);
	p.node_Coord_y=p.nodes(:,2);
	p.node_Coord_z=p.nodes(:,3);

	% define problem's dimension : 1D/2D/3D;
	% = each node's dof
	p.problem_dimension=size(p.nodes,2);

	% global degree of freedom number
	p.global_dof_num=p.problem_dimension*p.node_num;
	% define global Nodal displacement colum vector
	p.displacements=zeros(p.global_dof_num,1);
	p.node_forces=zeros(p.global_dof_num,1);

	% define boundary conditions
	% fixed dof
	p.fix_dof=[2 4 5 6 7 8 9 10 11 12]';
	% load dof and amplitude
	p.load_dof=3;

	% initial global stiffness matrix
	p.global_stiffness_matrix=zeros(p.global_dof_num);

	% compute all ElemStiffnessMatrix and assembly stiff matrix
	elem_global_dof_num=p.problem_dimension*size(p.elems,2);
	p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 3d truss node has 2 node,each node have 3 dof,
	% i-th node --> U_(3i-2) / U_(3i-1) / U_(3*i) dof(displacement of node)
	elem_dof=[connectivity(1)3-2 connectivity(1)3-1 connectivity(1)*3 ...
	connectivity(2)3-2 connectivity(2)3-1 connectivity(2)*3];

	node1_x=p.node_Coord_x(connectivity(1));
	node1_y=p.node_Coord_y(connectivity(1));
	node1_z=p.node_Coord_z(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));
	node2_y=p.node_Coord_y(connectivity(2));
	node2_z=p.node_Coord_z(connectivity(2));

	elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2+(node1_z-node2_z).^2);
	Cx=(node2_x-node1_x)/elem_length;
	Cy=(node2_y-node1_y)/elem_length;
	Cz=(node2_z-node1_z)/elem_length;

	T=[CxCx CxCy Cx*Cz;
	CxCy CyCy Cy*Cz;
	CxCz CyCz Cz*Cz];

	k_e=(p.Ep.A(i)/elem_length)[T -T;-T T];
	p.elemStiffs(:,:,i)=k_e;

	% assemble stiffness matrix
	p.global_stiffness_matrix(elem_dof,elem_dof)=...
	p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;

	end

	% apply boundary condition
	p.displacements(p.fix_dof)=0;
	p.node_forces(p.load_dof)=p.P;

	% solution KU=F
	p = solutionStruct(p);

	% calculate stress
	p.stress=zeros(p.elem_num,1);
	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 3d truss node has 2 node,each node have 3 dof,
	% i-th node --> U_(3i-2) / U_(3i-1) / U_(3*i) dof(displacement of node)
	elem_dof=[connectivity(1)3-2 connectivity(1)3-1 connectivity(1)*3 ...
	connectivity(2)3-2 connectivity(2)3-1 connectivity(2)*3];

	node1_x=p.node_Coord_x(connectivity(1));
	node1_y=p.node_Coord_y(connectivity(1));
	node1_z=p.node_Coord_z(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));
	node2_y=p.node_Coord_y(connectivity(2));
	node2_z=p.node_Coord_z(connectivity(2));

	elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2+(node1_z-node2_z).^2);
	Cx=(node2_x-node1_x)/elem_length;
	Cy=(node2_y-node1_y)/elem_length;
	Cz=(node2_z-node1_z)/elem_length;

	p.stress(i)=(p.E/elem_length)[-Cx -Cy -Cz Cx Cy Cz]p.displacements(elem_dof);
	end

	% show compute result
	disp('displacements of dofs :')
	disp([(1:p.global_dof_num)', p.displacements])

	disp('forces of dofs :')
	disp([(1:p.global_dof_num)',p.node_forces])

	disp('3d truss elem stress :')
	disp([(1:p.elem_num)',p.stress])

	figure
	axis equal

	% plot undeformed mesh
	draw_3Dtruss_mesh(p.nodes,p)

	% plot deformed mesh
	deformed_node_coordinate=cal_deformed_nodes(p,50.0);
	draw_3Dtruss_mesh(deformed_node_coordinate,p)

	title('位移变形')
	grid on


	function draw_3Dtruss_mesh(nodes,p)
	xy=zeros(size(p.elems,2)*p.elem_num,p.problem_dimension);
	cur1=1;
	cur2=0;
	for i =1:p.elem_num

	connectivity=p.elems(i,:);
	attach_node_num=size(connectivity,2);

	cur1=cur2+1;
	cur2=cur2+attach_node_num;

	coord=zeros(attach_node_num,p.problem_dimension);
	for j=1:size(connectivity,2)
	coord(j,:)=nodes(connectivity(j),:);
	end

	xy(cur1:cur2,:)=coord;
	end
	line(xy(:,1),xy(:,2),xy(:,3),"Marker","o","LineWidth",1,"Color",rand(1,3))
	end

	单位制:SI(m)
	displacements of dofs :
	1.0000 -0.0711
	2.0000 0
	3.0000 -0.2662
	4.0000 0
	5.0000 0
	6.0000 0
	7.0000 0
	8.0000 0
	9.0000 0
	10.0000 0
	11.0000 0
	12.0000 0


	forces of dofs :
	1.0e+03 *

	0.0010 0
	0.0020 -0.2232
	0.0030 -1.0000
	0.0040 0.2561
	0.0050 -0.1281
	0.0060 0
	0.0070 -0.7024
	0.0080 0.3512
	0.0090 0.7024
	0.0100 0.4463
	0.0110 0
	0.0120 0.2976


	3d truss elem stress :
	1.0e+03 *

	0.0010 -0.9482
	0.0020 1.4454
	0.0030 -2.8685


	% define model parameters
	% unit:SI(mm)
	p.E = 210000;
	p.A = 100;
	p.P = -10000;
	disp('单位制:SI(mm)')
	% define elems and nodes

	p.nodes=[4 4 3 ; % node's coordinates (x,y,z)
	0 4 0;
	0 4 6;
	4 0 3;
	8 -1 1];
	p.elems=[1 2 ;
	1 3;
	1 4;
	1 5];

	......

	% define boundary conditions
	% fixed dof
	p.fix_dof=[4:15]';
	% load dof and amplitude
	p.load_dof=2;

	......

	单位制:SI(mm)
	displacements of dofs :
	1.0000 -0.0003
	2.0000 -0.0015
	3.0000 0.0003

	forces of dofs :
	1.0e+04 *

	0.0004 0.0271
	0.0005 0
	0.0006 0.0203
	0.0007 0.1355
	0.0008 0
	0.0009 -0.1016
	0.0010 0
	0.0011 0.7968
	0.0012 0
	0.0013 -0.1626
	0.0014 0.2032
	0.0015 0.0813

	3d truss elem stress :
	1.0000 -3.3865
	2.0000 -16.9326
	3.0000 -79.6809
	4.0000 -27.2610

	deformed nodes:
	3.8791 3.3929 3.1075
	0 4.0000 0
	0 4.0000 6.0000
	4.0000 0 3.0000
	8.0000 -1.0000 1.0000

	% ! problem9.m --> based on <<MATLAB codes for Finite Element Analysis>>

	clear ;
	clc
	format default
	p=struct();
	Case=1;

	% define model parameters
	% unit:SI(mm)
	p.EI = 1;
	p.L = 1;
	% * distribute force
	p.P = -1;
	% define elems num of bornoulli beam
	p.elem_num=80;
	% define problem's dimension : 1D/2D/3D;
	% = each node's dof
	p.problem_dimension=2;

	% ? -------------------------------

	% node's coordinates (x)
	p.nodes=linspace(0,p.L,p.elem_num+1)';
	p.elems=zeros(p.elem_num,2);
	for i=1:p.elem_num
	p.elems(i,1)=i;
	p.elems(i,2)=i+1;
	end
	p.node_num=size(p.nodes,1);
	p.node_Coord_x=p.nodes(:,1);

	% global degree of freedom number
	p.global_dof_num=2*p.node_num;

	% define global Nodal displacement colum vector
	p.displacements=zeros(p.global_dof_num,1);
	p.node_forces=zeros(p.global_dof_num,1);

	% define boundary conditions
	if Case==1
	% Case 01: clamped at x=0
	% fixed dof
	p.fix_dof=[1 2]';
	end
	if Case==2
	% Case 02: clamped-clamped
	% fixed dof
	p.fix_dof=[1 p.elem_num2+1 2 p.elem_num2+2]';
	end
	if Case==3
	% Case 03: simply supported-simply supported
	p.fix_dof=[1 2*p.elem_num+1]';
	end

	% initial global stiffness matrix
	p.global_stiffness_matrix=zeros(p.global_dof_num);

	% compute all ElemStiffnessMatrix and assembly stiff matrix
	elem_global_dof_num=p.problem_dimension*size(p.elems,2);
	p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
	% i-th node --> w_(2i-1) / dw_(2i-2) / dof(displacement of node)
	elem_dof=[connectivity(1)2-1 connectivity(1)2 ...
	connectivity(2)2-1 connectivity(2)2 ];

	node1_x=p.node_Coord_x(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));

	elem_length=node2_x-node1_x;
	a=0.5*elem_length;
	k_e=(p.EI/(a^3))0.5[3 3a -3 3a;
	3a 4aa -3a 2aa;
	-3 -3a 3 -3a;
	3a 2aa -3a 4aa];
	p.elemStiffs(:,:,i)=k_e;
	p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.Pa/3)[3 a 3 -a])';

	% assemble stiffness matrix
	p.global_stiffness_matrix(elem_dof,elem_dof)=...
	p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
	end

	% apply boundary condition
	p.displacements(p.fix_dof)=0;

	% solution KU=F
	p = solutionStruct_bornoulli_beam(p);

	plot(p.node_Coord_x,p.displacements(1:2:2*p.node_num),'.')
	minDispalce=min(p.displacements(1:2:2*p.node_num));
	text(0.35,0.5*minDispalce,"挠度最小值:"+minDispalce)

	if Case==1
	title("Case 01: clamped at x=0")
	end
	if Case==2
	title("Case 02: clamped-clamped")
	end
	if Case==3
	title("Case 03: simply supported-simply supported")
	end


	function new_p=solutionStruct_bornoulli_beam(p)
	% free dof in structure
	free_dof=setdiff((1:p.global_dof_num)',p.fix_dof);

	% solute the equation :KU=F --> U=K\F
	U=p.global_stiffness_matrix(free_dof,free_dof)\p.node_forces(free_dof);

	p.displacements(free_dof)=U;
	% displacements=zeros(p.global_dof_num,1);

	p.node_forces_fu=p.global_stiffness_matrix*p.displacements;

	new_p=p;
	end

	% ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>

	clear;
	clc;
	format compact;

	p=struct();

	% define model paraeters
	p.E=1000000;
	p.L=10;
	t=p.L/1000;
	p.Iz=(1*t^3)/12;
	p.P=-1000;
	% define elems of bornoulli beam
	p.elem_num=20;
	% dof of each nornoulli beam nodes
	p.bornoulli_node_dof_num=2;
	p.spring_k=100;

	% -------------------------------

	% node's coordinate
	p.nodes=linspace(0,p.L,p.elem_num+1)';
	p.elems=zeros(p.elem_num,2);
	for i=1:p.elem_num
	p.elems(i,1)=i;
	p.elems(i,2)=i+1;
	end
	p.node_num=size(p.nodes,1);
	p.node_Coord_x=p.nodes(:,1);

	% global degree of freedom number:
	% ! bornoulli_beam node's dof + 1 spring node's dof
	p.global_dof_num=2*p.node_num;

	% define global Nodal displacement colum vector
	p.displacements=zeros(p.global_dof_num+1,1);
	p.node_forces=zeros(p.global_dof_num+1,1);

	% define boundary conditions
	% ! BC: clamped at x=0
	p.fix_dof=[1,2,p.global_dof_num+1]';

	% initial global stiffness matrix
	p.global_stiffness_matrix=zeros(p.global_dof_num+1);

	% compute all ElemStiffnessMatrix and assembly stiff matrix
	beam_dof_num=p.bornoulli_node_dof_num*size(p.elems,2);
	p.elemStiffs=zeros(beam_dof_num,beam_dof_num,p.elem_num);

	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
	% i-th node --> w_(2i-1) / dw_(2i-2) / dof(displacement of node)
	elem_dof=[connectivity(1)2-1 connectivity(1)2 ...
	connectivity(2)2-1 connectivity(2)2 ];

	node1_x=p.node_Coord_x(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));

	elem_length=node2_x-node1_x;
	a=0.5*elem_length;
	k_e=(p.Ep.Iz/(a^3))0.5[3 3a -3 3*a;
	3a 4aa -3a 2aa;
	-3 -3a 3 -3a;
	3a 2aa -3a 4aa];
	p.elemStiffs(:,:,i)=k_e;
	p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.Pa/3)[3 a 3 -a])';

	% assemble stiffness matrix
	p.global_stiffness_matrix(elem_dof,elem_dof)=...
	p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
	end

	% add spring nodes into global_stiffness_matrix
	p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])=...
	p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])+p.spring_k*[1,-1;-1,1];

	% apply boundary condition
	p.node_forces(p.fix_dof)=0;

	% solution KU=F
	p=solutionStruct(p);

	% ! display result
	disp('Bornoulli beam result !')

	disp('beam挠度:')
	omega_index=(1:2:p.global_dof_num)';
	disp([omega_index,p.displacements(omega_index)])

	disp('beam转角:')
	theta_index=(2:2:p.global_dof_num)';
	disp([theta_index,p.displacements(theta_index)])

	disp('beam支反力:')
	disp([p.fix_dof,p.node_forces(p.fix_dof)])

	plot(p.node_Coord_x,p.displacements(omega_index),'-')
	title('Problem10.m')

	Bornoulli beam result !
	beam挠度:
	1.0e+05 *
	0.0000 0
	0.0000 -0.1722
	0.0001 -0.6300
	0.0001 -1.2909
	0.0001 -2.0800
	0.0001 -2.9297
	0.0001 -3.7800
	0.0001 -4.5785
	0.0002 -5.2801
	0.0002 -5.8473
	0.0002 -6.2501
	0.0002 -6.4661
	0.0003 -6.4802
	0.0003 -6.2849
	0.0003 -5.8802
	0.0003 -5.2737
	0.0003 -4.4803
	0.0003 -3.5225
	0.0004 -2.4303
	0.0004 -1.2413
	0.0004 -0.0004

	beam转角:
	1.0e+05 *
	0.0000 0
	0.0000 -0.6588
	0.0001 -1.1450
	0.0001 -1.4738
	0.0001 -1.6600
	0.0001 -1.7188
	0.0001 -1.6650
	0.0002 -1.5138
	0.0002 -1.2800
	0.0002 -0.9788
	0.0002 -0.6250
	0.0002 -0.2338
	0.0003 0.1800
	0.0003 0.6012
	0.0003 1.0149
	0.0003 1.4062
	0.0003 1.7599
	0.0004 2.0612
	0.0004 2.2949
	0.0004 2.4462
	0.0004 2.4999

	beam支反力:
	1.0e+04 *
	0.0001 0.6000
	0.0002 1.2479
	0.0043 0.3750

	% ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>
	% ! chpater 7 using 2d-frame element

	clear ;
	clc
	format default
	p=struct();

	% define model parameters

	p.E = 210000; % SI(mm)
	p.A = 200;
	p.P = [15000 10e6]';
	p.Iz=2e8;

	% define elems and nodes

	p.nodes=[0 0 ; % node's coordinates (x,y)
	0 6000 ;
	6000 6000 ;
	6000 0;];
	p.elems=[1 2; 2 3; 3 4;];

	p.node_num=size(p.nodes,1);
	p.elem_num=size(p.elems,1);
	p.node_Coord_x=p.nodes(:,1);
	p.node_Coord_y=p.nodes(:,2);

	% define problem's dimension = each node's dof
	p.problem_dimension=3;

	% global degree of freedom number
	p.global_dof_num=3*p.node_num;
	% define global Nodal displacement colum vector
	p.displacements=zeros(p.global_dof_num,1);
	p.node_forces=zeros(p.global_dof_num,1);

	% define boundary conditions
	% fixed dof
	p.fix_dof=[1 2 3 10 11 12]';
	% load dof and amplitude
	p.load_dof=[4 6]';

	% initial global stiffness matrix
	p.global_stiffness_matrix=zeros(p.global_dof_num);

	% compute all ElemStiffnessMatrix and assembly stiff matrix
	elem_global_dof_num=p.problem_dimension*size(p.elems,2);
	p.elemStiffs=zeros(elem_global_dof_num,elem_global_dof_num,p.elem_num);

	for i=1:p.elem_num
	connectivity=p.elems(i,:);
	% elem's all dof : a 2d frame elem has 2 node,each node have 3 dof,
	% i-th node --> U_(3i-2) and U_(3i-1) and U_(3*i) dof(displacement of node)
	elem_dof=[connectivity(1)3-2 connectivity(1)3-1 connectivity(1)3 connectivity(2)3-2 connectivity(2)3-1 connectivity(2)3];

	node1_x=p.node_Coord_x(connectivity(1));
	node1_y=p.node_Coord_y(connectivity(1));
	node2_x=p.node_Coord_x(connectivity(2));
	node2_y=p.node_Coord_y(connectivity(2));

	elem_length=sqrt((node1_x-node2_x).^2+(node1_y-node2_y).^2);
	l=(node2_x-node1_x)/elem_length;
	m=(node2_y-node1_y)/elem_length;

	% coordinate transform matrix: U_l=L*U_g
	L=[l m 0 0 0 0;
	-m l 0 0 0 0;
	0 0 1 0 0 0;
	0 0 0 l m 0;
	0 0 0 -m l 0;
	0 0 0 0 0 1;];

	% in local coordinates, the stiffness matrix of the frame element is obtained by com-
	% bination of the stiffness of the bar element and the Bernoulli beam element
	E=p.E;
	A=p.A;
	Iz=p.Iz;

	k_e_local=zeros(elem_global_dof_num,elem_global_dof_num);

	k_e_local(1,1)=(E*A)/elem_length;
	k_e_local(4,4)=(E*A)/elem_length;

	k_e_local(1,4)=-(E*A)/elem_length;
	k_e_local(4,1)=-(E*A)/elem_length;

	k_e_local(2,2)=12EIz/(elem_length^3);
	k_e_local(2,3)=6EIz/(elem_length^2);
	k_e_local(2,5)=12EIz/(elem_length^3);
	k_e_local(2,6)=6EIz/(elem_length^2);

	k_e_local(3,2)=6EIz/(elem_length^2);
	k_e_local(5,2)=12EIz/(elem_length^3);
	k_e_local(6,2)=6EIz/(elem_length^2);

	k_e_local(3,3)=4EIz/(elem_length);
	k_e_local(3,5)=-6EIz/(elem_length*elem_length);
	k_e_local(3,6)=2EIz/elem_length;

	k_e_local(6,3)=2EIz/elem_length;
	k_e_local(5,3)=-6EIz/(elem_length*elem_length);

	k_e_local(5,5)=12EIz/(elem_lengthelem_lengthelem_length);
	k_e_local(5,6)=-6EIz/(elem_length*elem_length);
	k_e_local(6,5)=-6EIz/(elem_length*elem_length);

	k_e_local(6,6)=4EIz/elem_length;

	k_e_global=L' k_e_local L;
	p.elemStiffs(:,:,i)=k_e_global;

	% assemble stiffness matrix
	p.global_stiffness_matrix(elem_dof,elem_dof)=...
	p.global_stiffness_matrix(elem_dof,elem_dof)+k_e_global;

	end

	% apply boundary condition
	p.displacements(p.fix_dof)=0;
	p.node_forces(p.load_dof)=p.P;

	% solution KU=F
	p = solutionStruct(p);

	% show compute result
	disp('displacements of dofs :')
	disp([(1:p.global_dof_num)', p.displacements])

	disp('forces of dofs :')
	disp([(1:p.global_dof_num)',p.node_forces])

	figure
	axis equal

	% plot undeformed mesh
	draw_2Dtruss_mesh(p.nodes,p,'black')

	deform_scale=500;
	deformed_node_coordinate=[p.node_Coord_x+deform_scalep.displacements(3(1:p.node_num)-2), ...
	p.node_Coord_y+deform_scalep.displacements(3(1:p.node_num)-1)];
	% plot deformed mesh
	draw_2Dtruss_mesh(deformed_node_coordinate,p,'red')

	title('位移变形')
	grid on

有限元鹰

人物禀常格，有始必有终。 年时俯仰过，功名宜速崇。

有限元方法[Matlab]-笔记

第一章:matlab基础

第二章: 离散系统

笔记、例题

Matlab代码

problem1.m

outputDisplacementsReactions.m

第三章:1D 杆单元

笔记

例题1

matlab 代码

problem2.m

solution.m

shapeFunctionL2.m

例题2

matlab 代码

problem3Structure.m

outputDisplacementsReactionsStructure.m

shapeFunctionL2Structure.m

solutionStructure.m

第四章 2D桁架分析

笔记:单元刚度矩阵推导

例题1

手算过程

Abaqus 模拟结果

matlab 代码

problem4.m

例题2

Abaqus 结果

matlab 计算

第五章 3D 桁架分析

笔记 - 3d truss 单元刚度矩阵

例题 1

matlab 计算

例题2

Abaqus 结果

Matlab 计算

第六章 2D 欧拉-伯努利梁分析

单元刚度矩阵推导

例题 1

Abaqus 计算

Matlab 计算

例题2

matlab 计算

第7章 2d Frame单元静态分析

单元刚度矩阵推导

例题

Matlab 代码

第八章 3d Frame单元静态分析

第九章 analysis of grids

第十章 Timoshenko beams 静态分析

10.Timoshenko梁的静力分析-Formulation

公告

搜索

积分与排名

合集 (4)

随笔分类 (156)

随笔档案 (127)

常用链接

阅读排行榜

人物禀常格，有始必有终。年时俯仰过，功名宜速崇。