En este trabajo se introducen, en el contexto del Método de Elementos Finitos, dos alternativas posibles en relación con el concepto de acción repartida equivalente. La primera consiste en emplear pocos elementos, elevando el orden de dicha acción, mientras que la segunda se basa en emplear un mayor número de elementos dejando la acción en el orden más bajo posible. Se ilustran ambas situaciones mediante aplicaciones a los modelos de vigas de Timoshenko y Bernoulli-Euler, empleando estas acciones con diferentes órdenes, las cuales aproximan a la acción original, mediante polinomios ortogonales de Legendre en cada elemento. Como conclusión destacable, se indica que cuando se considera el menor número posible de elementos, es decir uno, para los casos de carga poco regular, ha bastado con utilizar acciones repartidas equivalentes de orden ligeramente superior al mínimo (orden cuatro), para obtener una excelente aproximación en los desplazamientos, giros y esfuerzos en el interior de los elementos.
ABSTRACT
In the context of the Finite Element Method, two possible alternatives dealing with the concept of equivalent distributed load are presented in the paper. The first consist in using few finite elements, by slightly increasing the order of the load, while the second applies the use of a greater number of elements leaving the load in the lowest possible order. Both situations are sampled with application to the Timoshenko and Bernoulli-Euler beam models, with different orders of load are used. These equivalent distributed loads are the result of applying Legendre orthogonal polynomial approximations, to the original load, in each element. The most noteworthy conclusion is that when the least possible number of finite elements is used (i.e., one) also for considering low level of regularity load cases only equivalent distributed loads of slightly higher than minimum order (four) were needed to obtain an excellent approximation when computing the deflections, rotations, bending moments and shear forces inside the elements.
Palabras claveElementos finitosacción repartida equivalentesolución nodal exactasplines generalizadospolinomios ortogonales de Legendremodelos de vigas de Timoshenko y Bernoulli-EulerKeywordsFinite elementsequivalent distributed loadsexact nodal solutiongeneralized splinesLegendre orthogonal polynomialsTimoshenko and Bernoulli-Euler beam modelsINTRODUCCIÓN
El tratamiento de modelos unidimensionales por elementos finitos, permite obtener, cuando se emplean determinadas funciones de forma, resultados exactos para las variables nodales y localizaciones óptimas para la aproximación de los esfuerzos en el interior de los elementos. En este sentido las referencias (1)(2)(3) y (4) son clásicas en la literatura sobre elementos finitos y en particular sobre la aplicación al modelo de Bernoulli-Euler (también denominado Euler-Bernoulli o Navier-Bernoulli).
Las cuestiones de la exactitud nodal y la idea de aproximación en el interior de los elementos, han sido tratadas en dos trabajos previos sobre splines generalizados y elementos finitos, (5) y (6), en los que se introdujo el concepto de acción repartida equivalente, con el propósito de optimizar los resultados derivados de la aplicación de los elementos finitos al modelo de Bernoulli-Euler.
En otras publicaciones posteriores, (7)(8) y (9), también se han realizado aplicaciones en la misma línea. En la primera al modelo de Timoshenko y en las otras dos al análisis de pilares con comportamiento lineal o no lineal del material, empleando el modelo de Bernoulli-Euler.
En este trabajo se desarrolla y concreta el concepto de acción repartida equivalente de cualquier orden, que aproxima la acción original en cada elemento finito mediante polinomios ortogonales de Legendre, con aplicación también a los modelos de Timoshenko y Bernoulli-Euler.
En relación con el concepto de exactitud nodal, que es básico para el desarrollo que se sigue y que fue comentado ampliamente en los trabajos (5) y (6), cabe destacar respecto al modelo de Timoshenko, la publicación (10) en la que se plantea la construcción de las funciones de interpolación mediante soluciones de la forma homogénea del sistema de Timoshenko. También en (11) se trata la cuestión de la exactitud nodal con aplicaciones a modelos de vigas. Pero quizás la referencia que puede considerarse pionera, en este tema, es la de P. Tong (1). No obstante, como ya se indicó en (5), la idea de fondo es bien conocida en problemas lineales de ingeniería estructural: Las cargas nodales equivalentes correspondientes a una acción repartida, son las opuestas de las cargas nodales de equilibrio o de contorno, relativas a la solución del problema con condiciones de tipo homogéneo en los nodos. Es decir, mediante el proceso de empotramiento en los nodos y aplicando el principio de superposición, por la linealidad del problema, se pasa de las cargas repartidas originales a cargas concentradas únicamente en los nodos. Se llega por tanto de forma intuitiva al hecho de que los movimientos de los nodos de la estructura con las cargas originales y los de la que tiene únicamente cargas concentradas (las correspondientes por el proceso indicado) en los nodos, son los mismos y por tanto dichos movimientos son exactos.
Desde el punto de vista matemático, tal y como se indicó en (6), la idea que lleva a la exactitud nodal se corresponde, aunque con una metodología distinta, con la del procedimiento de la función de Green, aplicado a un conjunto discreto de nodos. En este último aspecto, la metodología desarrollada en esta línea de trabajo guarda cierta relación con el método de Trefftz (12) y evoluciones posteriores del mismo en el campo de los elementos finitos (13), aunque el procedimiento expuesto en el trabajo se distingue por desarrollarse en dos etapas: primero se determinan los desplazamientos en los nodos y después se mejora la solución en el interior de los elementos.
Por otro lado, sobre el concepto de acción repartida equivalente tal y como se ha desarrollado y en conexión con la metodología de elementos finitos, los autores no han encontrado en la literatura aportaciones similares. Sin embargo ideas próximas, en el campo de las ecuaciones diferenciales sí las hay, pues los llamados métodos espectrales constituyen un ejemplo, en alguna medida, de lo indicado. En este sentido, en (14) se describen diversas técnicas de resolución de problemas de valor inicial y de contorno, consistentes en proyectar ortogonalmente los términos independientes de las ecuaciones sobre espacios de polinomios.
En trabajos previos se ha formulado dicho concepto, sobre la idea de que todos los problemas cuyas acciones originales generasen las mismas cargas nodales equivalentes en los elementos, tendrían la misma solución en los nodos. Al mismo tiempo se desarrollaron algoritmos basados en la interpolación, para obtener las soluciones en el interior de los elementos y que éstas se correspondieran con cargas de estructura sencilla. La acción repartida que allí se definió consistía en proyectar en cada elemento, la acción original, sobre el espacio de funciones formado por las soluciones de la ecuación diferencial, o sistema, del problema homogéneo. Por tanto, la acción repartida considerada era una función formada por una combinación lineal de soluciones del problema homogéneo y el orden mínimo de la misma, de acuerdo con la definición de acción repartida, que se introdujo en dichos trabajos previos, coincidía con el orden del operador diferencial de la ecuación. Por tanto, el orden de las acciones consideradas en las aplicaciones anteriores, por ejemplo a los problemas de Timoshenko y Bernoulli-Euler, era cuatro.
La aportación que aquí se realiza consiste en considerar otras acciones repartidas equivalentes de cualquier orden, con el propósito de llevar la aplicación de la idea a casos en los que se decida emplear un número muy reducido de elementos finitos, por ejemplo un solo elemento. Para ello se mejora la aproximación de dicha acción repartida, a la acción original, mediante desarrollos limitados de Fourier de la carga, en términos de los polinomios ortogonales de Legendre.
Estas dos cuestiones de considerar más elementos por un lado, manteniendo el orden de la acción repartida equivalente en el valor más bajo, y por otro, reducir el número de elementos, aumentando el orden de dicha acción en cada elemento, guarda cierto paralelismo con la versión h-p del método de elementos finitos tratada por I. Babuska y otros, en (15) y en otros trabajos previos, (16) y (17). Este enfoque para el modelo de Timoshenko se ha seguido en (18) empleando polinomios del mismo grado para aproximar el desplazamiento y el giro, pero con un enfoque diferente al utilizado en este trabajo, al no considerar dicho autor la exactitud nodal, que evita además el fenómeno de bloqueo en la solución, ni emplear tampoco la proyección de la acción original sobre el espacio de soluciones del problema homogéneo.
El interés por el estudio de estos modelos de Timoshenko, Bernoulli-Euler y otros próximos, como los de alto orden, y su tratamiento mediante diferentes técnicas numéricas, esencialmente la de elementos finitos, sigue vigente, como se pone de manifiesto por el gran número de trabajos que continúan apareciendo en la literatura. Algunas publicaciones recientes permiten observar diferentes tendencias. En (19) se emplean soluciones fundamentales para la viga de Timoshenko y por tanto se obtienen soluciones exactas pero de forma laboriosa y diferente a la metodología expuesta en este trabajo. En (20) los modelos Bernoulli-Euler, Timoshenko y otros de alto orden como los de Levinson y Reddy, son reformulados mediante teorías no locales de elasticidad que permiten también un cálculo por elementos finitos. El análisis de vibraciones de brazos de robot, hélices de helicópteros, etc, es un problema de gran interés que se modeliza dentro de la teoría de vigas de Timoshenko para los casos de poca esbeltez, y con tratamiento numérico por elementos finitos esencialmente, aunque también se están introduciendo otras técnicas con algunas ventajas computacionales, como la de cuadratura diferencial (21). Asimismo continúan las formulaciones del modelo de Timoshenko y su resolución empleando variables ficticias en desplazamiento o giro, transformándolo en un modelo análogo de Bernoulli-Euler y tratamiento posterior por elementos finitos, como se indica en (22). Por último se puede señalar que aunque el modelo de Timoshenko es una referencia constante, quizás sea en los modelos de alto orden donde los investigadores están poniendo mayor énfasis actualmente (23), (24) y (25).
PRELIMINARES
Las ecuaciones diferenciales, vectorial y escalar, respectivamente, para los modelos de Timoshenko y Bernoulli-Euler, son L(U) = F, L(U) = [–(ksAG(w′ – ψ))′, – (EIψ′)′ – ksAG(w′ – ψ)]T, con U = [w, ψ]T, F = [f, 0]T y L(w) = f con L(w) = (EIw″)″, como se expone en (5)(6)(7) y (10), donde ks, A y G son, respectivamente, el factor de corrección para el esfuerzo cortante, el área de la sección y el módulo de rigidez de cortante. Llamando H = EI y K = ksAG el sistema de Timoshenko puede también ponerse en la forma
La formulación variacional se obtiene, realizando la ponderación usual en el sistema original, con V = [v φ]T y tras el correspondiente proceso de integración por partes en el intervalo genérico [α, β], resulta
donde las formas, bilineal y lineales, son respectivamente
representando, las cuatro últimas, las acciones nodales de equilibrio o de contorno, que se corresponden en el extremo derecho x = β con los valores del esfuerzo cortante y momento flector respectivamente, y con los valores opuestos de dichos esfuerzos, en el extremo izquierdo x = α.
A partir de la expresión variacional [2] para el intervalo Ω = [a, b], y por ejemplo, para condiciones de contorno esenciales de tipo homogéneo (pieza biempotrada), puede formularse el problema: Obtener (w, ψ) ∈ H10 (Ω) × H10 (Ω) tal que
donde H1 (Ω) es el espacio de Sobolev H1 (Ω) = {g / g, g′ ∈ L2 (Ω)} y H10 (Ω) el formado por los elementos del espacio anterior que se anulan en los extremos del dominio. Cuando las condiciones esenciales son no homogéneas el problema se formula en los correspondientes trasladados de los espacios vectoriales de las funciones de ponderación. Por otra parte, la formulación variacional [2] se puede poner integrando nuevamente por partes, como
desapareciendo el último sumando si se toma V en el espacio de soluciones del problema homogéneo, y en dicho caso no dependiendo el primer miembro de los valores de U en el interior del intervalo. Para dicha situación donde L(V) = [0 0]T, los métodos de Galerkin: Bunov-Galerkin y Petrov-Galerkin generan las mismas ecuaciones discretas tal y como se indica en (6) y (7).
La ecuación de equilibrio para el elemento genérico [α, β] siguiendo lo expuesto en las referencias citadas es: Keue = fe + qe donde los desplazamientos nodales son
las cargas nodales equivalentes
las nodales de equilibrio
Los elementos de la matriz de rigidez local Ke = [keij] son keij = a(N
i, Nj) = li (Nj) que se pueden obtener mediante derivación de las funciones de forma. También en la forma usual mediante integración:
Las funciones de forma Ni = [N1iN2i]T, i = 1,...,4 verifican que li(Nj) = δij, i, j = 1,...,4 al constituir la mismas una base de Lagrange para la interpolación.
La aproximación de los desplazamientos y giros en cada elemento, viene dada por las soluciones del sistema homogéneo de Timoshenko, expresadas mediante las funciones de forma en función de los movimientos en los extremos por [10]
Las expresiones de N1i, N2i, i = 1,...,4 se incluyen en el Apéndice A para los casos donde H y K son constantes en cada elemento. En (7) se dan también para los casos de rigidez variable. La matriz de rigidez Ke deducida en (7) por derivación es
donde m = H / K = EI / (ks AG), h = β – α, λ = 1 + 3m / h2, γ = 1 – 6m / h2, expresada con igual notación que en (10) pero deducida de forma diferente.
Las ecuaciones de equilibrio local Keue = fe + qe y global KGuG = FG + QG, determinada ésta última mediante el correspondiente ensamblado de las ecuaciones locales, son exactas en el sentido de que son verificadas por cualquier solución del sistema L(U) = F tal y como se expone en (7).
Para cualquier problema, con el número necesario de coacciones (que impiden que sea un mecanismo), una vez considerados los desplazamientos y acciones conocidas (a partir de las condiciones esenciales y naturales dadas), resolviendo el correspondiente sistema (el reducido), se determinan los desplazamientos y giros desconocidos y después las acciones nodales de equilibrio globales desconocidas. Entrando luego en cada ecuación de equilibrio local se determina para cada elemento, el vector qe = Ke ue – fe, es decir, los valores exactos de los esfuerzos. Finalmente, de estos últimos los valores exactos de w′ y ψ′ en α+ y β– para cada elemento [α, β]. También de la relación K(w′ – ψ) = –(H ψ′)′, al ser conocidos w′ y ψ, quedan determinados asimismo los valores ψ″(α+) y ψ″(β–).
Observación: Para m = 0 y tomando ψ = w′ puede considerarse en el desarrollo, el caso de Bernoulli-Euler dentro del de Timoshenko, con las cargas nodales de equilibrio, como se expone en (5), dadas por
ACCIÓN REPARTIDA EQUIVALENTE. PROPIEDADES Y CONCEPTO DE ORDEN
Un estado de carga definido por f en el dominio [a, b] del problema, viene dado en general como la acción conjunta de cargas repartidas definidas por funciones continuas a trozos, y deltas de Dirac y dipolos para acciones puntuales. Estas cuestiones se describen con amplitud en (5) y (6).
Acciones equivalentes: Considerando la descomposición del dominio en la forma , se dice que dos acciones f y f son equivalentes respecto de dicha descomposición, si en cada [α, β] de la descomposición, f y f, generan las mismas cargas nodales equivalentes. O sea, si se verifica la igualdad de los siguientes productos escalares
definidos en el espacio de funciones de cuadrado integrable Lebesgue (14) con función peso ŵ(x) = 1, donde
Si f y f son equivalentes, las cargas nodales equivalentes globales también son las mismas, y ambas generan los mismos desplazamientos y giros en los nodos. También ambas cargas, f y f, generan los mismos esfuerzos en los extremos de cada elemento. Obsérvese que una vez determinados los valores nodales, que son exactos, se tiene que
y la expresión análoga para las f equivalentes. Ambas determinan por tanto los mismos esfuerzos en los extremos de cada elemento y consecuentemente los mismos valores de w′ y ψ′, en α+ y β– en [α, β], como se ha indicado. Recíprocamente, dada la carga f y determinados los valores de w y ψ en α y β y asimismo w′ y ψ′ en α+ y β– en cada subintervalo, pueden determinarse las funciones w, ψ y f interpolando w y ψ los valores indicados, y verificando al tiempo, el sistema de ecuaciones diferenciales [1]. Como consecuencia la acción f hallada es equivalente a la original f.
Observación: En la ecuación [15], el vector de cargas nodales de equilibrio, también se puede expresar en la forma
Propiedades de ortogonalidad e interpolación: De las [13] teniendo en cuenta que
integrando por partes y considerando la igualdad de cortantes y momentos en los extremos de cada elemento, resulta [18]
es decir, la función diferencia entre el cortante para la acción original f y el cortante para la acción equivalente f, es ortogonal al espacio de funciones engendrado por las N'1i, i = 1,...,4. Se puede demostrar fácilmente que engendran un espacio de dimensión tres.
De manera análoga, la diferencia, entre el momento flector para f y el momento flector para la acción equivalente f, es ortogonal al espacio de funciones engendrado por las N′′1i, i = 1,...,4. Estas engendran un espacio de dimensión dos. Además los cortantes en los extremos para ambas acciones coinciden como se ha indicado, luego se tiene junto con la ortogonalidad citada, la propiedad de interpolación en el sentido de Lagrange
Para los momentos, se verifican junto con la ortogonalidad, los siguientes resultados de interpolación en el sentido de Hermite, pues Q = –M′
Estas propiedades de ortogonalidad e interpolación, en cada elemento, para los momentos flectores y esfuerzos cortantes, para una acción equivalente a la acción original, son la base del buen comportamiento del método expuesto, para la aproximación de los esfuerzos y también para los desplazamientos y giros.
Acción repartida equivalente: De las infinitas acciones equivalentes a una acción dada f se puede elegir una f con gran regularidad y muy sencilla de calcular así como los desplazamientos, giros y esfuerzos asociados a ella. Dicha acción fue denominada en (5), (6) y (7) como acción repartida equivalente y es la definida como la proyección ortogonal de f en el espacio de funciones engendrado por las N′1i, i = 1,...,4 en cada subintervalo [α, β] por
Expresión que resulta del sistema de ecuaciones normales Ĝ{λi} = {fi}, donde Ĝ = (ĝij) (para distinguirla del módulo de rigidez de cortante G = E / [2(1 + v)]) es la matriz de Gram relativa a la base constituida por las N1i, i = 1,...,4, la cual tiene por elementos los productos escalares, para el peso ŵ(x) = 1
Un hecho importante es que se puede calcular la solución U = [wψ]T y los esfuerzos derivados de ella Q(x), M(x) relativos a f, como se indica en el Apéndice B, sin necesidad de calcular previamente f.
Propiedad de solución exacta I: Cuando la acción original f, en el elemento [α, β], es del espacio engendrado por las N1i, i = 1,...,4, entonces la solución U = [wψ]Tcoincide con la exacta en el elemento.
Esta propiedad es utilizada para la determinación de la solución exacta en (5), (6) y (7) y también en este trabajo (donde se incluye además la propiedad II). Dicha solución exacta se ha determinado por esta vía para realizar las comparaciones con las soluciones aproximadas del nuevo procedimiento propuesto.
A continuación se introduce el concepto de orden de la acción repartida equivalente.
Orden de la acción repartida equivalente: Se define este concepto como la dimensión k ≥ 4 del espacio Wk en el que se toma la acción repartida equivalente f. Siempre puede elegirse un conjunto arbitrario de funciones linealmente independientes
ampliando el espacio engendrado por las funciones N1i, i = 1,...,4 de soluciones del problema homogéneo, de forma que φi, i = 1,...,4 generen el mismo espacio que las N1i. La base del espacio Wk con k ≥ 4 puede construirse, para facilitar el proceso algorítmico, de modo que sea una base ortogonal respecto del producto escalar definido en L2ŵ([α, β]) para ŵ(x) = 1
es decir, que f sea la proyección ortogonal de f sobre Wk.
Los coeficientes de Fourier de f respecto de la base ortogonal B son
La carga f se puede poner en la forma fk (haciendo referencia explícita al orden) donde
Luego (g, N1j)ŵ = 0, j = 1,...,4. Véase que la información relevante de fk, respecto a los esfuerzos, está en y no en g, al generar la primera las mismas cargas nodales equivalentes que la f original, pues la carga adicional da, por la construcción realizada, cargas nodales equivalentes nulas y por tanto esfuerzos nulos. Por lo anterior g es una carga equivalente a la nula.
Obsérvese que para k > 4 no es un elemento del espacio engendrado por las funciones N1i, i = 1,...,4, soluciones de la homogénea, mientras que , sí lo es.
En la práctica la aproximación a la solución se realiza, en la mayor parte de los casos, empleando la carga f4 cuando la acción f tiene cierta regularidad. Por otra parte, en otros casos menos frecuentes, donde la acción repartida tiene poca regularidad, si se desea mantener, al tiempo, un número reducido de elementos finitos, se emplearía entonces acciones repartidas equivalentes de orden creciente f4, f5, f6,... pertenecientes respectivamente a los espacios W4 ⊂ W5 ⊂ W6 ⊂ ... Las soluciones relativas a dichas acciones van mejorando progresivamente en el interior del intervalo, de manera acorde con la construcción sucesiva de los espacios Wk. Se tiene en cuenta para ello el resultado básico de ortogonalidad: teorema de Pitágoras para las funciones de L2ŵ([α, β]) con ŵ(x) = 1, dado por
Se deduce de lo anterior que el error mínimo cuadrático disminuye al hacer crecer k con tal de que las funciones φi que se van añadiendo, sucesivamente, no sean ortogonales a f, es decir, que (f, φi)ŵ ≠ 0. Como se expone en los resultados numéricos, en los casos más singulares donde la aproximación con f4 no es la deseada, se ha comprobado que con valores de k poco superiores a 4, por ejemplo k = 5, k = 6 o poco más, se consiguen ya excelentes aproximaciones a los valores de la solución exacta en el interior de los intervalos. Además si f ∈ Wn, n ≤ k entonces fk = f y se tiene:
Propiedad de solución exacta II: Cuando la acción original f, en el elemento [α, β], es de alguno de los espacios Wn, 4 ≤ n ≤ k entonces la solución Uk = [wkψk]Tcoincide con la exacta en dicho elemento. Cuando n = 4 resulta la propiedad de solución exacta I.
Se puede resumir, lo anterior indicando que, dada f, la carga repartida equivalente fk se descompone en la suma de la carga f4 que es la repartida equivalente a f y de orden mínimo, y de la carga g que es una carga repartida equivalente a la carga nula. Y que g mejora ligeramente los resultados dados por f4, y solamente en el interior de los elementos, pues la información significativa sobre movimientos y esfuerzos en los nodos y en el interior la da ya f4.
APLICACIÓN DE LOS POLINOMIOS ORTOGONALES DE LEGENDRE
De acuerdo con lo anteriormente expuesto, cuando H = EI y K = ksAG son constantes, los espacios Wk, para cualquier k > 4, engendrados por las funciones {φ1,φ2,...,φk}, se toman por simplicidad y por las ventajas que conlleva la elección, como Wk = Pk – 1, pues ya W4 = P3 al constituir los polinomios de grado menor o igual que tres, el espacio de soluciones del problema homogéneo (en relación con los desplazamientos w). Como a su vez se eligen, por sencillez algorítmica, las funciones de manera que sean mutuamente ortogonales respecto al producto escalar ya definido, resulta que con dicha elección las funciones son precisamente los elementos de la sucesión de polinomios ortogonales de Legendre {pi(x)}∞i = 0 en [α, β]. Esta sucesión, como es conocido, forma un sistema completo de funciones en el intervalo (26), es decir se verifica
para cualquier función f de L2ŵ([α, β]) con ŵ(x) = 1.
En resumen , i = 0,..., k – 1, siendo fk = f4 + g
En el desarrollo se toma, para el elemento genérico [α, β], la sucesión de polinomios ortogonales pn(x), n = 0,1,2,... que resultan al transformar, los polinomios de Legendre clásicos Pn(t) definidos en [–1,1], con el cambio de variable, t = 2[x –(α + β) / 2] / (β – α). Dichos polinomios de Legendre verifican la relación de recurrencia (14):
que permite generarlos con gran facilidad hasta el grado deseado. Hay múltiples formas para generar estos polinomios como se expone, por ejemplo, en (26), (27) y (28). El cambio citado mantiene la condición de estandarización de los de Legendre clásicos, pues pn(β) = 1 y pn(α) = (–1)n. Además se verifica (26)
Como pues las bases {p0, p1, p2, p3} y {N11, N12, N13, N14} generan el mismo espacio de los polinomios cúbicos W4 = P3, luego para dicha función es indiferente emplear una u otra base. Como se ha expuesto en (5), (6) y (7) no es necesario su cálculo para la determinación de la solución correspondiente, cuando se emplea el método interpolatorio expuesto también en el apartado 5 de forma esquemática.
Con la notación usual de teoría de la aproximación (14)fk = ∏ŵ, k – 1f ∈ Pk – 1, con ŵ(x) = 1, donde ∏ŵ, k – 1 es el operador proyección ortogonal sobre el espacio de polinomios Pk – 1.
Por otra parte y puede ponerse cada fn, n ≥ 5 en la forma
para k tan grande como sea necesario.
De lo anterior resulta que dada f con norma finita , es claro que para n → ∞. En la mayoría de las aplicaciones, es suficiente con emplear acciones repartidas equivalentes de orden relativamente bajo, como f4, f5, f6...
En el apartado 5 se indica el proceso algorítmico para la determinación de las funciones wk,ψk, Mk, Qk, correspondientes a la acción repartida fk. Dichas funciones se pueden determinar de manera progresiva a partir de las w4,ψ4, M4, Q4 correspondientes a la acción repartida f4 añadiendo sucesivamente las soluciones incrementales Δwn, Δψn, ΔMn = HΔψ'n, ΔQn = –HΔψ″n relativas a las cargas incrementales αnpn con n ≥ 4. Así para f5 = f4 + α4p4 se tiene
Para f6 = f5 + α5p5 de manera análoga
y así sucesivamente, pudiéndose poner también para
Propiedades osculadoras y de ortogonalidad para las soluciones incrementales relativas a los sumandos de la carga <italic><overline>g</overline></italic>:
La solución incremental en [α, β] es la correspondiente al polinomio de grado n, Δfn = αnpn con n ≥ 4. Dado que la carga αnpn es una función ortogonal a las N1i, i = 1,...,4, es decir al espacio P3 = W4, entonces genera cargas nodales equivalentes nulas. Como la correspondiente solución incremental Δwn, Δψn verifica que Δwn(α) = Δψn(α) = 0 y Δwn(β) = Δψn(β) = 0, la pieza sometida a dicha carga incremental, al tener movimientos nulos, equivale al caso biempotrado, y de la ecuación de equilibrio [15] se deduce que los cortantes y momentos son también nulos en los extremos. Obsérvese que la solución incremental tiene movimientos y esfuerzos no nulos en el interior de los elementos.
Ahora se deducen propiedades en cada elemento mediante la integración por partes.
La solución Δwn, Δψn del sistema de Timoshenko relativa a dicha carga Δfn = αnpn verifica:
de donde se deduce que Δψn ∈ Pn + 3, Δwn ∈ Pn + 4. La acción Δfn = αnpn como se ha indicado para g es también una carga repartida equivalente a la carga nula, pues se cumple la igualdad siguiente de productos escalares
Pero también, por la propiedad de los polinomios ortogonales, de ser ortogonal a cualquier polinomio de grado inferior (14), resulta
Además de tenerse las propiedades osculadoras o de interpolación deducidas de la ecuación de equilibrio [15], para los giros, con un contacto de la gráfica de Δψn con el eje x de orden 2 en los extremos del intervalo, es decir
se tiene del proceso de integración por partes de [32] que los siguientes productos escalares son nulos
Dado que H es constante, las relaciones anteriores conducen a que Δψ′″n es ortogonal, en el dominio [α, β], al espacio de polinomios Pn – 1, del mismo modo Δψ″n a Pn – 2, Δψ′n a Pn – 3 y finalmente Δψn es ortogonal a Pn – 4. Si, por ejemplo n = 4, entonces la solución en giros Δψn ∈ Pn + 3 = P7 y las funciones Δψn, Δψ′n, Δψ″n, Δψ′″n , o más ampliamente
tienen integral nula en [α, β], es decir, compensan áreas en dicho intervalo y las funciones Δψn, Δψ′n, Δψ″n y Δψ′″n se comportan, en este sentido, como los polinomios ortogonales de grados 1, 2, 3 y 4, no olvidando, que a su vez, se cumplen las propiedades de interpolación mencionadas.
De KΔw′n = KΔψn – (HΔψ′n)′ se deducen de [15], para los desplazamientos incrementales Δwn, propiedades osculadoras análogas, con un contacto con el eje x de orden 1 en los extremos del intervalo, pues
Por otra parte de e integrando por partes queda
En resumen Δw′n es ortogonal en el intervalo a Pn – 4 y Δwn es ortogonal a Pn – 5. Si por ejemplo, n = 4, entonces Δw′n tiene integral nula en [α, β] y si n = 5, las funciones Δwn, Δw′n, Δw′nx tienen integral nula en el intervalo y sus gráficas por tanto compensan áreas en el citado intervalo.
Para la viga de Bernoulli-Euler se deducen propiedades análogas, obsérvese que en este caso ψ = w′ y la ecuación (HΔw″n)″ = αnpn se puede expresar en la forma de Timoshenko (HΔψ′n)″ = αnpn.
Estas propiedades de interpolación u osculadoras y de ortogonalidad para los sumandos de g se ilustran gráficamente en el apartado de resultados numéricos. En ellos puede verse que la contribución de dichos sumandos a la solución es en general muy pequeña en relación con la correspondiente a la carga f4. De ahí que en la mayoría de los casos, como se viene indicando, con dicha carga se obtenga una buena aproximación a la solución exacta, empleando incluso, muy pocos elementos.
DETERMINACIÓN DE LOS MOVIMIENTOS Y ESFUERZOS PARA UNA ACCIÓN REPARTIDA EQUIVALENTE DE CUALQUIER ORDEN <italic>k</italic> ≥ 4
La resolución del problema de contorno definido por el sistema Timoshenko, con condiciones de contorno esenciales y naturales dadas en los nodos relativos a la partición del dominio [a,b], se realiza dando los pasos que se indican para las dos vías posibles: interpolación, para la acción f4 y proyección ortogonal, para la acción fk, con k ≥ 4.
Movimientos y esfuerzos para una acción repartida equivalente de orden 4. Método de interpolación
1) Se determinará, empleando el método de elementos finitos, el vector de desplazamientos nodales globales. 2) Para cada uno de los elementos, teniendo en cuenta la ecuación de equilibrio local dada por [15] y siendo conocidos del paso anterior los desplazamientos nodales (desplazamientos verticales y giros), se calculará el vector de cargas nodales de equilibrio. 3) De dichas cargas nodales de equilibrio y de los desplazamientos nodales ya determinados en el paso primero, se hallarán los valores de las derivadas primeras de los desplazamientos y giros en los extremos de cada elemento, y mediante la expresión interpolatoria [B.3] del Apéndice B, se determinan los movimientos y esfuerzos en el interior de cada elemento para la acción repartida equivalente f4. Obsérvese que dicha carga se calcula, por este método, en caso necesario, al final del proceso.
A continuación se describe el procedimiento de proyección ortogonal para acciones repartidas equivalentes de cualquier orden, que incluye el caso de la carga f4, siendo por tanto una vía alternativa para el mismo.
Movimientos y esfuerzos para una acción repartida equivalente de cualquier orden <italic>k</italic> mayor o igual a 4. Método de proyección ortogonal
Se procede como en el caso anterior determinando el vector de desplazamientos nodales globales.
Para cada elemento [α,β] se proyecta la carga f sobre el espacio Wk = Pk – 1 aplicando [24], es decir para i = 0,..., k – 1 quedando así hallada, fk, acción repartida equivalente de orden k ≥ 4. Esto puede hacerse de manera incremental como se indica después.
Se resuelve ahora el sistema de Timoshenko [1] para la carga fk en cada elemento. Sean, w̃k(x) y ψ̃k(x), una solución arbitraria de Hψ′″ = fk, Kw′ = Kψ – Hψ″ dada por las expresiones
La solución que en los extremos del intervalo toma los valores w(α), ψ(α), w(β), ψ(β), es
Obsérvese que se está sumando a la solución de la homogénea una particular de la completa con valores nulos en los extremos.
Los esfuerzos en el interior del intervalo quedan determinados por Mk = Hψ′k, Qk = –HΔψ″k o igualmente por
Resumiendo, dada f y fijado el orden k ≥ 4 de la acción repartida equivalente, quedan determinadas en cada intervalo [α,β], por el procedimiento indicado, las siguientes funciones:
que aproximan, respectivamente, los desplazamientos, giros, momentos flectores, esfuerzos cortantes y carga repartida, de la solución exacta del problema de Timoshenko (también de Bernoulli-Euler, con m = 0 y w′ = ψ).
Procedimiento incremental
En la práctica pueden hallarse las funciones wk(x), ψk(x), Mk(x), Qk(x), fk(x) de manera progresiva a partir de w4(x), ψ4(x), M4(x), Q4(x), f4(x), pues fn + 1 = fn + Δfn con Δfn = αnpn y
siendo
No es necesaria la corrección de valores en los extremos como en [40] si para el cálculo de las primitivas se emplea [44]
pues para n ≥ 4, resulta Δwn(α) = Δψn(α) = 0 por construcción y Δwn(β) = Δψn(β) = 0 por la ortogonalidad de pn con los polinomios de grado inferior a n.
RESULTADOS NUMÉRICOS
Se ha realizado una variedad de ejemplos, en la que se aborda el cálculo en régimen lineal de la deformada y de las leyes de esfuerzos de una viga de rigidez constante en toda su longitud y sección rectangular. Las dimensiones de la pieza son 9m de largo y sección de 0.2m de ancho y 1m de canto. Se han considerado para el coeficiente de Poisson, módulo de elasticidad y factor de corrección de cortante los valores: v = 0.2, E = 3 × 107kN / m2 y ks = 5/6, respectivamente.
Los diferentes ejemplos parten de cargas repartidas en toda la longitud de la viga y posteriormente se incluyen cargas puntuales aplicadas en puntos interiores de la misma. Para la aplicación del método de cálculo se han empleado 1 y 2 elementos finitos, con objeto de ponernos en situaciones muy desfavorables a la hora de comparar los resultados obtenidos mediante la acción repartida equivalente con los de la solución exacta. Los ejemplos se realizan tanto para el modelo de Timoshenko como para el de Bernoulli-Euler, siendo el segundo un caso particular del primero cuando m = 0 con ψ = w′.
EJEMPLO 1. VIGA EMPOTRADA-APOYADA CON CARGA REPARTIDA Y PUNTUAL, EMPLEANDO UN ELEMENTO FINITO.
Se analiza en este caso una viga empotrada en el extremo izquierdo y apoyada en el derecho, con una distribución de carga como se muestra en la Figura 1.
Viga empotrada-apoyada.
En las Figuras 2 y 3 se puede apreciar que tanto para el desplazamiento como para el giro, los resultados obtenidos con las acciones equivalentes de orden 4 y 5 son prácticamente iguales a las de la solución exacta. En las Figuras 4 y 5 se puede apreciar que los errores máximos obtenidos para el momento y cortante son del 2 % y del 10 % respectivamente.
Comparación de desplazamientos: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de giros: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de momentos flectores: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de esfuerzos cortantes: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
En la Figura 6 se puede ver el comportamiento de las acciones repartidas equivalentes f4 y f5 y cómo estas recogen el efecto de la carga puntual.
Representación de la acción equivalente: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Estos últimos gráficos sobre las acciones repartidas equivalentes, tienen solamente un propósito ilustrativo en relación con la teoría desarrollada en este trabajo, tal y como se señala en el apartado 4.
EJEMPLO 2. VIGA EMPOTRADA-APOYADA CON CARGAS REPARTIDAS DE DISTINTO SIGNO, EMPLEANDO UNO Y DOS ELEMENTOS FINITOS
Se analiza una viga empotrada en el extremo izquierdo y apoyada en el otro extremo, sometida a una carga repartida con cambio de signo, con el propósito de aumentar el nivel de irregularidad, tal como se muestra en la Figura 7. Este ejemplo se realiza empleando uno y dos elementos.
Viga empotrada-apoyada.
Los resultados obtenidos para la acción equivalente de orden 4 y empleando un solo elemento no se comportan del mismo modo que en el ejemplo anterior al haber introducido una carga con menor regularidad, debido al cambio de signo. Sin embargo, al elevar el orden de la acción repartida equivalente los resultados obtenidos mejoran considerablemente.
En las Figuras 8 a 12, se representan los resultados obtenidos para las acciones repartidas equivalentes sucesivas desde la de orden 4 hasta la de orden 7, destacándose que las de orden 5 y 6 son iguales debido a que α5 = 0. Para los desplazamientos (Figura 8) el error máximo es del orden de 17 % para la acción repartida de orden 4, para la de orden 5 los errores son inferiores al 1,5 % y para la de orden 7 es inferior al 1 %. En la gráfica se puede ver además que empleando 2 elementos y manteniendo el orden de la acción repartida en 4 se obtienen mejores resultados, ya que los errores son inferiores al 1 %.
Comparación de desplazamientos: (a) Timoshenko y (b) Bernoulli-Euler.
Para los giros (Figura 9), el error máximo obtenido para las acciones equivalentes de orden 4, 5 y 7 y un solo elemento son del 12 %, 2,5 % y 1,5 % respectivamente. Empleando 2 elementos y acción repartida equivalente de orden 4, el error es inferior al 1,5 %. En el caso del momento (Figura 10), el error máximo para la acción repartida equivalente de orden 4 es del 12 % bajando al 4 % para la de orden 5 y al 3 % para la de orden 7 empleando un elemento, y del 3 % con 2 elementos y orden 4.
Comparación de giros: (a) Timoshenko y (b) Bernoulli-Euler.
Comparación de momentos flectores: (a) Timoshenko y (b) Bernoulli-Euler.
Finalmente para el caso del cortante (Figura 11), el error máximo para la acción de orden 4 es del 29 % y del 17 % para la de orden 5, decreciendo al 14 % para la de orden 7. Cuando se emplean 2 elementos y orden 4, el error máximo es del 15 %. Aunque cabe mencionar que este valor del error se da en el entorno del punto donde se tiene cambio de signo en la carga repartida y que en el resto de los puntos es inferior al 10 %. Como puede verse el decrecimiento del error es más lento para el cortante debido a la singularidad que introduce el cambio de signo de la carga repartida.
Comparación de esfuerzos cortantes: (a) Timoshenko y (b) Bernoulli-Euler.
En la Figura 12 se aprecia como las acciones repartidas equivalentes recogen el cambio de signo de la carga original.
Representación de la acción equivalente: (a) Timoshenko y (b) Bernoulli-Euler.
Adicionalmente en las Figuras 13, 14 y 15 se ilustran las propiedades osculadoras y de ortogonalidad para las soluciones incrementales relativas a los sumandos de g (solo para el modelo de Timoshenko). Asimismo se pueden comparar los valores de los desplazamientos y giros (Figura 13), así como de los momentos flectores y esfuerzos cortantes (Figura 14) para la solución exacta, solución equivalente de orden 4 y las soluciones incrementales de orden superior, empleando un solo elemento.
Propiedades osculadoras y de ortogonalidad (Timoshenko): (a) desplazamiento y (b) giro. Un elemento finito.
Propiedades osculadoras y de ortogonalidad (Timoshenko): (a) momento flector y (b) cortante. Un elemento finito.
Representación de las sucesivas cargas repartidas equivalentes (Timoshenko). Un elemento finito.
En las gráficas se puede ver claramente que los resultados provenientes de la acción repartida equivalente de orden 4 son los más significativos y la contribución de los sumandos de g es en general pequeña. Concretamente, se observa en este caso cómo a medida que aumenta el orden de la acción repartida equivalente las soluciones incrementales decrecen notablemente a partir de Δw4 para el desplazamiento y de manera similar, para el giro y el momento. Sin embargo, para el cortante el decrecimiento es más lento y empieza a ser significativo a partir de ΔQ6.
Al sumar a la w4, la incremental Δw4 resulta la w5, que coincide prácticamente con la solución exacta, tal y como puede observarse en la Figura 8. Se ha incluido la gráfica de la w7 = w4 + Δw4 + ... + Δw6 con propósito ilustrativo, y de manera análogo para el giro y los esfuerzos.
Téngase en cuenta, de acuerdo con lo expuesto en el apartado 4, que , con
EJEMPLO 3. VIGA BIEMPOTRADA CON CARGA PUNTUAL, EMPLEANDO UN ELEMENTO FINITO
En este ejemplo se analiza una viga biempotrada con una carga puntual de 150kN en el centro del vano, como se indica en la Figura 16. Se trata de ver un caso con mayor singularidad al considerar únicamente una acción puntual. En las Figuras 17, 18, 19 y 20, se pueden comparar los valores de los desplazamientos, giros, momentos flectores y esfuerzos cortantes de la solución exacta con los de las acciones repartidas equivalentes de distintos órdenes.
Viga biempotrada, carga puntual.
Comparación de desplazamientos: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de giros: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de momentos flectores: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Comparación de esfuerzos cortantes: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
Se puede ver que para el caso de los desplazamientos (Figura 17) el error máximo es del orden del 14 % para la acción repartida equivalente de orden 4, pasando al 3 % para la de orden 5 y al 1,6 % para la de orden 7. En los giros (Figura 18) el error máximo obtenido para la de orden 4 es del 8 %, del 2,5 % para la de orden 5 y del 1,7 % para la de orden 7. En el caso del momento (Figura 19) el error máximo para la acción equivalente de orden 4 es del 19 %, para la de orden 5 del 12 % y para la de orden 7 del 8,5 %. Como se puede ver el error disminuye lentamente por la singularidad de la carga, lo que se ha puesto de manifiesto empleando acciones repartidas equivalentes sucesivas hasta la de orden 11.
Finalmente para el caso del cortante (Figura 20) se produce un fenómeno análogo al de Gibbs en el punto donde se aplica la carga, ya que origina un salto en la ley de cortantes. En la Figura 21 donde se representa la acción repartida equivalente, se puede apreciar cómo las gráficas de las cargas repartidas equivalentes tratan de ir recogiendo información de la carga puntual, adoptando localmente una forma de campana y elevando su ordenada en el punto de aplicación al tiempo que se estrecha su base.
Representación de la acción equivalente: (a) Timoshenko y (b) Bernoulli-Euler. Un elemento finito.
CONCLUSIONES
En el trabajo se ha expuesto un procedimiento que optimiza los resultados obtenidos por el método de los elementos finitos, aplicados a los modelos de Timoshenko y Bernoulli-Euler, basado en el concepto de acción repartida equivalente de cualquier orden k ≥ 4. La acción original se ha aproximado mediante polinomios ortogonales de Legendre en cada elemento. Como conclusión destacable, se puede señalar que cuando se ha empleado el menor número posible de elementos, o sea uno, para los casos de carga poco regular, ha bastado con emplear acciones repartidas equivalentes de orden algo superior al mínimo (cuatro), para obtener una excelente aproximación en los desplazamientos, giros y esfuerzos en el interior de los elementos. Por otro lado para los casos de carga de mayor regularidad, con muy pocos elementos y acción repartida equivalente de orden mínimo, ha sido suficiente para conseguir una excelente aproximación en los resultados.
AGRADECIMIENTOS
Los autores expresan su reconocimiento a los revisores por las modificaciones sugeridas y el tiempo dedicado, que han permitido mejorar la configuración final del artículo. Asimismo, agradecen el apoyo recibido mediante la beca JAE-Predoc otorgada a E. Mercedes López (JAEPre098), correspondiente al Programa Junta para la Ampliación de Estudios del CSIC y cofinanciada por el Fondo Social Europeo.
APÉNDICE A
Las funciones de forma para la interpolación en los elementos finitos [α,β] en el modelo de viga de Timoshenko son N1i, N2i, i = 1,...,4 con z = x – α, m = H / K y h = β – α. Para m = 0 se transforman Ni, i = 1,...,4 que son las funciones de interpolación de Hermite para los elementos finitos de la viga de Bernoulli-Euler
siendo las N1i, N2i, i = 1,...,4
APÉNDICE B
De Hψ′″= f, w′ = ψ – mψ″ se deduce (7) que la solución equivalente U = (z = x – α) es
donde
Los desplazamientos, giros, flectores y cortantes equivalentes y asimismo la acción repartida equivalente, f en cada intervalo [α,β] recogidos en una única expresión son
siendo la matriz B la que se indica en (7) como A–1
REFERENCIASTongP.Exact solution of certain problems by finite-element methodAIAA Journal196971178180doi: http://dx.doi.org/10.2514/3.5067FilhoF. V.Comments on Computation of Streess Resultants from the Element Stiffness Matrices.AIAA Journal196863571572doi: http://dx.doi.org/10.2514/3.55382CraigA. P.JuneC. D.StricklinJ. A.Computation of Stress Resultants from the Element Stiffness MatricesAIAA Journal19664610951096doi: http://dx.doi.org/10.2514/3.3614BarlowJ.Optimal Stress Locations in Finite Element MethodInternational Journal for Numerical Methods in Engineering1976102243251doi: http://dx.doi.org/10.1002/nme.1620100202RomeroJ. L.OrtegaM. A.Acciones equivalentes y solución en desplazamientos interpolada en la viga de Benoulli-EulerInformes de la Construcción199849454527doi: http://dx.doi.org/10.3989/ic.1998.v49.i454.907RomeroJ. L.OrtegaM. A.Splines generalizados y solución nodal exacta en el método de elementos finitosInformes de la Construcción1999514644185doi: http://dx.doi.org/10.3989/ic.1999.v51.i464.872RomeroJ. L.OrtegaM. A.CorralesJ. M.Estudio y Resolución del Modelo de Viga de Timoshenko. Algoritmo de Acciones Equivalentes2002, 3-6 de junioV Congreso de Métodos Numéricos en IngenieríaMadridOrtegaM. A.Análisis del pandeo de pilares en regimen no lineal mediante splines generalizados2004ETS Ingenieros de Caminos – Universidad Politécnica de MadridMadridTesis DoctoralRomeroJ. L.OrtegaM. ANavarroF. J.Método de acciones equivalentes en el análisis del pandeo de pilares con comportamiento lineal o no lineal del materialVI Congreso de Métodos Numéricos en IngenieríaGranada2005, 4-7 de julioReddyJ. N.On locking-free shear deformable beam finite elementsComput. Methods Appl. Mech. Engrg.19971491-4113132doi: http://dx.doi.org/10.1016/S0045-7825(97)00075-3HlavacekI.KrizekM.On exact results in the finite element methodApplications of Mathematics2001466467478doi: http://dx.doi.org/10.1023/A:1013716729409HerreraI.Trefftz Method: A General theoryNumer. Meth. Partial Differ. Eq.2000166561580doi: http://dx.doi.org/10.1002/1098-2426(200011)16:6<561::AID-NUM4>3.3.CO;2-MQing-Hua QinTrefftz Finite Element Method and Its ApplicationsApplied Mechanics Reviews2005585318337doi: http://dx.doi.org/10.1115/1.1995716FunaroD.Polynomial Approximation of Differential Equations (Lecture Notes in Physics)1992Springer-VerlagBerlinGuoB.BabuškaI.The h-p version of finite element method, Part 1: The basic approximation resultComput. Mech.1986112141doi: http://dx.doi.org/10.1007/BF00298636BabuškaI.SzaboB. A.KatzI. N.The p-Version of Finite Element MethodSIAM J. Numer. Anal1981183515545doi: http://dx.doi.org/10.1137/0718033BabuškaI.SzaboB. A.Basic Mathematical ConceptsLecture notes on finite element analysis19832LikangL.Discretization of the Timoshenko Beam Problem by the p and the h-p Versions of the Finite Element MethodNumer. Math.1990571413420doi: http://dx.doi.org/10.1007/BF01386420AntesH.Fundamental solution and integral equations for Timoshenko beamsComputer and Structures2003816383396doi: http://dx.doi.org/10.1016/S0045-7949(02)00452-2ReddyJ. N.Nonlocal theories for bending, buckling and vibration of beamsInternational Journal of Engineering Science2007452-8288307doi: http://dx.doi.org/10.1016/j.ijengsci.2007.04.004FelixD. H.RossiR. E.BambillD. V.Análisis de la vibración libre de una viga Timoshenko escalonada, centrífugamente rigidizada, mediante el método de cuadratura diferencialRev, Int. Mét. Num. Cálc. Dis. Ing2009252111132FalsoneG.SettineriD.An Euler-Bernoulli-like finite element method for Timoshenko beamsMechanics Research Communications20113811216doi: http://dx.doi.org/10.1016/j.mechrescom.2010.10.009ChallamelN.Higher-order shear beam theories and enriched continuumMechanics Research Communications2011385388392doi: http://dx.doi.org/10.1016/j.mechrescom.2011.05.004GhugalY. M.SharmaR.A refined shear deformation theory for flexure of thick beamsLatin American Journal of Solids and Structures201182183195WangX. D.ShiG.Boundary Layer Solutions Induced by Displacement Boundary Conditions of Shear Deformable Beams and Accuracy Study of Several Higher-Order Beam TheoriesJournal of Engineering Mechanics20121381113881399doi: http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000440DavisP. J.Interpolation & Approximation1975Dover Publications, Inc.New YorkSzegoG.Orthogonal Polynomials1975American Mathematical SocietyProvidence, Rhode IslandFreudG.Orthogonal Polynomials1971Pergamon PressOxford