几何非线性matlab有限元程序_大变形_大变形有限元_matlab有限元_几何非线性_matlab变形

% Problem : Geometric Nonlinear Analysis of Membranes % Clamped boundary condition is used. %-------------------------------------------------------------------------- % Code written by : Amit Patil % E-mail :Amit Patil, aspatil07@gmail.com %-------------------------------------------------------------------------- clear clc disp('Please wait Programme is under Run') %-------------------------------------------------------------------------- % % Geometrical and material properties of membrane %-------------------------------------------------------------------------- a = 1 ; % length of the membrane (along X-axes) b = 1 ; % breadth of the membrane (along Y-axes) E = 10920; % elastic modulus nu = 0.3; % Poisson's ratio t = 0.001 ; % membrane thickness %Number of mesh element in x and y direction Nx=10; Ny=10; %-------------------------------------------------------------------------- % Input data for nodal connectivity for each element %-------------------------------------------------------------------------- [coordinates, nodes,nel,nnode] = MeshRectangular(a,b,Nx,Ny) ; % for node connectivity counting starts from 1 towards +y axis and new row again start at y=0. %-------------------------------------------------------------------------- % Transverse uniform pressure on membrane %--------------------------------------------------------------------------see pdf for more details. Pfinal = -100 ; totalstep=400; deltaP=Pfinal/totalstep; % Setting force and displacement zero, as this is a equilibrium configuration. P=0; u=zeros(nnode,1); v=zeros(nnode,1); w=zeros(nnode,1); displacement=zeros(3*nnode,1); % The nodal displacement vector is a={u,v,w} nnel=4; % number of nodes per element ndof=3; % number of dofs per node sdof=nnode*ndof; % total system dofs edof=nnel*ndof; % degrees of freedom per element for loadstep=1:401 if loadstep==1 deltaP=-0.01; else deltaP=Pfinal/totalstep; end % For first iteration the increment in load is very small, as it is imp to get first solution correct. P=P+deltaP %-------------------------------------------------------------------------- % Boundary condition %-------------------------------------------------------------------------- typeBC = 'c-c-c-c' ; for updateintforce=1:100 % This is Newton-Raphson iterative method for nonlinear analysis. [stiffness,tangentstiffness,force,bcdof]=SDF(P,E,nu,t,coordinates,nodes,nel,nnel,ndof,sdof,edof,displacement,typeBC,loadstep,updateintforce); totaldof=1:sdof; activedof=setdiff(totaldof,bcdof); deltadisplacement=zeros(sdof,1); residual=stiffness(activedof,activedof)*displacement(activedof)-force(activedof); residualnorm=norm(residual) deltadisplacement(activedof) = -tangentstiffness(activedof,activedof)\residual; displacement=displacement+deltadisplacement; norm(deltadisplacement); u = displacement(1:3:sdof) ; v = displacement(2:3:sdof) ; w = displacement(3:3:sdof) ; if residualnorm<10^-3 break; end end % Maximum transverse displacement format long minw = min(w) store(loadstep+1,1)=P; store(loadstep+1,2)=minw; end %-------------------------------------------------------------------------- % Deformed Shape %-------------------------------------------------------------------------- x = coordinates(:,1) ; y = coordinates(:,2) ; f3 = figure ; set(f3,'name','Postprocessing','numbertitle','off') ; plot3(x,y,w,'.') ; title('membrane deformation') ; plot(-store(:,1),-store(:,2),'LineWidth',2); title('Equilibrium path') ; xlabel('-P') ylabel('-w') %-------------------------------------------------------------------------- % Contour Plots %-------------------------------------------------------------------------- PlotFieldonDefoMesh(coordinates,nodes,w,w) title('Profile of w on deformed Mesh') ; PlotFieldonMesh(coordinates,nodes,u) title('Profile of u on membrane') PlotFieldonMesh(coordinates,nodes,v) title('Profile of v on membrane') PlotFieldonMesh(coordinates,nodes,w) title('Profile of w on membrane')