Ramakrishnan

Dr. Ramakrishnan Vaidyanathan

279 Reputation

5 Badges

5 years, 53 days

Social Networks and Content at Maplesoft.com

I have retired as Professor-Mechanical in Sri Venkateswara College of Engineering and Technology under Anna University affiliated colleges in Tamil Nadu, India. I have 19 years of Industrial and 20 years of teaching experience. I am learning Maple for the past four and half years hoping to make at least one appreciable maple presentation.

MaplePrimes Activity


These are replies submitted by Ramakrishnan

@Preben Alsholm 

Thank you very much

@Preben Alsholm 

I am just adding my another attempted doc for experts to help me state 'follow the rules'. Maple will not accept the pi as a number in the axis range for a plot.

I tried to enter 4*pi as one of the axis range in plot prperties. Maple says it is not a number and in my opinion it has the right to say so because of so many restrictions it has with alphabets pi as a number as well as NAN.

sin(2*x)/cos(x+3*Pi*(1/2)) = 1

sin(2*x)/sin(x) = 1

(1)

"(->)"

[[x = (1/3)*Pi], [x = -(1/3)*Pi]]

(2)

sin(2*x)

sin(2*x)

(3)

"->"

 

cos(x+3*Pi*(1/2))

sin(x)

(4)

``


Download Try_Axis_range.mwTry_Axis_range.mw

@Markiyan Hirnyk 

Contrary assumptions means we have made impossible assumptions like less tha a parameter 10 (say X) and greater than a vale ten (or say Y, both X and Y unknown to me). So this is not a bug, a clear statement of what maple feels about our assumptions. Experts have maade rule for maple to follow and maple follows meticulously!). My attempt in a document would explain probably with less confusion and more clarity. Thanks for not throwing stones at me.

Regards;

NULL

We may see pi value will not be considered unless we evaluate its value.

evalb(-1 < -2)

false

(1)

``

evalb(5 > Pi)

Pi < 5

(2)

evalb(5 < evalf(Pi))

false

(3)

a := -5*Pi; b := -4*Pi

evalb(a < b)

-5*Pi < -4*Pi

(4)

The pi value is to be evaluated for maple to recognize, else it takes as undeclared variable (value uknown not zero) and hence declares we have made contrary assumptions.


Download Maple_suggests.mwMaple_suggests.mw

@tomleslie 

Thank you very much. At later stage while doing dynamic analysis, i shall follow the partial derivative terms. I understood the importance of clear communication with maple. Thanks again for the useful tip.I have gone through and saw the maple error free statements for your commands. Thanks again for the useful tips.

Ramakrishnan V

@Preben Alsholm 

It just says 'i dont understand!! In the limits, you give only limits, no recursive calculations, first multiplying and then dividing!! We do it separately and maple comes forward to give best answers for better problems. Following is my sincere attempt to gat help and maple obliges.

Hope my document gives clearer perspective.

``

solve([sin(2*x)/cos(x+3*Pi*(1/2)) = 1, x > -4*Pi, x < -2.5*Pi], x, allsolutions, explicit);

{x = -11.51917307}

(1)

solve([sin(2*x)/cos(x+3*Pi*(1/2)) = 1, x > -4*Pi, x < -(5*.5)*Pi], x, allsolutions, explicit);

{x = -11.51917307}

(2)

solve([sin(2*x)/cos(x+3*Pi*(1/2)) = 1, x > -4*Pi, x < -(5/2)*Pi], x, allsolutions, explicit);

Error, (in assume) contradictory assumptions

 

a := 5/2

5/2

(3)

solve([sin(2*x)/cos(x+3*Pi*(1/2)) = 1, x > -4*Pi, x < a*Pi], x, allsolutions, explicit);

{x = arctan(-3^(1/2))-2*Pi}, {x = arctan(-3^(1/2))}, {x = arctan(-3^(1/2))+2*Pi}, {x = arctan(3^(1/2))-4*Pi}, {x = arctan(3^(1/2))-2*Pi}, {x = arctan(3^(1/2))}, {x = arctan(3^(1/2))+2*Pi}

(4)

``


Download Maple_suggests.mwMaple_suggests.mw

 

@EugeneKalentev 

I have found out the answer  thanks to Professor Robert Lopez.

Since l__e was redefined in the line where q__0 was defined, % is taken. The doubt was clarified by Professor Robert Lopez.

My sincere thanks go on record for him for the immediate responses from him whenever i send a mail.

