The stochastic finite element method solver for a simple tension test  

Prof. Marcin Kami?ski, Ph.D., D.Sc.  

Chair of Mechanics of Materials, Technical University of ??d?,  

Al. Politechniki 6, 93-590 ??d?, POLAND,  

tel/fax 48-42-6313551, email: Marcin.Kaminski@p.lodz.pl  

http://kmm.p.lodz.pl/pracownicy/Marcin_Kaminski/index.html  

This Maple script enables to solve numerically the simple tension test of the linear elastic bar by the constant force P applied at its both edges. The solution is provided using the stochastic perturbation-based finite element method derived from the Taylor expansion of all random parameters in the problem. The approach illustrated below is adequate to the 10th order expansion, where the expected values, standard deviations and coefficients of variation of the tensioned edge displacement are derived analytically and can be computed according to the 2nd, 4th, 6th, 8th and 10th order expansions. The plot3d option is utilized to make a visualization of those moments at the particular nodal point of the mesh, whereas the entire methodology can be linked with the other FEM Maple programs as well. Theoretical considerations are provided in: Computers & Structures, Volume 85, Issue 10, May 2007, Generalized perturbation-based stochastic finite element method in elastostatics, by Marcin Kami?ski, Elsevier Ltd.  

Lodz, POLAND, October 2006   

Matrices and problem initialization  

 

Deterministic problem solution for the rank 11 matrices, 0th order eqns,  linear finite elements with a single degree of freedom. 

>

C:=matrix([[1,-1,0,0,0,0,0,0,0,0,0],[-1,2,-1,0,0,0,0,0,0,0,0],[0,-1,2,-1,0,0,0,0,0,0,0],[0,0,-1,2,-1,0,0,0,0,0,0],[0,0,0,-1,2,-1,0,0,0,0,0],[0,0,0,0,-1,2,-1,0,0,0,0],[0,0,0,0,0,-1,2,-1,0,0,0],[0,0,0,0,0,0,-1,2,-1,0,0],[0,0,0,0,0,0,0,-1,2,-1,0],[0,0,0,0,0,0,0,0,-1,2,-1],[0,0,0,0,0,0,0,0,0,-1,1]]): kk:=Y*A/L: K:= kk*C:  evalm(K);

 

(1.1)

 

>	ff:=vector([-1,0,0,0,0,0,0,0,0,0,1]): pp:=P: rhs0:=pp*ff:

 

Random quantity definition.  

>

b:=Y:

 

Solutions for up to the 10th order linear equations systems  

 

>	soln:=evalm(leastsqrs(K,rhs0,'optimize')); assign(q=soln);

 

(2.1)

 

1st order solution to the problem  

>	K1:=diff(kk,b)*C: evalm(K1): f1:=diff(pp,b)*ff:

 

>	rhs1:=f1-multiply(K1,q): evalm(rhs1):

 

>	q1:=evalm(leastsqrs(K,rhs1,'optimize'));

 

(2.2)

 

2nd order solution to the problem  

>	K2:=C*diff(kk,b,b): evalm(K2): f2:=diff(pp,b,b)*ff:

 

>	rhs2:=f2-2*multiply(K1,q1): evalm(rhs2):

 

>	q2:=evalm(leastsqrs(K,rhs2,'optimize'));

 

(2.3)

 

3rd order solution to the problem  

>	K3:=diff(kk,b,b,b)*C: evalm(K3): f3:=diff(pp,b,b,b)*ff:

 

>	rhs3:=f3-3*multiply(K1,q2):evalm(rhs3):

 

>	q3:=evalm(leastsqrs(K,rhs3,'optimize'));

 

(2.4)

 

4th order solution to the problem 

>	K4:=diff(kk,b,b,b,b)*C: evalm(K4): f4:=diff(pp,b,b,b,b)*ff:

 

>	rhs4:=f4-4*multiply(K1,q3):evalm(rhs4):

 

>	q4:=evalm(leastsqrs(K,rhs4,'optimize'));

 

(2.5)

 

5th order solution to the problem 

>	K5:=diff(kk,b,b,b,b,b)*C: evalm(K5): f5:=diff(pp,b,b,b,b,b)*ff:

 

>	rhs5:=f5-5*multiply(K1,q4):evalm(rhs5):

 

>	q5:=evalm(leastsqrs(K,rhs5,'optimize'));

 

(2.6)

 

6th order solution to the problem 

>	K6:=diff(kk,b,b,b,b,b,b)*C: evalm(K6): f6:=diff(pp,b,b,b,b,b,b)*ff:

 

>	rhs6:=f6-6*multiply(K1,q5):evalm(rhs6):

 

>	q6:=evalm(leastsqrs(K,rhs6,'optimize'));

 

(2.7)

 

7th order solution to the problem  

>	K7:=diff(kk,b,b,b,b,b,b,b)*C: evalm(K7): f7:=diff(pp,b,b,b,b,b,b,b)*ff:

 

>	rhs7:=f7-7*multiply(K1,q6):evalm(rhs7):

 

>	q7:=evalm(leastsqrs(K,rhs7,'optimize'));

 

(2.8)

 

8th order solution to the problem 

