% FEATool Multiphysics axisymmetric stress-strain with custom equation parsing. close all, clear all, clc test_case = 2; % 0 = Single triangle test. % 1 = Pressurized hollow cylinder. % 2 = Pressurized hollow sphere. % 3 = Disc with fixed edge and central point load. sfun = 'sflag1'; p = 20e4; % Load force. nu = 0.3; E = 200e9; % Geometry and grid. switch( test_case ) case 0 % Single triangle test. fea.grid.p = [ 3 4 3 ; 0 0 1 ]; fea.grid.c = [1:3]'; fea.grid.a = zeros(3,1); fea.grid.s = 1; fea.grid.b = gridbdr( fea.grid.p, fea.grid.c, fea.grid.a ); case 1 % Pressurized hollow cylinder. a = 1; b = 2; fea.grid = rectgrid( 10, 10, [a b;0 1] ); case 2 % Pressurized hollow sphere. a = 1; b = 2; fea.grid = ringgrid(12,72,a,b); fea.grid = delcells( fea.grid, selcells( fea.grid, '(x<=eps) + (y<=eps)') ); case 3 % Disc with fixed edge and central point load. p = 10; nu = 1/3; E = 1e3; a = 0; b = 10; fea.grid = rectgrid( 20, 20, [a b;-0.5 0.5] ); end % Equation definitions. eqn_r = '- c1*( (1-nu)*ur_r + nu*wz_r + nu/r*( u_r - ur_t - wz_t ) - (1-nu)/r^2*u_t ) - c1*c2*( uz_z + wr_z ) = Fr'; eqn_z = '- c1*( nu*ur_z + nu/r*u_z + (1-nu)*wz_z ) - c1*c2*( uz_r + wr_r ) = Fz'; % fea struct. fea.sdim = { 'r' 'z' }; fea.dvar = { 'u' 'w' }; fea.sfun = { sfun sfun }; fea.eqn = parseeqn( { eqn_r eqn_z }, fea.dvar, fea.sdim ); % Constants and expressions. c1 = E / ((1+nu)*(1-2*nu) ); c2 = (1-2*nu) / 2; fea.expr = { 'c1' [] [] {[num2str(c1*2*pi),'*r']} ; 'c2' [] [] {c2} ; 'nu' [] [] {nu} ; 'Fr' [] [] {0} ; 'Fz' [] [] {0} }; % Boundary conditions. switch( test_case ) case 0 fea.bdr.d = { [] [] [] ; 0 [] [] }; fea.bdr.n = { [] [num2str(p*2*pi),'*r'] [] ; [] [] [] }; case 1 fea.bdr.d = { [] [] [] [] ; 0 0 0 0 }; fea.bdr.n = { [] [] [] [num2str(p*2*pi),'*r'] ; [] [] [] [] }; case 2 fea.bdr.d = { [] [] [] [0] ; [] [] [0] [] }; fea.bdr.n = { [num2str(p*2*pi),'*r*-nr'] [] [] [] ; [num2str(p*2*pi),'*r*-nz'] [] [] [] }; case 3 fea.bdr.d = { [] [] [] [] ; [] [] [] [] }; fea.bdr.n = { [] [] [] [] ; [] [] [] [-p] }; % Set point constraints. [~,ix] = min( (fea.grid.p(1,:)-b).^2 + (fea.grid.p(2,:)-0).^2 ); fea.pnt(1).index = ix; fea.pnt(1).type = 'constr'; fea.pnt(1).dvar = 2; fea.pnt(1).expr = 0; end % Solve. fea = parseprob( fea ); fea.sol.u = solvestat( fea, 'icub', 5 ); % Postprocessing. switch( test_case ) case 0 syms r z real positive B = [ -1 0 1 0 0 0 ; (4-r-z)/r 0 (-3+r)/r 0 z/r 0 ; 0 -1 0 0 0 1 ; -1 -1 0 1 1 0 ]; C = E/((1+nu)*(1-2*nu))*[ 1-nu nu nu 0 ; nu 1-nu nu 0 ; nu nu 1-nu 0 ; 0 0 0 (1-2*nu)/2 ]; K = simplify( B'*C'*B*r ); A = eval( 2*pi*int( int( K, z, 0, 4-r ), r, 3, 4 ) ); A([2 4],:) = 0; A(:,[2 4]) = 0; A(2,2) = 1; A(4,4) = 1; f = zeros(6,1); f([1 5]) = 3*pi*p; u_ref = A\f; err_norm = norm( fea.sol.u - u_ref([1 3 5 2 4 6]) ) case 1 n = 20; r = linspace(a,b,n); z = 0.5*ones(1,n); u_ref = p * a^2*(1+nu)*(b^2+r.^2*(1-2*nu)) ./ (E*(b^2-a^2)*r); u = evalexpr( 'u', [r;z], fea ); subplot(1,2,1) postplot( fea, 'surfexpr', 'u' ) title('r-displacement') subplot(1,2,2), hold on plot(u_ref,r,'r-') plot(u,r,'b.') legend('exact solution','computed solution') xlabel('r') grid on case 2 n = 20; r = linspace(a,b,n); z = zeros(1,n); u_ref = 1./(2*E*(b^3-a^3)*r.^2) .* (2*(p*a^3)*(1-2*nu)*r.^3+p*(1+nu)*b^3*a^3); u = evalexpr( 'u', [r;z], fea ); subplot(1,2,1) postplot( fea, 'surfexpr', 'sqrt(u^2+w^2)', 'arrowexpr', {'u' 'w'} ) title('computed displacement') subplot(1,2,2), hold on plot(u_ref,r,'r-') plot(u,r,'b.') title( 'radial displacement') legend('exact solution','computed solution') xlabel('r') grid on case 3 subplot(2,2,1) postplot( fea, 'surfexpr', 'u' ) title('computed r-displacement') subplot(2,2,2) D = E*1^3/(12*(1-nu^2)); u_ref = [num2str(-p),'/(8*pi*',num2str(D),')*(',num2str((3+nu)/(1+nu)-1),'-2*log(r/',num2str(b),'))*r*z']; postplot( fea, 'surfexpr', u_ref ) title('exact r-displacement') subplot(2,2,3) postplot( fea, 'surfexpr', 'w' ) title('computed z-displacement') subplot(2,2,4) D = E*1^3/(12*(1-nu^2)); u_ref = [num2str(-p),'/(16*pi*',num2str(D),')*(',num2str((3+nu)/(1+nu)),'*(',num2str(b^2),'-r^2)+2*r^2*log(r/',num2str(b),'))']; postplot( fea, 'surfexpr', u_ref ) title('exact z-displacement') end