- Home
- Documents
*1 09 - finite element method 09 - finite element method - density growth - implementation*

prev

next

of 16

View

223Download

0

Embed Size (px)

09 - finite element method09 - finite element method -density growth - implementation

finite element methodintegration point basedstaggered solutionloop over all time stepsglobal newton iterationloop over all elementsloop over all quadrature pointslocal newton iteration to determinedetermine element residual & partial derivativedetermine global residual and iterational matrixdeterminedetermine state of biological equilibrium

finite element methodintegration point basedstaggered solutionloop over all time stepsglobal newton iterationloop over all elementsloop over all quadrature pointslocal newton iterationdetermine element residual & tangentdetermine global residual and tangentdeterminedetermine state of biological equilibriumnlin_femcnst_lawnlin_femelement1upd_denscnst_lawnlin_femnlin_femelement1

finite element methodnlin_fem.m%%% loop over all load steps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for is = (nsteps+1):(nsteps+inpstep); iter = 0; residuum = 1;%%% global newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while residuum > tol iter=iter+1; R = zeros(ndof,1); K = sparse(ndof,ndof); e_spa = extr_dof(edof,dof); %%%%% loop over all elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for ie = 1:nel [Ke,Re,Ie] = element1(e_mat(ie,:),e_spa(ie,:),i_var(ie,:),mat); [K, R, I ] = assm_sys(edof(ie,:),K,Ke,R,Re,I,Ie); end%%%%% loop over all elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u_inc(:,2)=dt*u_pre(:,2); R = R - time*F_pre; dofold = dof; [dof,F] = solve_nr(K,R,dof,iter,u_inc); residuum= res_norm((dof-dofold),u_inc); end%%% global newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% time = time + dt; i_var = I; plot_int(e_spa,i_var,nel,is);end %%% loop over all load steps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

finite element methodintegration point baseddiscrete residual discrete residual residual of mechanical equilibrium/balance of momentumcheck in matlab!

finite element methodintegration point basedlinearization stiffness matrix / iteration matrix linearization of residual wrt nodal dofscheck in matlab!

finite element methodelement1.mfunction [Ke,Re,Ie]=element1(e_mat,e_spa,i_var,mat)%%% element stiffness matrix Ke, residual Re, internal variables Ie %%%% Ie = i_var;Re = zeros(8,1);Ke = zeros(8,8);nod=4; delta = eye(2);indx=[1;3;5;7]; ex_mat=e_mat(indx); indy=[2;4;6;8]; ey_mat=e_mat(indy); %%% integration points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%g1=0.577350269189626; w1=1;gp(:,1)=[-g1; g1;-g1; g1]; w(:,1)=[ w1; w1; w1; w1]; gp(:,2)=[-g1;-g1; g1; g1]; w(:,2)=[ w1; w1; w1; w1];wp=w(:,1).*w(:,2); xsi=gp(:,1); eta=gp(:,2);%%% shape functions and derivatives in isoparametric space %%%%%%%%%%%%%N(:,1)=(1-xsi).*(1-eta)/4; N(:,2)=(1+xsi).*(1-eta)/4;N(:,3)=(1+xsi).*(1+eta)/4; N(:,4)=(1-xsi).*(1+eta)/4;dNr(1:2:8 ,1)=-(1-eta)/4; dNr(1:2:8 ,2)= (1-eta)/4;dNr(1:2:8 ,3)= (1+eta)/4; dNr(1:2:8 ,4)=-(1+eta)/4;dNr(2:2:8+1,1)=-(1-xsi)/4; dNr(2:2:8+1,2)=-(1+xsi)/4;dNr(2:2:8+1,3)= (1+xsi)/4; dNr(2:2:8+1,4)= (1-xsi)/4;JT=dNr*[ex_mat;ey_mat]';%%% element stiffness matrix Ke, residual Re, internal variables Ie %%%%

- finite element methodelement1.m%%% loop over all integration points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for ip=1:4 indx=[2*ip-1; 2*ip]; detJ=det(JT(indx,:)); if detJ
finite element methodassm_sys.mfunction [K,R,I]=assm_sys(edof,K,Ke,R,Re,I,Ie)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% assemble local element contributions to global tangent & residual %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% input: edof = [ elem X1 Y1 X2 Y2 ] ... incidence matrix%%% Ke = [ ndof x ndof ] ... element tangent Ke%%% Re = [ ndof x 1] ... element residual Re%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output: K = [ 4 x 4 ] ... global tangent K%%% R = [ fx_1 fy_1 fx_2 fy_2] ... global residual R%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[nie,n]=size(edof);I(edof(:,1),:)=Ie(:);t=edof(:,2:n);for i = 1:nie K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke; R(t(i,:)) =R(t(i,:)) +Re;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