Thanks to you also for the attempt. Since I have not used dimensions in the command mode, problem was not from units and dimensions. Cheers.

Ramakrishnan V

@Lenin Araujo Castillo 

Thanks.I will try maplesim. Hope to use this later. I am learning but difficult without a book it seems. Some basic learning from class room is essential. Thanks for your response. Regards.

Ramakrishnan V

@Markiyan Hirnyk 

Maple 2015.2 is the one i use and Windows10 is in my PC. I shall fill up the product(s) data in my questionnaire hereafter. Thanks for letting me know my mistakes,however small it may look from one end. Regards.

Ramakrishnan V

@Markiyan Hirnyk 

Hope the attached equation and solution for the heat transfer problems would help understand the method.There may be unprofessional use of commands and manual mode since i am only learning 

 

# Determine the temperature distribution in a fin of circular cross section shown in Fig. Base of the fin is maintained at 100 0C and convection occurs at the tip:

 


restart

h := .2:    W/m2K; k := 2:W/cmK:   D__1 := 1:cm; D__2 := 1:  cm;

``phi__f := 20:0C;
A__1 := 22*D__1^2/(7*4) = 11/14``

A__2 := 22*D__2^2/(7*4) = 11/14``

l__1 := 4 cm:p := (22/7)*D__1 = 22/7 cm
l__2 := 4
 = 4cm:p := (22/7)*D__2 = 22/7NULL cm

Discretisation:

Number of linear one dimensional elements = 2;
Total Number of nodes = 3;
Global stiffness matrix is 3 x 3
For element 1
Stiffness matrix for conduction:
                                            Vector[row](2, {(1) = 1, (2) = 2})
K__1a := k*A__1/l__1.(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1})) = Matrix([[11/28, -11/28], [-11/28, 11/28]])  Vector(2, {(1) = 1, (2) = 2})

 

Stiffness matrix for convection:
                                            Vector[row](2, {(1) = 1, (2) = 2})
K__1b := (1/6)*h*p*l__1.(Matrix(2, 2, {(1, 1) = 2, (1, 2) = 1, (2, 1) = 1, (2, 2) = 2})) = Matrix([[.838095238000000, .419047619000000], [.419047619000000, .838095238000000]])  Vector(2, {(1) = 1, (2) = 2})

                                          Vector[row](2, {(1) = 1, (2) = 2})
K__1 := K__1a+K__1b = Matrix([[1.23095238085714, 0.261904761428571e-1], [0.261904761428571e-1, 1.23095238085714]]) Vector(2, {(1) = 1, (2) = 2})
``

For element 2
Stiffness matrix for conduction:
                                            Vector[row](2, {(1) = 2, (2) = 3})
K__2a := k*A__1/l__1.(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1})) = Matrix([[11/28, -11/28], [-11/28, 11/28]])  Vector(2, {(1) = 2, (2) = 3})

 

Stiffness matrix for convection:
                                            Vector[row](2, {(1) = 2, (2) = 3})
K__2b := (1/6)*h*p*l__1.(Matrix(2, 2, {(1, 1) = 2, (1, 2) = 1, (2, 1) = 1, (2, 2) = 2})) = Matrix([[.838095238000000, .419047619000000], [.419047619000000, .838095238000000]])  Vector(2, {(1) = 2, (2) = 3})

Stiffness matrix for convection:

 

  K__2c := h*A__2.(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})) = Matrix([[0., 0.], [0., .157142857100000]])  Vector(2, {(1) = 2, (2) = 3})

NULL
                                          Vector[row](2, {(1) = 1, (2) = 2})
K__2 := K__2a+K__2b+K__2c = Matrix([[1.23095238085714, 0.261904761428571e-1], [0.261904761428571e-1, 1.38809523795714]]) Vector(2, {(1) = 1, (2) = 2})

NULL                                                                       Vector[row](3, {(1) = 1, (2) = 2, (3) = 3})
K__g := Matrix(3, 3, {(1, 1) = K__1(1), (1, 2) = K__1(3), (1, 3) = 0, (2, 1) = K__1(2), (2, 2) = K__1(4)+K__2(1), (2, 3) = K__2(3), (3, 1) = 0, (3, 2) = K__2(2), (3, 3) = K__2(4)}) = Matrix([[K__1(1), K__1(3), 0], [K__1(2), K__1(4)+K__2(1), K__2(3)], [0, K__2(2), K__2(4)]])``  Vector(3, {(1) = 1, (2) = 2, (3) = 3}) 

Load vector:
for element 1
     