>	K8:=diff(kk,b,b,b,b,b,b,b,b)*C: evalm(K8): f8:=diff(pp,b,b,b,b,b,b,b,b)*ff:

 

>	rhs8:=f8-8*multiply(K1,q7):evalm(rhs8):

 

>	q8:=evalm(leastsqrs(K,rhs8,'optimize'));

 

(2.9)

 

9th order solution to the problem 

>	K9:=diff(kk,b,b,b,b,b,b,b,b,b)*C: evalm(K9): f9:=diff(pp,b,b,b,b,b,b,b,b,b)*ff:

 

>	rhs9:=f9-9*multiply(K1,q8):evalm(rhs9):

 

>	q9:=evalm(leastsqrs(K,rhs9,'optimize'));

 

(2.10)

 

10th order solution to the problem 

>	K10:=diff(kk,b,b,b,b,b,b,b,b,b,b)*C: evalm(K10): f10:=diff(pp,b,b,b,b,b,b,b,b,b,b)*ff:

 

>	rhs10:=f10-10*multiply(K1,q9):evalm(rhs10):

 

>	q10:=evalm(leastsqrs(K,rhs10,'optimize'));

 

(2.11)

 

Computations of the response probabilistic moments 

 

Input probabilistic data for the Gaussian variable b  

>	Y:=209e9: A:=1E-4: L:=1.0E-1: sigb:=alfab*Y: P:=10E5:

 

>	mi2b:=sigb^2:

 

>	mi4b:=3*sigb^4:

 

>	mi6b:=1*3*5*sigb^6:

 

>	mi8b:=1*3*5*7*sigb^8:

 

>	mi10b:=1*3*5*7*9*sigb^10:

 

Computations of probabilistic moments for the output  

Expected values in various orders of the perturbation analysis 

>	E2q:=evalm(q+1/(2!)*eps^2*q2*mi2b);

 

(3.1)

 

>	E4q:=evalm(q+1/(2!)*eps^2*q2*mi2b+1/(4!)*eps^4*q4*mi4b);

 

(3.2)

 

>	E6q:=evalm(q+1/(2!)*eps^2*q2*mi2b+1/(4!)*eps^4*q4*mi4b+1/(6!)*eps^6*q6*mi6b);

 

(3.3)

 

>	E8q:=evalm(q+1/(2!)*eps^2*q2*mi2b+1/(4!)*eps^4*q4*mi4b+1/(6!)*eps^6*q6*mi6b+1/(8!)*eps^8*q8*mi8b);

 

(3.4)

 

>	E10q:=evalm(q+1/(2!)*eps^2*q2*mi2b+1/(4!)*eps^4*q4*mi4b+1/(6!)*eps^6*q6*mi6b+1/(8!)*eps^8*q8*mi8b+1/(10!)*eps^10*q10*mi10b);

 

(3.5)

 

Standard deviations in various orders of perturbation analysis 

>	stdev2:=sqrt(eps^2*q1[11]*q1[11]*mi2b); evalm(stdev2);

 

(3.6)

 

(3.6)

 

>	stdev4:=sqrt(eps^2*q1[11]*q1[11]*mi2b+eps^4*mi4b*(1/4*q2[11]*q2[11]+1/3*q1[11]*q3[11])); evalm(stdev4);

 

(3.7)

 

(3.7)

 

>	stdev6:=sqrt(eps^2*q1[11]*q1[11]*mi2b+eps^4*mi4b*(1/4*q2[11]*q2[11]+1/3*q1[11]*q3[11])+eps^6*mi6b*(1/36*q3[11]*q3[11]+1/24*q4[11]*q2[11]+1/60*q5[11]*q1[11])); evalm(stdev6);

 

(3.8)

 

(3.8)

 

 

Plotting of probabilistic moments for the tensioned edge displacement  

 

    Plotting standard deviations 

>	plot3d([stdev2,stdev4],eps=0.8..1.2,alfab=0..0.3,title='standard_deviations_of_q',axes=BOXED,labels=[eps,alfab,stdev],orientation=[-35,60]);

 

 

>	plot3d(stdev6,eps=0.8..1.2,alfab=0..0.3,title='standard_deviations_of_q',axes=BOXED,labels=[eps,alfab,stdev],orientation=[-35,60]);

 

 

     Plotting coefficients of variation 

>	E2:=E2q[11]:E4:=E4q[11]:E6:=E6q[11]:E8:=E8q[11]:E10:=E10q[11]:

 

>	alfa2:=stdev2/E2:alfa4:=stdev4/E4:alfa6:=stdev6/E6:

 

>	plot3d([alfa2,alfa4],eps=0.8..1.2,alfab=0..0.3,title='coefficient_of_variation',axes=BOXED,labels=[eps,alfab,coefalfa],orientation=[-35,60]);

 

 

>	plot3d(alfa6,eps=0.8..1.2,alfab=0..0.3,title='coefficient_of_variation',axes=BOXED,labels=[eps,alfab,coefalfa],orientation=[-35,60]);

 

 

    Expectations for the displacements  

>	plot3d([E2,E4,E6,E8,E10],eps=0.8..1.2,alfab=0..0.3,axes=BOXED,title='Expectations_of_q');

 

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.