finite element methodintegration point basedconstitutive equations stress calculation @ integration point level constitutive equations - given calculatecheck in matlab!

finite element methodintegration point basedconstitutive equations tangent operator / constitutive moduli linearization of stress wrt deformation gradientcheck in matlab!

finite element methodcnst_law.mfunction [A,P,var]=cnst_law(F,var,mat)%%% determine tangent, stress and internal variable %%%%%%%%%%%%%%%%%%%emod = mat(1); nue = mat(2); rho0 = mat(3); psi0 = mat(4); expm = mat(5); expn = mat(6); dt = mat(7);xmu = emod / 2.0 / (1.0+nue);xlm = emod * nue / (1.0+nue) / ( 1.0-2.0*nue ); F_inv = inv(F); J = det(F); delta = [1 0; 0 1];%%% update internal variable %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[var,facs,fact]=upd_dens(F,var,mat); %%% first piola kirchhoff stress %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%P = facs * (xmu * F + (xlm * log(J) - xmu) * F_inv'); %%% tangent %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i=1:2for j=1:2for k=1:2for l=1:2 A(i,j,k,l) = xlm * F_inv(j,i)*F_inv(l,k) ... - (xlm * log(J) - xmu) * F_inv(l,i)*F_inv(j,k) ... + xmu * delta(i,k)*delta(j,l); A(i,j,k,l) = facs * A(i,j,k,l) + fact * P(i,j)*P(k,l); end, end, end, end%%% determine tangent, stress and internal variable %%%%%%%%%%%%%%%%%%%

finite element methodintegration point basedconstitutive equationscheck in matlab! residual of biological equilibrium / balance of mass discrete density update

finite element methodupd_dens.mfunction [var,facs,fact]=upd_dens(F,var,mat)%%% update internal variable density %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tol = 1e-8; var = 0.0;xmu = emod / 2.0 / (1.0+nue);xlm = emod * nue / (1.0+nue) / ( 1.0-2.0*nue ); J = det(F);C = F'*F;I1 = trace(C);psi0_neo = xlm/2 * log(J)^2 + xmu/2 *(I1 - 2 - 2*log(J)); rho_k0 = (1+var)*rho0; rho_k1 = (1+var)*rho0; iter = 0; res = 1;%%% local newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while abs(res) > tol iter=iter+1; res =((rho_k1/rho0)^(expn-expm)*psi0_neo-psi0)*dt-rho_k1+rho_k0; dres= (expn-expm)*(rho_k1/rho0)^(expn-expm)*psi0_neo*dt/rho_k1-1; drho=- res/dres; rho_k1 = rho_k1+drho; end%%% local newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rho = rho_k1; var = rho / rho0 - 1;facs = (rho/rho0)^expn;facr = 1/dt - (expn-expm) * (rho/rho0)^(expn-expm) / rho *psi0_neo;fact = expn / rho * (rho/rho0)^( -expm) / facr;%%% update internal variable density %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

finite element methodex_frame.mfunction [q0,edof,bc,F_ext,mat,nel,node,ndof,nip] = ex_frame%%% input data for frame example %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%emod = 1000; nue = 0.3; rho0 = 1.0; psi0 = 1.0;expm = 3.0; expn = 2.0; dt = 1.0;mat = [emod,nue,rho0,psi0,expm,expn,dt];

xbox(1) = 0.0; xbox(2) = 2.0; nx = 8; ybox(1) = 0.0; ybox(2) = 1.0; ny = 4;[q0,edof] = mesh_sqr(xbox,ybox,nx,ny); [nel, sizen] = size(edof);[ndof,sizen] = size(q0); node = ndof/2; nip = 4;

%%% dirichlet boundary conditionsbc(1,1) = 2*(ny+1)*(0 ) +2; bc(1,2) = 0;bc(2,1) = 2*(ny+1)*(nx ) +2; bc(2,2) = 0; bc(3,1) = 2*(ny+1)*(nx/2+0) +1; bc(3,2) = 0; bc(4,1) = 2*(ny+1)*(nx/2+1) ; bc(4,2) = -ybox(2)/50;

%%% neumann boundary conditionsF_ext = zeros(ndof,1);%%% input data for frame example %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

finite element methodex_frame.m