F__1 := (1/2)*h*p*l__1*`&phi;__f`.(Vector(2, {(1) = 1, (2) = 1})) = Vector[column]([[25.1428571400000], [25.1428571400000]]) Vector(2, {(1) = 1, (2) = 2})

``

for element 2

F__2 := (1/2)*h*p*l__2*`&phi;__f`.(Vector(2, {(1) = 1, (2) = 1}))+h*A__2*`&phi;__f`.(Vector(2, {(1) = 0, (2) = -1})) = Vector[column]([[25.1428571400000], [21.9999999980000]]) Vector(2, {(1) = 2, (2) = 3})

F__g := eval(Vector(3, {(1) = F__1(1), (2) = F__1(2)+F__2(1), (3) = F__2(2)})) = Vector[column]([[F__1(1)], [47.6666666657143], [F__2(2)]])`` Vector(3, {(1) = 1, (2) = 2, (3) = 3}) 
Tip insulation does not contribute to any load vector since heat transferred is zero.
Assembling matrices give the system equation

K__g, (Vector(3, {(1) = `&phi;__1`, (2) = `&phi;__2`, (3) = `&phi;__3`})) = F__g = Matrix(3, 3, {(1, 1) = K__1(1), (1, 2) = K__1(3), (1, 3) = 0, (2, 1) = K__1(2), (2, 2) = K__1(4)+K__2(1), (2, 3) = K__2(3), (3, 1) = 0, (3, 2) = K__2(2), (3, 3) = K__2(4)}), (Vector(3, {(1) = `#msub(mi("&phi;",fontstyle = "normal"),mi("1"))`, (2) = `#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`, (3) = `#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`})) = (Vector(3, {(1) = F__1(1), (2) = HFloat(47.66666666571429), (3) = F__2(2)}))

 

# Base of the fin is at 1000C;  phi__1 := 100:

Multiply the first column by `&phi;__1` to get the modified column.

We know phi__1 := 100: 

   = Substituting and eliminating the first row and column and transferring the modified first column to RHS

``

 

F__g[2] := -phi__1*K__g[2, 1]+F__g[2]

HFloat(47.66666666571429)

(1)

(Matrix(2, 2, {(1, 1) = K__g[2, 2], (1, 2) = K__g[2, 3], (2, 1) = K__g[3, 2], (2, 2) = K__g[3, 3]})).(Vector(2, {(1) = `&phi;__2`, (2) = `&phi;__3`})) = (Vector(2, {(1) = F__g[2], (2) = F__g[3]})) = (Vector(2, {(1) = 2.4619*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`+0.262e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`, (2) = 0.262e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`+1.3881*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`})) = (Vector(2, {(1) = F__g[2], (2) = F__g[3]}))``

  Solving for `&phi;__1`,`&phi;__2` 

NULL

solve({0.262e-1*phi__2+1.3881*phi__3 = 22.000, 2.4619*phi__2+0.262e-1*phi__3 = 45.0476}) = {phi__2 = 18.13287427, phi__3 = 15.50674929}NULL

 

0.262e-1*16.98+1.3881*20.05 = 28.276281

``

NULL

  # Determine the temperatures at the nodal interfaces for the two layered wall shown in Fig. The left face is supplied with heat flux of `q"` := 5 W/cm2;

  # and right face is maintained at phi__3 := 20 0C;

``

restart

 

``

````

``

`q"__1` := 5:  W/m2K; k__1 := .2:W/mK;  = W/mKA := 1:cm2;
We shall consider two one dimensional linear elements
For element 1
l__12 := 4: cm;

                                                       Vector[row](2, {(1) = 1, (2) = 2})
[K] = K__1 := k__1*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))/l__12 = Matrix([[0.500000000000000e-1, -0.500000000000000e-1], [-0.500000000000000e-1, 0.500000000000000e-1]]) 
Vector(2, {(1) = 1, (2) = 2})

[F]1= F__1 := `q"`*A.(Vector(2, {(1) = 1, (2) = 0})) = `q"`.(Vector(2, {(1) = 1, (2) = 0})) Vector(2, {(1) = 1, (2) = 2})NULL

or element 2
l__23 := 2: cm;

                                                      Vector[row](2, {(1) = 2, (2) = 3})
[K] = K__2 := k__2*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))/l__23 = Matrix([[0.300000000000000e-1, -0.300000000000000e-1], [-0.300000000000000e-1, 0.300000000000000e-1]]) Vector(2, {(1) = 2, (2) = 3})

 

[F]1= F__2 := `q"`*A.(Vector(2, {(1) = 0, (2) = 0})) = `q"`.(Vector(2, {(1) = 0, (2) = 0})) Vector(2, {(1) = 2, (2) = 3})

Assembling the matrices, we get the global system equation,

Global Stiffness matrix

           Vector[row](3, {(1) = 1, (2) = 2, (3) = 3})
K__g := Matrix(3, 3, {(1, 1) = K__1(1), (1, 2) = K__1(3), (1, 3) = 0, (2, 1) = K__1(2), (2, 2) = K__1(4)+K__2(1), (2, 3) = K__2(3), (3, 1) = 0, (3, 2) = K__2(2), (3, 3) = K__2(4)})
 = Matrix([[K__1(1), K__1(3), 0], [K__1(2), K__1(4)+K__2(1), K__2(3)], [0, K__2(2), K__2(4)]])``

Global load matrix
 = F__g := Vector(3, {(1) = F__1(1), (2) = F__1(2)+F__2(1), (3) = F__2(2)}) = Vector[column]([[F__1(1)], [F__1(2)+F__2(1)], [F__2(2)]])``

K__g.(Vector(3, {(1) = `&phi;__1`, (2) = `&phi;__2`, (3) = `&phi;__3`})) = F__g = (Vector(3, {(1) = 5.0000-0.500e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`, (2) = -5.0000+0.800e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`-0.300e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`, (3) = -0.300e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`+0.300e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`})) = (Vector(3, {(1) = F__1(1), (2) = F__1(2)+F__2(1), (3) = F__2(2)}))``

We know phi__3 := 20;Substituting and solving for `&phi;__1`,`&phi;__2` 

``

solve({-0.500e-1*u__1+0.800e-1*u__2 = .6, 0.500e-1*u__1-0.500e-1*u__2 = 5})

{u__1 = 286.6666667, u__2 = 186.6666667}

(2)

`&phi;__1` = 286.67; `&phi;__2` = 186.670C

NULL

``


Download 1Dheattransfersolutions.mw1Dheattransfersolutions.mw

maple.

@Markiyan Hirnyk 

I am sorry if I am not confusing. I know a reduction method in which the second order differential equations are converted to lower order and solvd. Galerkin and Ritz are the contributers. I attach with a doc detailing this methods with complete answers. Hope this will be useful. Ramakrishnna V

Consider a cantilever beam subjected to uniformly distributed load q__0as shown in Fig. Calculate the displacement and compare with exact solution.

  restart
The governing differential equation is  I*E*(diff(v(x), x, x, x, x))-q__0 = 0; Boundary conditions are v(0) = 0; (diff(v(x), x))*0 = 0; (diff(v(x), x, x))(L) = 0; (diff(v(x), x, x, x))(L) = 0;
Boundary conditions v(0) = 0; (diff(v(x), x))*0 = 0These conditions enforce zero displacement and slope at fixed end;
Boundary conditions (diff(v(x), x, x))(L) = 0; (diff(v(x), x, x, x))(L) = 0These conditions enforce zero bending moment and shear force at free end;
Let us find the approximate solution using single continuous trial function.
Step 1: Assumption of trial solution (guess solution)
"v(x) = conjugate(v)(x) = `c__0` + `c__1`*x + `c__2`*x^(2) + `c__3`*x^(3) + `c__4`*x^(4)"
The constants c__0, c__1 are determined using boundary conditions v(0) = 0; (diff(v(x), x))*0 = 0; c__0 = c__1 and c__1 = 0
The constants c__2, c__3 are determined in terms of c__4 using boundary conditions (diff(v(x), x, x))(L) = 0; (diff(v(x), x, x, x))(L) = 0;
c__2 = -6*L^2*c__4-3*L*c__3; c__3 = -4*c__4*L
Substituting the values for c__0, c__1, c__2, c__3, c__4 in the trial solution, we get
v(x) = c__4*(6*L^2*x^2-4*L*x^3+x^4) [Note that finding the trial function satisfying all the boundary conditions is thus cumbursome]
Step 2: Find the domain residual.
Substituting the trial solution in the governing equation, we get
    diff(c__4*(6*L^2*x^2-4*L*x^3+x^4), x, x, x, x) = 24*c__4

 

The governing equation is
R__d = I*E*(diff(v(x), x, x, x, x))-q__0 and I*E*(diff(v(x), x, x, x, x))-q__0 = 24*EI*c__4-q__0;
Step 3: Minimise the residual.24*EI*c__4-q__0 = 0"(->)"[[c__4 = (1/24)*q__0/EI]]NULL 
The final solution is v(x) = 0.4167e-1*q__0*(6*L^2*x^2-4*L*x^3+x^4)/EI

``

The following differential equation is available for a physical phenomenon.
  diff(y(x), x, x)+50 = 0, 0 <= x and x <= 10
Trial function is y = a__1*x*(10-x)
Boundary conditions are y(0) = 0; y(10) = 0
 Find the solution byGalerkin method

``

  unassign('a__1')

Step 1: Assumption of trial solution (guess solution)

y(x) = ((D@@2)(y))(x) and ((D@@2)(y))(x) = a__1*x*(10-x) and a__1*x*(10-x) = -a__1*x^2+10*a__1*x

diff(-a__1*x^2+10*a__1*x, x, x) = -2*a__1NULL

Step 2: Find the residue by substituting in the LHS of the governing equation
The governing equation is
diff(y(x), x, x)+50 = 0"->"R := -2*a__1+50 = -2*a__1+50:
The condition for Galerkin's method
Weighing function = Trial function (solution guessed)

w := y(x) = y(x) 

   

w*R = y(x)*(-2*a__1+50) 

int((-2*a__1+50)*(-a__1*x^2+10*a__1*x), x) = (-2*a__1+50)*(-(1/3)*x^3*a__1+5*a__1*x^2)int((-2*a__1+50)*(-a__1*x^2+10*a__1*x), x = 0 .. 10) = 0 = (500/3)*(-2*a__1+50)*a__1 = 0"(->)"[[a__1 = 0], [a__1 = 25]]

Answer considered is a__1 := 25 = 25``

   Hence the final solution is y(x) = a__1*x(10-x) = y(x) = 25*x(10-x)    ````

restartNULL

``

NULL


Download Approximate_solutions.mwApproximate_solutions.mw

@Kitonum 

Your suggestion is very useful and working when i call back eta later. Thanks. Ramakrishnan V

@Preben Alsholm 

Thanks

@tomleslie 

Dear maple user, I am sorry for not giving a staight answer to the question raised. The reason is i could not follow the f(x) in the equation and unclear Boundary Conditions. I am also not supposed to give a numerical answer to a question where the answer may be in functional form.

All I have given is a similar problem for which answer is available and can be used to understand more about alternative answers. Thanks for your comment and I confess that I have no full answer and  i know clearly that they two are not exactly same.

Probably I shall try to put the equation provided by tomleslie and get an answer with my boundary conditions.

Cheers.

heattransfer1Dfea_answer.mw

Example 1: Find the normal temperatures in a composite wall shown in Fig. The wall is maitained at 1000 C at the left face and convection mode of heat transfer occurs between the right face and surrounding fluid. Thermal conduction  

 

restart:

K__1 := 0.6e-1:W/cm 0C;

Convection heat transfer coefficient between wall and fluid, h := .1: W/m2 0C; Field variable is temperature phi; Temperature at node 1, phi__1 := 100: 0C;

Temperature of the fluid, phi__f := 25:  0C:       

Area of cross section, A := 1:``  cm2 ; Length of element 1, l__1 := 2:cm:  Length of element 2, l__2 := 4:cm;

``

Solution:``

 

Finite Element Equation for elements

For element 1 (nodes 1-2), the finite element equation is

For element 1 (nodes 2-3), the finite element equation is

                           Vector[row](2, {(1) = 1, (2) = 2})

(Vector(2, {(1) = F__1, (2) = F__2})) = K__1*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))/l__1.(Vector(2, {(1) = u__1, (2) = u__2}))     Vector(2, {(1) = 1, (2) = 2}) 

``

                        Vector[row](2, {(1) = 2, (2) = 3})

(Vector(2, {(1) = F__2, (2) = F__3})) = K__2*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))*(Vector(2, {(1) = u__2, (2) = u__3}))/l__2    Vector(2, {(1) = 2, (2) = 3})

``

 

Stiffness Matrices for elements

Stiffness matrix for element 1

Stiffness matrix for element 2

K__1*A/l__1 = 0.3000000000e-1

M1 := K__1*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))/l__1 = Matrix([[0.300000000000000e-1, -0.300000000000000e-1], [-0.300000000000000e-1, 0.300000000000000e-1]])

``

K__2*A/l__2 = 0.5000000000e-1

for conduction only

M2__1 := K__2*A*(Matrix(2, 2, {(1, 1) = 1, (1, 2) = -1, (2, 1) = -1, (2, 2) = 1}))/l__2 = Matrix([[0.500000000000000e-1, -0.500000000000000e-1], [-0.500000000000000e-1, 0.500000000000000e-1]])

for convection only

h*A = .1NULL

M2__2 := h*A*(Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})) = Matrix([[0., 0.], [0., .100000000000000]])NULL

 for both conduction and convection combined,

M2 := M2__1+M2__2 = Matrix([[0.500000000000000e-1, -0.500000000000000e-1], [-0.500000000000000e-1, .150000000000000]])NULL

``

NULL

 

Combining the matrices we get global stiffness matrix,

 

Kglobal := Matrix(3, 3, {(1, 1) = M1(1), (1, 2) = M1(3), (1, 3) = 0, (2, 1) = M1(2), (2, 2) = M1(4)+M2(1), (2, 3) = M2(3), (3, 1) = 0, (3, 2) = M2(2), (3, 3) = M2(4)}) = Matrix([[M1(1), M1(3), 0], [M1(2), M1(4)+M2(1), M2(3)], [0, M2(2), M2(4)]])``

``

Nodal Load Matrices for the elements

Load matrix for element1 (nodes 1-2) is

Load matrix for element2 (nodes 2-3) is

                         Nodes

L1 := Vector(2, {(1) = 0, (2) = 0}) = Vector[column]([[0], [0]])   Vector(2, {(1) = 1, (2) = 2})

``

 

h*phi__f*A = 2.5NULL

NULL

                                 Nodes

L2 := h*`&phi;__f`*A*(Vector(2, {(1) = 0, (2) = 1})) = Vector[column]([[0.], [2.50000000000000]])   Vector(2, {(1) = 2, (2) = 3})

Combining the element load matrics, we get global load matrix as

NULL 

F := Vector(3, {(1) = L1(1), (2) = L1(2)+L2(1), (3) = L2(2)}) = Vector[column]([[L1(1)], [L1(2)+L2(1)], [L2(2)]])NULL

``

 

      

Nodal Displacement Matrices for elements

Nodal Displacement matrix is

phi := Vector(3, {(1) = `&phi;__1`, (2) = `&phi;__2`, (3) = `&phi;__3`}) = Vector[column]([[phi__1], [phi__2], [phi__3]])

``

 

  Finite element equation for the system is [Kglobal.[ϕ]=[F];``

  Matrix([[M1(1), M1(3), 0], [M1(2), M1(4)+M2(1), M2(3)], [0, M2(2), M2(4)]])     =      

 

Substituting the nonzweo boundary condition phi__1 = 100 , deleting  the row and column containing `&phi;__1` and transferring the column containing node to rhs, we get 

``

(Matrix(2, 2, {(1, 1) = m[1, 1], (1, 2) = m[1, 2], (2, 1) = m[2, 1], (2, 2) = m[2, 2]})).(Vector(2, {(1) = `&phi;__2`, (2) = `&phi;__3`})) = (Vector(2, {(1) = m[1, 1], (2) = m[2, 1]}))+`&phi;__1`*(Vector(2, {(1) = -M1(2), (2) = 0})) = (Vector(2, {(1) = 0.800e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`-0.500e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`, (2) = -0.500e-1*`#msub(mi("&phi;",fontstyle = "normal"),mi("2"))`+.1500*`#msub(mi("&phi;",fontstyle = "normal"),mi("3"))`})) = (Vector(2, {(1) = HFloat(3.0), (2) = HFloat(2.5)}))

G := Matrix(2, 3, {(1, 1) = m[1, 1], (1, 2) = m[1, 2], (1, 3) = m[1, 3], (2, 1) = m[2, 1], (2, 2) = m[2, 2], (2, 3) = m[2, 3]}) = Matrix([[0.8e-1, -0.5e-1, 3.0], [-0.5e-1, .15, 2.5]])````

"By Gaussian Elimination,"

 

Matrix([[0.8e-1, -0.5e-1, 3.0], [0., .1187500000, 4.375000000]])

(1)

`` solve({.1188*phi__3 = 4.375, 0.8e-1*phi__2-0.5e-1*phi__3 = 3}) = {phi__2 = 60.51662458, phi__3 = 36.82659933}``

``

NULL


Download heattransfer1Dfea_answer.mw

@tomleslie 

Thank you for prompt response. I give the full document here. In case you could find time, please do the needful to help me learn more abou the intricacies in using these command. Thanks in advance.

First 7 8 9 10 11 Page 9 of 11